Contact map fingerprints of protein-ligand unbinding trajectories reveal mechanisms determining residence times computed from Scaled Molecular Dynamics
The binding kinetic properties of potential drugs may significantly influence their subsequent clinical efficacy. Predictions of these properties based on computer simulations provide a useful alternative to their expensive and time-demanding experimental counterparts, even at an early drug discovery stage. Herein, we perform Scaled Molecular Dynamics (ScaledMD) simulations on a set of 27 ligands of HSP90 belonging to more than 7 chemical series in order to estimate their relative residence time. We introduce two new techniques for the analysis and the classification of the simulated unbinding trajectories. The first technique, which helps in estimating the limits of the free energy well around the bound state and the second one, based on a new contact map fingerprint, allows the description and the comparison of the paths that lead to unbinding.Using these analyses, we find that ScaledMD's relative residence time generally enables the identification of the slowest unbinders. We propose an explanation for the underestimation of the residence times of a subset of compounds and we investigate how the biasing in ScaledMD can affect the mechanistic insights that can be gained from the simulations.
Introduction<!>Dataset overview<!>System preparation<!>Scaled molecular dynamics: theory<!>Scaled molecular dynamics: simulations<!>Analysis of the unbinding process<!>Early unbinding<!>Late unbinding<!>An example of the unbinding events analysis<!>Early unbinding: unbinding time and distance cutoff<!>Late unbinding: contact maps and exit paths<!>Clearly identifiable Transition States are crossed early<!>Most ligands follow nearly equally the two main exit pathways in ScaledMD simulations<!>Estimation of Relative Residence Time<!>Loop-in Helix<!>ScaledMD distinguishes short versus medium τ ligands within series and underestimates τ for ligands with long τ<!>Non-specific interactions at TS with unstructured α-helix 3 leads to τ underestimation<!>Conclusion
<p>During the past decade, different binding kinetic properties of potential drugs have emerged as one of the key factors that characterize their subsequent clinical efficacy. [1][2][3] It is thus advisable, if not necessary, to obtain good estimates of these properties even at the early drug discovery stage, for example when selecting chemical series found through virtual screening approaches. In this context, experiments to measure these properties cannot always be performed and in silico predictions of these properties provide a useful alternative. The interest of these kinetic properties in drug discovery, as well as the various methods that have been proposed to evaluate them, have been reviewed by several groups. [4][5][6][7][8][9][10] The molecular factors that are believed to influence binding kinetics include, among others, binding site accessibility, ligand and protein conformational fluctuations, as well as electrostatic and hydrophobic effects. 11 More recently, ligand desolvation has been pointed out as a key factor influencing residence time [12][13][14] (abbreviated as τ ). While several of these factors can be relatively straightforward to calculate, the actual estimation and evaluation of drug binding kinetics using computational methods remains challenging, indeed more so than drug binding thermodynamics, as it requires a sufficient sampling of the a priori unknown unbinding paths.</p><p>As members of the Kinetics for Drug Discovery (K4DD) consortium, 15 we were interested in the development and evaluation of molecular simulation approaches to compute drug-binding kinetic properties in the context of hit-to-lead or lead optimization projects. It led some of us to develop the τ -Random Acceleration Molecular Dynamics (τ -RAMD) 16 method for the evaluation of τ of a drug binding to its pharmacological target. It has been coupled with ML analysis 17 and MD-IFP (Molecular Dynamics Interaction Fingerprints) 18 for providing robust estimates of the residence time and insights into the features important for residence time using Machine Learning on the ligand's exit trajectories and applied to systems of varying levels of complexity. 19,20 Related less com-putationally intensive approaches that rely on the availability of a training set, such as COM-BINE analysis, 21 demonstrate remarkable performance by extracting protein-specific multilinear relationships between dissociation rate k of f and Lennard-Jones and Coulombic perresidue descriptors. Other approaches have the potential to rank and eventually prioritize new leads by estimating τ . Some of them, based on Metadynamics, 22 are able to accurately order 23 or even quantitatively estimate [24][25][26] τ , albeit only on small sets of compounds. These approaches rely on Collective Variables (CVs) whose number and definition are often system-specific and generally require trial and error iterations before being considered appropriate. 23 Nevertheless, several more general protocols have been proposed recently such as using generic electrostatics-like CVs, 27 Path CVs, 28 Machine Learning techniques methods 29 for finding optimal CVs and efficiently computing kinetics constants using them, 30 or Ratchet&Pawl simulations for finding CVs and reconstructing the full free energy landscape. 26 Recently, a method based on Steered Molecular Dynamics and on the evaluation of the desolvation energy during unbinding 31 has shown remarkable performance on a set of adenosine A 2A antagonists. One useful characteristic of this method is that it leads to simulation times that are approximately independent of the computed residence time of the ligand considered.</p><p>In another approach, the use of scaled potentials during Molecular Dynamics (MD) simulations, coined ScaledMD or smoothed potential MD, on multiple replicas of the system, combined with a statistical treatment for calculating the confidence in the estimations, has been shown to perform well for ranking τ in congeneric series on several targets in both a retrospective and a prospective manner. [32][33][34][35] The above-mentioned procedure has been integrated within the Biki Life Sciences suite 36 in the Biki-Netics package.</p><p>Our initial goal here was to evaluate the performance of Biki-Netics in quantitatively predicting τ and qualitatively ranking compounds within series and across series in the context of a given industrially-relevant drug discovery project. For this purpose, we performed a series of ScaledMD simulations on a set of 27 ligands belonging to several chemical series and binding to human HSP90 α N-terminal domain (N-HSP90) (Figure 1) in two different conformations. We report our results in light of those obtained with τ -Random Acceleration Molecular Dynamics (τ -RAMD) 16 using the same input structures, topologies and forcefield parameters, and of those obtained with non-equilibrium Targeted MD simulations 37 on a common subset of our dataset. The qualitative comparison between the computed relative residence time and the experiment shows that, within a chemical series, compounds are generally correctly ranked by ScaledMD. Quantitatively, two distinct trends emerge in the ScaledMD results: the τ estimation for the fastest unbinders (defined as τ τ lig3 , see below) is correct, but τ is almost systematically underestimated for the slowest unbinders (τ > 10 min). We analyze the unbinding events with two focuses: the early unbinding event, which is the crossing of the transition state (TS) before diffusing out of the protein, and the late unbinding event, where the ligand separates from the protein. Since no energy profile of the dissociation pathway was generated, we locate approximately the TS by correlating the computed τ and the time needed to reach a certain distance cutoff from the bound state. For describing the late unbinding event, we developed an ad hoc protein-ligand contact fingerprint (coined as cF P ) in which are listed the last residues close to the ligand during its exit. The clustering of the cF P highlights the two main exit pathways already reported in several other studies, 16,35,37 the front route (1) and the back route (2 in Figure 1), and we discuss their mechanistic relevance. We then dissect the early unbinding event and find that what is key for the exit time is the distance to the actual binding shell, defined as the residues in contact with the ligand in the bound complex. Through an in-depth examination of the unbinding trajectories in a molecular matched pair of compounds, 1 and 2, we relate the underestimation of τ to the interactions of the ligand at the TS with parts of the protein that are left unstructured by the scaling potential. We also compare several aspects of our work with results obtained by other teams 32,35 on close analogs of the compounds considered here. Overall, our results highlight the strengths and limitations of ScaledMD in the evaluation of protein-ligand binding kinetics and their mechanistic analysis in the perspective of prioritization of small molecules at the early stage of drug discovery process. Moreover, the protein-ligand contact fingerprint defined and used herein (cF P ) provides an efficient way for the mechanistic description of ligand unbinding process and a potential measure for direct comparison between different approaches.</p><!><p>A set of 27 inhibitors is considered in the present study (Figure 2). It is diverse in terms of chemical structure: the Tanimoto similarity of the ligands to their nearest neighbor in the set computed using Morgan fingerprint 38 is 0.46 ± 0.25. They display different bound protein-ligand conformations (both 'helix' and 'loop' conformations are represented, Figure S1) and a wide range of kinetic and thermodynamic properties, as summarized in the kinetic plot in Figure 2 (See the Supplementary Information -SI -for the structures, experimental data, and respective names in ref. 16 and 37 of the compounds considered in this work). The ligands can be clustered into 8 chemical series, seven based on their central heterocyclic moiety and one containing three singletons. Two of the ligands, compounds 7 and 21, can be clustered either in Cluster 1 or Cluster 3.</p><p>The two N-HSP90 conformations considered differ by the secondary structure of residues 105 to 113: these residues from the sequence subset C form the L1 loop that covers the binding site in the conformation represented in Figure 1 that will be called 'loop' in the text 39 (the sequence subsets A to G are defined in Figure 1). These residues are involved in the α 3 he- lix spanning residues 101 to 124 in the 'helix' conformation of L1 (see Figure S1 for a visual comparison of the two conformations). In four of the chemical clusters, all compounds in the series bind to the same N-HSP90 conformation; for the other three clusters, the N-HSP90 conformation is different depending on the compound. For all but two inhibitors (7, 22), Xray crystal structures were used as the starting point for system preparation. The starting structures of the complexes for the remaining two compounds were prepared as described in ref. 16.</p><!><p>The input structures and topologies were obtained from ref. 16 in Amber format and were then converted to Gromacs format using Acpype, 40 the Amber 11 utility ptraj and inhouse scripts. A detailed step-by-step description of the conversion process and all scripts can be found on the KBbox website. 41,42 All the ligands were protonated using Epik 43,44 at pH 7.5 and parametrized for the GAFF force-field using Antechamber. 45 RESP partial charges 46,47 were fitted from electrostatic potentials generated at the HF/6-31G*(1d) level using GAMESS. 48 The Amber14 force field 49 was used for the protein, a 10Å buffer of  TIP3P 50 water molecules was built around the complex with tLeap, 51 and Na + or Clions were added to ensure system neutrality and 150mM salt concentration to mimic physiological conditions. All protein-ligand systems were then minimized, heated, and relaxed for 10ns before the ScaledMD simulations.</p><!><p>In ScaledMD, the conformational space sampled in MD is increased by scaling down linearly the potential energy of the system by a factor λ.</p><p>For the dissociation of a protein-ligand complex with ∆G ± diss , ∆H ± diss and ∆S ± diss , the activation free energy, enthalpy and entropy of the dissociation process, we have:</p><p>Where A is the pre-exponential factor. The relative residence time τ comp between ligand l 2 and reference ligand l 1 can be expressed as: 33</p><p>∆S ± diss (λ) can be written as :</p><p>With ζ(l, λ) the deviation of the activation entropy from linearity with respect to λ.</p><p>As described in SI, under simple approximations τ comp can then be rewritten as</p><p>In cases where the smoothing of the Potential Energy Surface by the scaling leads a different ensemble of TS for the two ligands, the contributions ζ(l, λ) might be different so that e (ζ(l 2 ,λ)−ζ(l 1 ,λ))/R generally differ from 1.</p><!><p>In our study, each of the prepared proteinligand systems was used as the starting point of circa 20 independent MD simulations (with the exception of compound 16, for which only 7 simulations led to an unbinding event) using GROMACS 4.6.1 52 modified for ScaledMD 53 as implemented in the Biki suite (version 1.0.7). 36 The importance of using multiple replicas in MD simulations has been previously discussed in ref. 54 and thus a sufficiently large number of replicas was used. It is however important to note that some simulations did not led to a complete dissociation and ligand exit from the binding site and were discarded from the analysis.</p><p>To prevent protein unfolding during proteinligand dissociation, it is common practice in ScaledMD to restrain all protein atoms outside the binding site. In our systems, in order for the ScaledMD simulations to be comparable between different ligands, the same set of unrestrained protein residues was used throughout the simulations of all compounds and defined as follows. The protein residues closer than 6Å to all ligand heavy atoms in their bound structure were considered first. The selection was then expanded if needed in order to include single residues located between two selections. It led to the following residues being unrestrained : 47-56 (similar to subset B), 93-114 (similar to subset C), 132-138 (similar to subset D) and 183-186. Despite the fact that the initial 6Å selection criterion is the same as in ref. 32 (it is 5Å in ref. 35), in the present study, the final selection of unrestrained residues is larger, reflecting the greater diversity of the bound conformations in our dataset. The whole subset C is thus left free of restraints during the simulations. It allows for the protein to sample the same conformational space regardless of the ligand-bound starting conformation, and to represent the induced-fit known to occur for the helix binders, as subset C is in loop conformations when in apo form. Nevertheless, leaving subset C unrestrained translates into the sampling of several unfolded conformations of α-helix 3 and L1 during the ScaledMD simulations, which might have affected negatively the performance of the unbinding times estimations. Positional harmonic restraints with a force constant of 50 kJ/mol/nm 2 were imposed on the backbone of all restrained protein residues. As in ref. 32, the unbinding time t of f was defined as the first time when the ligand is surrounded by two full solvation shells, defined as a union of spheres of 6Å radius, centered on the ligand atoms and free of protein atoms. For each ligand, the unbinding time was computed from the average unbinding time over all replicas and its uncertainty was estimated by a 1000-fold bootstrapping (see Table S3 in SI).</p><p>One of the crucial parameters in ScaledMD simulations is the choice of the scaling factor λ. Similarly to ref. 32, different λ values were tested in order to identify the optimal one as the best balance between a reasonable computing time and unbinding times ranging from a few nanoseconds to tens of nanoseconds. For this evaluation, we simulated eight ligands (2, 3, 8-13) covering a large range of measured residence times (from four seconds to more than 2 hours), bound to two different conformations of N-HSP90 and belonging to different chemical series.</p><p>Three different values of λ were tested, namely 0.40, 0.42 and 0.45, leading to 480 simulations in total with t of f values ranging from 6ns to 300ns depending on the ligand and the scaling factor. A value of λ = 0.50 was also considered, but simulating the ligands with the smallest τ with this scaling factor led to mean t of f greater than 150ns. This is at least three times more than with λ = 0.45, suggesting that for the compounds with the largest τ in the set, the unbinding times would be in the 500ns-1µs range, which translates into a calculation time of 5 to 10 days for simulating the compounds with the largest τ on the computational facility we have used. This duration was considered to be too long for practical use in the context of an industrial virtual screening where tens to hundreds of compounds at the very least would need to be evaluated during a screening campaign. λ = 0.40 was found to be optimal for this set of molecules as the correlation coefficient between the computed and measured residence times was the best and the error was the smallest (See Figure S4 and text in SI). This value of λ was used for ScaledMD simulations for the 19 remaining ligands. All the related analyses can be found in SI. All simulations were performed on Intel processors under Linux on the Marconi facility at Cineca (Italy).</p><!><p>In order to characterize the features of the ligand unbinding trajectories, three complementary approaches were followed and implemented using MD Traj 55 and contact_map 56 python libraries. The first two approaches focus on characterizing the early stage of the unbinding process and the third describes the latest stages of the unbinding trajectories. The notations used in the text are in Table 1.</p><!><p>Before we describe the process followed for the characterization of the early unbinding, we need to define the metrics used. First, we define the initial contact shell S as the set of the protein residues that are close to the ligand l in the bound state (colored black in Figure 3). These are residues whose heavy atoms are found closer than 3.5Å to the heavy atoms of the ligand at the starting point of the simulation. Then, d(t) is the shortest distance between any heavy atom of the ligand and the nearest of the residues in S at time t. Finally, we define the center of mass (COM) metric d COM (t) as the distance between the COM of the ligand at time t and at t = 0, which corresponds to the ligand bound state. The above-mentioned metrics are illustrated in Figure 3.</p><p>For the characterization of the early unbinding using the first approach, for each ligand and for each exit trajectory, we perform the following steps: first, for each timeframe t of the whole trajectory, the distance d(i ∈ S, t) between the heavy atoms of the ligand and those of all residues of S is computed. Then, given a specific cutoff distance d c , the first passage time t f pt (d c ) is determined as the first time where</p><p>The second approach for characterizing the early unbinding relies on d COM (t), the displacement of the COM of each ligand l with respect to its initial position at time t along the trajectory, as described in Figure 3</p><!><p>As discussed in ref. 35, despite the challenge of characterizing the ligand exit pathways in ScaledMD simulations, these simulations might give indications on relevant structural aspects of their unbinding mechanisms. Herein, we propose a simple procedure for describing these exit pathways using as a metric a protein-ligand contact fingerprint (cF P ) which is built as described in Algorithm 1. As illustrated in Figure 3, this way of describing the late ligand unbinding process has two features: first, it ensures that the nearest distance between each of the nearest residues and the ligand during the last frames of the late unbinding is kept in the cF P , and second, by limiting the total number of residues selected in the process, it ensures that the subsequent comparison between cF P will not be skewed by the longest exit trajectories. Qualitatively, this fingerprint highlights the nearest residues encountered by the ligand during its unbinding from the protein. Importantly, this is done without the need for explicit consideration of time.</p><p>A distance function δ CM (cF P A , cF P B ) was devised for comparing the similarity of two unbinding trajectories A and B in contact map space. We define cF P A and cF P B as the contact map fingerprints of the two trajectories, N A and N B as the sets of the nearest residues, n(N A ) and n(N B ) are the cardinalities of the sets considered, so that the nonzero values of e.g. cF P A can be written as (res i,A , d i,A ) 1≤i≤n(N A ) : for example, in Figure 3, res 1,A = 31 and d 1,A = 4.44. The distance function δ CM (cF P A , cF P B ) is expressed by Eq. 5, using the union, intersection and symmetric difference (∆) between the two sets. This function simplifies to a normalized Euclidean distance when both trajectories have all their nearest residues in common, and d f ar is the contribution to δ CM for each of the nearest residues which are not in common in the two trajectories. A rare case happens when, for both trajectories, the ligand exits immediately after having reached d c and stays beyond d f ar for all protein residues. In this case, both n(N A ) = 0 and n(N B ) = 0 and by convention δ CM (cF</p><p>(Eq. 5) The function δ CM defined in Eq. 5 is used in order to perform cluster analysis of ensembles of late unbinding trajectories. First, δ CM is used to build a distance matrix between the cF P and a hierarchical clustering is performed using the complete method, with a threshold of 95% of the maximum distance between all clustered trajectories. Finally, the distance matrix and the cF P are sorted according to the similarity between the cF P determined at the clustering step.</p><!><p>As described in the Methods section, a set of analyses has been devised in order to best describe the early and late unbinding events and  3: a: Definition of the metrics used for the analysis of the early exit events. The residues of the initial contact shell S are indicated in black. The distance between the nearest heavy atoms of ligand and residues of S is indicated as d(t). The COM metric, d COM (t), measures the distance between the COM of the ligand at a given t and at t = 0. b, c: Illustration of the contact map fingerprint generation process. The fingerprint is built by scanning the exit event in reverse starting from the frame corresponding to ligand exit from the protein (defined here as frame 5, in blue) until the first frame after the first passage time t f pt (frame 2 in this example). Residues closer than d F P = 4.5Å relative to the ligand in each frame are added to the contact map fingerprint along with the nearest distance to the ligand among all frames.</p><p>to enable comparison between the different ligands. Below, we describe the outcome of the trajectory analysis for compound 8 using this set of tools. All the related results are depicted in Figure 4.</p><!><p>In Figure 4b    the initial bound state, and that the ligand is free to diffuse quickly up to the unbound state.</p><p>On the other hand, for a low enough value of d c such as 3Å, the ligand remains within the bound state energy well and the correlation between t of f and t f pt (d c ) is negligible. It is thus expected that the corresponding correlation coefficient increases sharply starting from a certain value of d c . In order to better characterize the limits of this basin, the correlation coefficient R(d, l) between t of f and t f pt (d c ) over all replicas was computed for several values of d c .</p><p>It was indeed found that for high enough values of distance d c , t f pt (d c ) and t of f were correlated and that the slope of this correlation R 2 (d c ) was close to 1. The same procedure was performed for relating the unbinding time t of f and the time t COM f pt (d c ) for the ligand COM to reach the distance d c from its initial position, which allows the computation of the squared correlation coefficient R 2 COM (d c ). The trends for all compounds of the set are displayed on Figure 5 and will be discussed later.</p><!><p>The cF P of the 20 unbinding trajectories with d c = 4.5Å is represented in Figure 4a. The cF P bits are related with the sequence subsets of HSP90 that are highlighted in Figure 1, and the cF P bits themselves are darker when the associated distance is higher. The distance matrix between the trajectories is displayed in Fig- ure 4d. It indicates the existence of three main exit clusters, as confirmed by the dendrogram in Figure 4c. Its coloring relates to the exit cluster with the unbinding and exit times displayed in Figure 4b. It shows for example that the earliest exit event belongs to the red trajectory cluster, which corresponds to direct exits of the ligand with only a few interactions outside subsets C and D residues. The cF P associated with the green cluster displays interactions with subsets A, C, D, E and F but not with subset B; it corresponds to the back route (route 2) in Figure 1, where the ligand goes below α-helix 3 (subset C) during its exit. In the case of the blue cluster, interactions with subsets B, C and D indicate the front route (route 1).</p><!><p>The evolution of R 2 , R 2 COM and of the mean RMSD of the protein (calculated for heavy atoms) with d(t) (resp. d COM (t)) is shown for all ligands in Figure 5. First, all mean RMSD plots have a very similar shape: they reach a plateau around d = 3.5Å to 4Å regardless of the bound conformation of the ligand. This plateau is lower for the compounds with the lowest τ comp as in this case the ligand leaves early, so that d rises quickly, before the unrestricted protein residues have fully sampled the conformational space available to them. Let us then consider the correlation coefficient R 2 , in green. It is close to zero when d is small : small displacements of the ligand out of the binding site do not end up in an exit event. However, in 18 of the 27 cases, R 2 then rises sharply to a plateau whose value is generally between 0.8 and 1.0. It means that beyond a certain distance threshold d ± , the ligand becomes free to diffuse rapidly out of the protein: the time to reach this threshold is highly correlated with t of f . In this case, a kinetically relevant TS is located before, but close to d ± . For 9 ligands, 7, 11, 14, 19, 20, 23, 24, 25, and 27, R 2 rises slowly and, in some cases, it stays well below 0.8, so that no specific TS event can be observed. Finally, the trend concerning R 2 COM , in blue, is less systematic, the plateau is lower than in the case of R 2 for the same ligand, and, in most cases, the plateau is reached for higher values of d. This last finding might come from the ligand keeping its interactions with some of its binding shell residues while it moves away from its bound position, either because these residues move along with the ligand, or because the ligand is large enough for its COM to shift significantly while some of its atoms are still close to the binding shell. It also demonstrates that residue-specific information describes more accurately the kinetically relevant unbinding event than the displacement of the ligand from its starting point, which is described by R 2 COM . Moreover, for several large compounds with rotatable bonds, such as 14, 18, 25, 26, and for some of the other com-pounds with the shortest τ values, such as 7, 11, 21, 27, the correlation coefficient R 2 COM is close to zero for the whole range of distances considered, either due to the limitation of the COM metric for the most complex compounds, or to a low or inexistent unbinding free energy barrier (see below for 11).</p><p>The comparison of these plots for two pairs of similar compounds highlights some consequences of the scaling on the description of the unbinding process. A first example is the pair of compounds 14 and 15. They share a common substructure and differ mainly by the quinoline moiety of 14 which interacts with the hydrophobic pocket formed under α-helix 3 (both are helix-binders). The RMSD plots indicate that this hydrophobic packing discussed in Kokh et al. 16 for these two compounds does not prevent destabilization of the protein during the unbinding process, as the RMSD plateau of compound 14 is similarly high and is reached sooner than for compound 15. The variation of R 2 with d for the two compounds is markedly different and this difference can be related to the differences in structure and illustrates a limit of the definition of the metric. Indeed, in the case of the long and buried compound 14, the ligand stays close to some residues of the initial binding shell until late in the unbinding process, which is much less the case for compound 15 and leads to an earlier increase of R 2 . A second matched pair is composed of compounds 7 and 21, with the methyl in 7 (loop binder) changed to a phenyl group in compound 21 (helix binder). They have a very similar K D , but compound 21 has only a slightly longer residence time.</p><p>Overall, the analysis shown in Figure 5 enables (for two thirds of the ligands) to pinpoint the threshold where the kinetically relevant TS has been crossed and after which the ligand is free to diffuse. This information can be of great value for identifying the interactions that most influence the kinetics of the unbinding process, i.e. the ones disappearing around the plateau of R 2 , especially in cases where the trends in R 2 and R 2 COM are different. The unbinding free energy profiles of compounds 9, 11 and 12 were computed in ref. 37 (SI Figure 2) using Ther-   modynamic Integration (TI). The results of this study and the present study are consistent: for compounds 9 and 12, d COM ± is found to be lying within the Transition State Region identified using TI, while in the case of 11, no free energy barrier was identified using TI along the unbinding pathway and no d COM ± can be identified, as R 2 COM stays well below 0.8 over the whole distance range considered. Finally, the influence of the potential scaling on d ± is illustrated in Figure S5 (top), where the squared correlation coefficients R 2 between the unbinding time t of f and the first passage time t f pt (d) is plotted against the distance d for the three values of λ used for the ligands of the training set. Two situations are observed: for ligands 3 and 8, d ± (λ = 0.45) > d ± (λ = 0.42), which indicate that the TS is displaced depending on the scaling. For the six other cases, d ± is similar for the different values of λ, suggesting a smaller effect of the scaling on the description of the TS.</p><!><p>The cF P of all unbinding trajectories ordered by cluster is shown in Figure S12. (Similar cF P and dendrograms representations together with the correponding distance matrices were produced for all ligands, see SI). 27 clusters were identified by the clustering analysis, the end of all exit trajectories were visually reviewed together cluster by cluster and annotated as Front-or Back-like trajectories at the cluster level. These visualizations confirm that unbinding trajectories belonging to the same cluster display similar features. In a few cases, the limitation of the cF P in describing topology makes it unable to assign the front or the back route e.g. where a ligand with a front-like cF P rotates and goes behind the unfolded α-helix 3 just when exiting.</p><p>In most trajectories, the ligand interacts with the beginning of subset C (bright green), which is the loop before α-helix 3 which surrounds the binding pocket. As discussed previously for ligand 8, other frequent combinations of interaction hotspots during ligand exit are subsets A, E and F on one hand (back route), and subsets B and D on the other (front route). It is interesting to note that for all ligands in the set, between c.a. 25% and 70% of all exit trajectories involve the back exit route (see individual heatmaps in SI). It is a greater proportion than in ref. 35 where the back route is found, but more anecdotally. This might be due to the fact that in our work the whole subset C is left free of restraints and adopts unfolded conformations during the ScaledMD simulations that open space for the ligands to exit by the back route. This further emphasizes the importance of choosing wisely the restraints used in ScaledMD simulations in order to balance between the ability to compare a large set of ligands using a larger set of unrestrained residues and the ability to obtain structurally relevant insights from the simulations by restricting the set of unrestrained residues. Indeed, in τ RAMD simulations, 16 the ligand exit channel that goes under the α-helix 3 is only open when α-helix 3 has a pure helical conformation and when the transient hydrophobic sub-pocket under this helix is open, which means that loop-binding compounds were not using the back route. Wolf et al. 37 have also compared those two routes using ligand 2 and found that the front route was favored in terms of work required to pull the ligand into the solvent. Nevertheless, the exit trajectories corresponding to the two groups display roughly the same mean exit time within statistical error in ScaledMD. This might be surprising at first, but it is less so considering the fact that the fingerprinting considers only residues that are beyond the d c = 4.5Å distance cutoff to the initial contact shell S. At this distance from S, as suggested by Figure 5, in most of the cases the ligand has already crossed the limits of the kinetic basin and the exact interactions with its surrounding do not have a notable effect on its exit time. From this analysis, it is tempting to question the relevance of studying the latest stages of the exit pathways when using ScaledMD for understanding the specifics of the Structure-Kinetics Relationships, as the interactions which have the most effect on exit time during simulations are those that are broken around the threshold in R 2 COM (d), which is often closer than d c to the initial contact shell residues.</p><!><p>Relative residence times, τ comp , were computed as the ratio of the mean unbinding time during the ScaledMD simulations to the mean unbinding time of ligand 3. τ comp were related to the measured residence times relative to ligand 3 τ rel (λ) = (t of f /t of f,lig3 ) λ and plotted in Figure 6.</p><p>The correlation obtained is modest when considering the whole set, but can be analyzed as the superposition of two trends. The first trend is related to the subset of compounds for which τ comp ∼ τ rel in Figure 6. The second trend is represented by a set of compounds, generally with larger relative residence time (τ rel (λ) > 2, i.e. t of f > 7.9.10 −2 s), where the slope between the relative residence time and τ comp is lower than 1. This trend is followed by larger and more rigid compounds such as 2, 9, 10, 14, the pair 15 and 16 from Cluster 7, 18, from Cluster 3 and the pair 23 and 24 from Cluster 5. In Schuetz et al., 35 compound 6, which is the largest of the dataset and has the longest measured residence time, was also underestimated in ScaledMD unbinding simulations performed with very similar settings. Similarly to their studies, this compound was the only compound whose τ is significantly underestimated in this dataset, it was rightfully considered as an outlier and discarded from analysis. As discussed in the Methods section, τ comp might be lower than one in cases where the smoothing of the Potential Energy Surface by the scaling favors a different ensemble of TS. A similar effect with respect to the scaling factor is illustrated in Figure S4, in SI, where the squared correlation coefficients R 2 (d) between the unbinding time t of f and the first passage time t f pt (d) is plotted</p><!><p>Figure 6: τ comp vs τ rel for λ = 0.4. Vertical error bars show standard deviation, and horizontal error bars show standard error. The color of the symbols corresponds to the chemical series (see Fig. 2) and their shape to the bound HSP90 conformation. The bold line indicates the line where τ comp and τ rel are equal and the dashed line shows the trend followed by a subset of the ligands with higher residence time.</p><p>against the distance d for the three values of λ used for the ligands of the training set. For compounds 3 and 8, the profiles are markedly different between λ = 0.45 and the two lower values, which suggests that the TS is located much further from the binding site with this scaling factor than with the lower ones.</p><p>In the next sections, we will describe and discuss the two trends in τ comp found from the analysis of several chemical series in this dataset: Clusters 1, 2, 3, 7 and 8 (see SI for description of Clusters 4-6). Based on the analysis of the pair from Cluster 8, we will propose an interpretation of the systematic underestimation of τ in one of the two subsets.</p><!><p>Figure S6 summarises the structures as well as the τ rel and τ comp values for the compounds of Cluster 1 and 3. Cluster 1 (in light blue) contains five compounds which have short residence times. They are close (structurally and in terms of kinetic constants) to the set of four HSP90 resorcinol ligands studied with the same methodology in ref. 32. In our case, however, the relative computed residence times are not linearly related to the relative measured unbinding constants. Their relative computed residence times are all very close and they do not allow ranking the compounds when the error estimated by bootstrapping is taken into account. This difference might be due to the fact that the compounds in this set are more diverse structurally than in reference 32. This structural diversity in the ligand structures also leads to interactions with a diverse set of residues in the bound state, leading to a larger set of unconstrained residues used in our work. Three of the Cluster 1 compounds (5, 8 and 11, see also Table S2) have also been simulated by Wolf et al. 37 with a very similar outcome: the three ligands are among the fastest unbinders but they are not differentiated quantitatively. Cluster 3 comprises seven compounds (including compounds 7 and 21, which are shared with Cluster 1) coming from a rescaffolding of the compounds from Cluster 1 and 2. More specifically, the resorcinol ring has been fused with the central heterocyclic ring to form an hydroxybenzopyrazole, which is modified on two different positions. In this series, the group in R2 position binds to the same hydrophobic pocket as the diverse groups in Cluster 2. On the other hand, the distal part of the group in R1 points towards the outside of the binding site, so that it does not interact much with N-HSP90 in the bound state. One large subseries of five compounds (17-20 and 22) is of N-methyl-amides bound to a diverse set of substituted phenyl or benzyl groups in R1. For 17, 19 and 20, the kinetic estimation is roughly correct, while it is underestimated for compounds 18 and 22 which are N-phenyl substituted compounds very close structurally to 17. These series are also close structural analogs of the four compounds in the main series studied in Schuetz et al., 35 which displays a similar trend, where the ligand with a long measured residence time is underestimated.</p><p>The two last compounds (7 and 21) are common between Cluster 1 and Cluster 3. They have the lowest τ and K D within Cluster 3, probably because they neither accept the hydrogen bond formed with the Thr184 side-chain that is in the large Cluster 3 subseries, nor donate an hydrogen bond to the Gly97 backbone C=O bond, as for 5 and 8 from Cluster 1, which are better binders than 7 and 21 (see Figure 7). They differ structurally only by a phenyl group in the R2 position. This structural modification leads to a modest difference in τ which is not captured by the simulations. Taken together, Clusters 1 and 3, which are close struc- turally, illustrate a general observation: neither the shortest τ compounds (with a measured relative residence time of less than 0.7) nor the longest τ ones (greater than than 0.7) are ranked correctly within each group. However, the longest τ compounds are clearly separated from the shortest τ compounds in Figure S6.</p><p>The relative residence times for compounds from Cluster 2 and 7 are displayed in Figure S7. In this case, the computed residence time for the shortest τ compounds (those with measured relative τ less than two) and for the longest τ ones (the rest of the set) is similar, but compounds can be ranked quite well within each set. Cluster 2 contains six compounds that are structurally very close as they represent single point variations, from R1=H (compound 3) and bromine (compound 6) to bulky amides and sulfonamides. For the latter compounds, the different R1 groups are located in the hydrophobic pocket beneath helix 3. 3 and 6 bind to N-HSP90 in its loop conformation, they are differentiated by the ScaledMD simulations and follow the expected quantitative trend, but τ comp for 6 is overestimated with respect to 3, which is used as a reference for the plot. It might be that the force-field parametrization of the bromine atom of 3 in GAFF, which does not include its polarization, leads to an underesti-mation of τ comp with respect to 6. Compounds 12 and 13, which bind to the helix conformation of N-HSP90, differ only by the position of a methyl group on a sulfonamide and are nearly indistinguishable by their relative τ comp . Compounds 9 and 10 are correctly differentiated, but their τ comp are underestimated with respect to τ rel . Interestingly, the τ of 9 is also underestimated by τ RAMD. 16 The two compounds in Cluster 7 only differ by the reduction of the ketone 16 to the alcohol 15 and their bound structures are superimposable. The difference in kinetic behaviour associated with this very subtle structural modification is correctly reproduced by the simulation but, in this case too, τ comp underestimates τ rel (Figure S7). Interestingly, this pair simulated using TMD 37 had equal unbinding mean work within the error bars.</p><!><p>The structures of 1 and 2, the two molecules that belong to Cluster 8, as well as their τ rel and τ comp are given in Figure S8. Their central heterocyclic core is unique in this dataset, but it is decorated similarly to compounds from Clusters 1 and 2. Their only structural difference is in the R1 position, which changes from a methoxy group in compound 1 to a methylmorpholine in compound 2. These structural differences are associated with a large difference in τ rel , which is not reflected by the simulations. Experimentally, the presence of the morpholine substituent lowers k of f by a factor of more than 20 compared to the methoxy: compound 2 is the ligand in our dataset with the lowest k of f , to the point that only an upper limit has been measured. On the other hand, the relative computed residence time of the slowest compound 2 is underestimated as it is found equal to 1 within the estimated error. In the bound state, the crystal structure indicates that the morpholino group of 2 is located outside the binding pocket, with no protein residue closer than 3Å. From Figure 5, the kinetically relevant unbind-ing event happens similarly early for the two compounds. Hence, snapshots of all unbinding trajectories between d = 3 and 3.5Å were extracted for further analysis. Clearly, for both ligands, the main event in these short extracts is the breaking of the hydrogen bond between the central NH of the ring and Thr184, while auxilliary polar interactions, e.g. with the two phenolic OH groups, are often maintained at this point of the unbinding process. However, it is notable that the morpholino group of ligand 2 interacts with residues from α-helix 3, much more than the smaller methoxy of ligand 1. As the conformational space spanned by subset C in the snapshots is large for both compounds, there is no indication of a specific interaction between alpha-helix 3 and 2. In order to better quantify the protein-ligand interactions taking place around the TS, we extracted from circa 20 unbinding trajectories of compounds 1 and 2 the part where d is between 3 and 3.5Å.</p><p>For each of these trajectories we extracted 20 regularly spaced snapshots of the system and computed the mean of the cF P of those snapshots. This procedure was performed on ScaledMD and τ RAMD trajectories from ref. 16 and the results are depicted in Figure 8. The mean cF P of compounds 1 and 2 simulated by ScaledMD, on one hand, and by RAMD, on the other hand, are very similar. However the cF P obtained for ScaledMD simulations and the one obtained for τ RAMD simulations for the same compound are markedly different: whereas the ligand simulated with τ RAMD interacts specifically with its surroundings, the ScaledMD fingerprint is blurred, reflecting the variety of the interactions occuring during the simulations. Moreover, a distinct signature differenciated the unbinding of 2 with respect to 1 using τ RAMD, which was not the case with ScaledMD: 2 in its neutral or protonated form interacts specifically with Ser113 and Gln133 during unbinding, and also with Asp102 when in protonated form. This difference in unbinding patterns translates into the difference in computed residence times obtained with τ RAMD, which correctly predicts a longer τ for 2 than for 1. These observations lead to a plausible interpretation of the underestima- ) in these simulations. The contributions to the activation entropy can be separated between those arising strictly from the protein and from the ligand degrees of freedom, and from contributions which are specific to the protein-ligand interactions, such as the restrictions of the conformational space due to a set of specific H-bonds or steric clashes. The degrees of freedom (DoF) specific to 1 and 2 are very similar, even if the morpholino group, as a non-aromatic ring, can have a DoF of its own. What differs more are the number and the relative population of the possible states in the complete space of the protein-ligand degrees of freedom. Indeed, while the interactions of 1 and 2 with the core of the binding site are very similar in the snapshots, the size of the conformational space spanned by residues of subset C in interactions with the morpholino group of 2 is large for λ = 0.40. It is expected to be much smaller when λ is increased, at least to the point where α-helix 3 keeps its secondary structure. The α-helix 3 behaviour would then be similar to what is observed in conventional or τ RAMD simulations. More generally, we could speculate that the non-linearities in the entropy contribution are responsible for the residence time underestimation in cases where, around the unbinding TS, the ligand interacts non-specifically with parts of the protein that are rendered unstructured by the scaling potential. This hypothesis is testable, but it is beyond the scope of the present study.</p><p>Overall, several trends summarize our findings in terms of residence time estimation. First, the compounds with short residence times were not ranked correctly even when in congeneric series, as seen in Clusters 1, 4 and 6. Second, in the largest chemical series displaying a wide range of τ ( resorcinols from Cluster 2, indazoles from Cluster 3), the trends in predicted residence times are globally correct despite the presence of outliers. Third, in two cases of pairs of compounds with point modifications, the relative residence times are well reproduced while this is not the case for a third pair (Cluster 8) and compound 2, whose residence time, the longest of the whole dataset, is underestimated in the ScaledMD simulations. The very long measured τ of ligand 2 and the underestimation of its residence time might be due to a strong interaction between the morpholino group and the α-helix 3 at the TS that is correctly described with τ RAMD, but not with ScaledMD, where these interactions are blurred, at least with the simulation settings used here.</p><!><p>In this work, we have evaluated the ability of ScaledMD to qualitatively and quantitatively reproduce the residence times of a diverse and sizeable set of ligands binding to two different conformations of their protein target in a context reasonably close to a realistic drug discov-ery project. In order to describe the structural aspects of the unbinding events, we have developed simple procedures for (1) estimating the limits of the basin around the bound state using statistics on the first passage time t f pt (d), and (2) describing and comparing the paths which lead to unbinding using contact-map fingerprints. These procedures are not specific to ScaledMD simulations, they can be used to analyze unbinding trajectories generated with other protocols. For many ligands, the evolution of the squared correlation coefficients R 2 between the unbinding time t of f and the first passage time t f pt (d) can be used to locate the kinetically relevant TSs at 3.5 to 4.5Å from the ligand's binding shell residues. The clustering of the late unbinding parts of the trajectories from our simulations revealed that the two exit routes from the N-HSP90 binding site that were reported previously by several other groups 16,35,37 were indeed followed by most ligands. What is more, each ligand followed both routes in sizeable proportions, and both led to comparable mean exit times.</p><p>In order to estimate the quantitative character of the τ prediction by ScaledMD, we have related τ comp to the measured normalized residence times relative to ligand 3 τ rel,lig3 (λ) = (t of f,ligand /t of f,lig3 ) λ . We found that τ comp allows the identification of the slowest unbinders of the set, and that it gives generally good trends within chemical series even when the ligands bind to different protein conformations. For a subset of the ligands which includes many of the slowest unbinders, the residence time is consistently underestimated in the simulations. For some of them, this underestimation might be due to specific structural or force-field issues, as 9 and the quinazolines 23 and 24 are also underestimated by τ RAMD. On the other hand, the apparent consistency of the underestimation and the fact that it was observed in other works using ScaledMD 35 required a proper investigation.</p><p>We have compared the structure of two homologous compounds, 1 whose τ is correctly predicted, and 2 whose τ is underestimated, around the unbinding TSs, and found a plausible explanation of this underestimation. The bias induced by the scaling of the interactions, together with the choice of residues to keep restrained during the simulations, increase the mobility of α-helix 3 up to a point where its secondary structure elements are no longer conserved, leading to an incorrect description of the interaction taking place and an artefactually low τ comp for ligands that interact with unstructured residues (here α-helix 3) around the unbinding TS.</p><p>Our work is complementary to other recent studies on the unbinding kinetics of N-HSP90 ligands. 16,35,37 It confirms the main results found by Schuetz et al. 35 using ScaledMD on a smaller dataset, such as the identification of two main exit routes or the underestimation of the τ for the slowest binders. With respect to TMD, 37 ScaledMD residence time estimations seem to be less sensitive to the initial protein conformation. On the other hand, the scaling of all interactions in ScaledMD leads to a greater perturbation of the whole system compared to TMD or in τ RAMD 16 where external forces are applied on the ligand. The perturbations due to the scaling makes relevant a detailed analysis of the unbinding trajectories generated by these techniques, for example when comparing possible exit routes. The quantitative performance we obtained in predicting τ using ScaledMD is not on par with τ RAMD's, but it is still reasonable. Both techniques fundamentally rely on one parameter, the magnitude of the random force for τ RAMD and the scaling factor λ for ScaledMD. Nevertheless, the restraints that need to be imposed on parts of the system when using ScaledMD correspond to another set of parameters, whose careful selection is of great importance when setting up ScaledMD simulations. It is actually recommended by Biki to select the restrained residues in a series-specific manner, however, as the comparison of ligands across series was one of the objectives of this work, a common, broader set of unrestrained residues was needed. Second, a higher scaling factor necessarily improves the description of the physics in ScaledMD, at the cost of much longer simulations. Reanalysis of the training dataset used for λ determination actually suggests that this selection might have been biased by the slowest unbinders whose τ comp were underestimated for the reasons explained above. Third, the selective application of the scaling to only a subset of the system is a recent and appealing variation of all-atoms ScaledMD. 57 Selectively scaledMD is a simple restraint-free approach which might alleviate some of the difficulties raised above.</p><p>In a drug design perspective, the use of unbinding simulations might be of use at the lead optimization stage, when high-quality structural information is available to perform highquality low-throughput comparisons within chemical series. It is also tempting to insert unbinding simulations at some point of a virtual screening with the objective of weeding out false positives. Our results suggest that ScaledMD might be useful in this respect as it allowed us to classify correctly most of the compounds with a long τ in our dataset using τ comp . Among the approaches discussed above for the evaluation of protein-ligand binding kinetics, the ones that perform best and are amenable to medium to high throughput calculations will be a welcome addition to the Computer-Aided Drug Design practitioner's toolbox.</p>