Estimation of Protein-Ligand Unbinding Kinetics Using Non-Equilibrium Targeted Molecular Dynamics Simulations
Steffen Wolf, Marta Amaral, Maryse Lowinski, Francois Vallée, Djordje Musil, Jörn Güldenhaupt, Matthias K. Dreyer, Jörg Bomke, Matthias Frech, Jürgen Schlitter, Klaus Gerwert
1 Estimation of Protein-Ligand Unbinding Kinetics Using Non-Equilibrium Targeted Molecular Dynamics Simulations
Steffen Wolf* ,†,‡, , Marta Amaral §,||, È , , Maryse Lowinski ^ , Francois Vallée ^ , Djordje Musil || , Jörn Güldenhaupt † , Matthias K. Dreyer Ð , Jörg Bomke Ç , Matthias Frech || , Jürgen Schlitter † , Klaus Gerwert † † Department of Biophysics, Ruhr-University Bochum, 44780 Bochum, Germany ‡ Institute of Physics, Albert-Ludwigs-University Freiburg, 79104 Freiburg, Germany § Instituto de Biologia Experimental e Tecnológica, 2780-157 Oeiras, Portugal || Molecular Interactions and Biophysics, Merck KGaA, 64293 Darmstadt, Germany È Sanofi-Aventis Deutschland GmbH, Biologics Research / Protein Therapeutics, 65926 Frankfurt am Main, Germany ^ Sanofi IDD-BioStructure and Biophysics, 94400 Vitry-sur-Seine, France Ð Sanofi-Aventis Deutschland GmbH, R&D Integrated Drug Discovery, 65926 Frankfurt am Main, Germany Ç Molecular Pharmacology, Merck KGaA, 64293 Darmstadt, Germany Abstract
We here report on non-equilibrium targeted Molecular Dynamics simulations as tool for the estimation of protein-ligand unbinding kinetics. Correlating simulations with experimental data from SPR kinetics measurements and X-ray crystallography on two small molecule compound libraries bound to the N-terminal domain of the chaperone Hsp90, we show that the mean non-equilibrium work computed in an ensemble of trajectories of enforced ligand unbinding is a promising predictor for ligand unbinding rates. We furthermore investigate the molecular basis determining unbinding rates within the compound libraries. We propose ligand conformational changes and protein-ligand nonbonded interactions to impact on unbinding rates. Ligands may remain longer at the protein if they exhibit strong electrostatic and/or van der Waals interactions with the target. In the case of ligands with rigid chemical scaffold that exhibit longer residence times however, transient electrostatic interactions with the protein appear to facilitate unbinding. Our results imply that understanding the unbinding pathway and the protein-ligand interactions along this path is crucial for the prediction of small molecule ligands with defined unbinding kinetics.
Introduction
While rational drug design traditionally focuses on the optimization of binding affinity of compounds to target proteins, optimization of target binding kinetics is emerging as a new paradigm in drug discovery.
Often, drugs with optimized binding kinetics exhibit better efficacy profiles and reduced off-target toxicity, and thus are more likely to pass later clinical phases. However, while the prerequisites for the rational design of high affinity drugs are well investigated, the rational optimization of kinetic parameters of small molecules is in its early stages. Molecular determinants believed to be important in the modulation of binding kinetics include ligand molecular size, hydrophobic effects, electrostatic interactions, and conformational fluctuations.
Recent reports further highlight the importance of protein-bound water molecules and of protein internal electrostatic interactions. However, the exact contribution and extent of each of these properties still needs to be further elucidated. In order to gain a systematic understanding of the impact of different molecular discriminants on binding kinetics, and thus help to establish a knowledge basis necessary for rational design of compounds with desired kinetics, we performed a combined experimental and theoretical analysis on the dynamics of unbinding of two series of compounds with different chemical scaffolds (see Figure 1) bound to the ATP-binding N-terminal domain of the chaperone heat shock protein 90 (Hsp90, Figure 1C), which is a well-known target for anti-cancer drugs. and preexisting data sets, we included a total of 26 compounds in the present analysis, which are listed in Table 1. Additionally, we determined by X-Ray crystallography the structures of one further protein-ligand complex (see Tables 1 and S2), and measured ligand binding kinetics and affinities of three further compounds via surface plasmon resonance (SPR). In detail, we investigated fourteen compounds with resorcinol backbone (compounds see Figure 1; amongst them the Hsp90 inhibitor Ganetespib ), eleven compounds with N-heterocycle functionalities (compounds ), and the macrocyclic lactam Hsp90 inhibitor 17-DMAG . Figure 1C displays an overview of the N-terminal domain of Hsp90 with bound compound . The binding site is located close to the protein surface, and exhibits two different conformations of the adjacent amino acids 102-114. These residues either form a helix conformation (helix 3) or a loop conformation, which was proposed to affect unbinding kinetics. Table 1 . List of compounds, dynamic properties and protein conformations of investigated compounds. Error bars denote the standard error of the mean (SEM) for N=2-4 measurements. n.d.: not determined. compound Ref. for kinetics and affinity k off / s -1 K D / M -1 k on / M -1 s -1 helix 3 conf. PDB ID with ref. 1a
2 <1.00E-04 <1.00E-09 (n.d.) loop (17-DMAG )
1 3.00E-04 4.57E-09 (n.d.) loop (out)
1 3.30E-04 ±2.1E-05 4.60E+05 ±4.0E-11 2.15E+05 ±5.40E+04 helix (Ganetespib) 1 5.70E-04 1.00E-09 (n.d.) loop
1 1.70E-03 ±4.6E-04 2.30E-08 ±4.4E-09 7.00E+04 ±7.50E+03 helix
1 1.79E-03 ±4.7E-06 3.81E-09 ±3.5E-10 4.72E+05 ±4.10E+04 helix
1 4.20E-03 ±5.2E-04 2.40E-08 ±1.0E-09 1.80E+05 ±2.50E+04 helix
Modelled from M5
1 6.40E-03 ±4.3E-04 8.70E-08 ±2.1E-09 7.70E+04 ±1.20E+04 helix
1 1.40E-02 ±2.2E-03 2.66E-08 ±2.6E-09 5.20E+05 ±1.30E+05 loop
1 1.40E-02 ±1.5E-03 4.30E-07 ±6.1E-08 3.30E+04 ±1.30E+03 helix
1 3.38E-02 ±1.13E-03 7.11E-08 ±4.327E-9 4.79E+05 ±1.65E+04 loop
2 6.34E-02 ±3.5E-03 5.14E-07 ±6.0E-09 1.23E+05 ±5.23E+03 loop
2 1.74E-01 ±2.2E-02 2.36E-07 ±1.9E-08 7.42E+05 ±1.53E+05 helix
1 2.10E-01 ±3.3E-02 1.80E-07 ±1.2E-08 1.20E+06 ±2.10E+05 loop
2 2.54E-01 ±1.8E-02 9.00E-07 ±1.7E-08 2.80E+05 ±1.47E+04 loop (here) 7.10E-05 ±4.0E-06 7.74E-09 ±4.0E-10 (n.d.) helix
2 1.36E-04 ±3.8E-06 8.48E-09 ±6.9E10 1.62E+04 ±1.78E+03 helix
2 1.89E-04 ±7.1E-05 4.66E-08 ±2.5E-08 4.77E+03 ±1.35E+03 helix (here) 1.94E-04 ±9.0E-06 1.01E-08 ±5.0E-10 (n.d.) helix
2 2.78E-04 ±4.65E-06 1.72E-07 ±1.2E-07 3.06E+03 ±2.09E+03 helix
2 2.85E-04 ±4.9E-05 3.61E-08 ±5.7E-09 7.77E+03 ±1.48E+02 helix
2 4.85E-04 ±1.39E-04 3.95E-09 ±1.7E-09 1.33E+05 ±2.37E+04 helix
2 7.65E-04 ±5.0E-05 2.40E-10 ±8.9E-11 3.58E+06 ±1.11E+06 helix
2 9.89E-04 ±1.3E-04 9.50E-08 ±4.5E-09 1.04E+04 ±9.04E+02 helix (here) 3.697E-03 ±1.5E-05 3.285E-8 ±5.61E-9 1.17E+05 ±1.49E4 helix
2 2.56E-02 ±1.4E-02 2.47E-08 ±6.5E-09 1.26E+06 ±8.78E+05 helix Figure 1.
Structure of the N-terminal domain of Hsp90 and of investigated ligand scaffolds. Compound names with two letters denominate alternative protonation states. A: resorcinol compounds ( ) and additional compound (17-DMAG). B: N-heterocycle series ( ) with fluorenamide and benzamide scaffolds and additional compounds and . C: overview of the N-terminal domain of HSP90 in complex with compound . Protein in cartoon, compound in sticks, helix 3 in red, alternate loop conformation in yellow. To assess the molecular mechanisms of unbinding in Hsp90, we performed non-equilibrium targeted molecular dynamics (TMD) simulations.
This method uses a constant velocity constraint as an additional force f c in the simulations to push the ligand out of the binding site. f c is calculated via a Lagrange multiplier with regards to the ligand’s center of mass (COM) and is updated each time step to move the ligand to a position that is in agreement with the preset constant velocity. The constraint force is applied in such a way that the distance between anchor group COM and ligand COM is increased with the preset constant velocity (see Fig. S1). The distance vector acts like a radial vector in a spherical coordinate system, while the ligand is free to move and change direction on the surface of the sphere. This leaves the ligand the freedom to perform diffusion perpendicular to the distance vector, conformational changes and rotations. The ligand thus has the possibility to probe different unbinding pathways, although this choice is limited by the ratio between constraint velocity and diffusion on the sphere. We focus our analysis on the contributions to unbinding kinetics, as unbinding events are easier to calculate than binding events. As we almost exclusively use protein/ligand crystal structures with positions of protein-internal water molecules being resolved, we have an excellent structural basis for carrying out such simulations. As the simulations are carried out under non-equilibrium conditions, i.e., non-stationary with a finite time and velocity, we do not obtain the free energy along the pulling coordinate, but a non-equilibrium work
ΔG ≤ W due to W containing irreversible work caused by dissipation, i.e., friction effects. We find that this non-equilibrium variable is a better predictor for unbinding kinetics than free energy profiles derived from stationary free energy calculations. Methods Chemistry and analytical data of 1j and
Information on the synthesis of chemical compounds is provided in refs. 7,20, on in ref. 25 and their analytical data in Table S1. LC/MS spectra of the products were recorded on an Agilent 1100 HPLC system (1100 high pressure gradient pump, 1100 diode array detector) interfaced to an Agilent 1100 mass spectrometer detector using a Chromolith SpeedROD RP 18e50-4.6 column. Polar gradient: Water (0.05% HCOOH) - acetonitrile (0.04% HCOOH) were used as eluent in mixtures as follows: 0 min, 4% ACN; 2.8 min, 100% ACN; 3.3 min, 100% ACN; Gradient: 5.5 min; Flow-rate: 2.4 ml/min; UV detection: 220 nm. H NMR spectra were recorded at 300 K unless otherwise specified using a Bruker Avance DPX 300, AV 400, DPX 500 spectrometer (TMS as an internal standard). 1H NMR chemical shifts are reported in parts per million (ppm). 1 H NMR data is reported as chemical shift (dH), relative integral, multiplicity (s = singlet, d = doublet, t = triplet, q = quartet, dd = doublet of doublets, ddd = doublet of doublet of doublets, dt = doublet of triplets, td = triplet of doublets, tt = triplet of triplets, qd = quartet of doublets) and coupling constant (J Hz). Both compounds have a purity ≥95%. Information on the synthesis of chemical compounds - and - can be found in patent WO2006087077 and in refs. 7,20,25; for compounds - , and in published patent applications WO2010106290, WO2006123061, WO2008049994 and ref. 19, for compounds and in patent WO2006091963, and in patent WO2005028434. Crystallization and Structure Determination for compound 2d:
A hexa-histidine tagged N-terminal fragment of Hsp90 (18-223) (NP_005339) was expressed and purified as described in ref.
20. Crystallization conditions are also described in ref. 20. Datasets were collected in-house on a Rigaku HF-007 rotating anode generator and a MAR CCD detector and in the synchrotron. Diffraction data were processed with either XDS or MOSFLM. The structures were solved by the molecular replacement method using one set of coordinates of N-HSP90 available in the Protein Data Bank (PDB code: 1YER). The structures were refined using either CNX7, REFMAC5 or BUSTER8 program packages, ligands were placed manually, and the structural models were manually rebuilt using either TURBO-FRODO ( ) or COOT9 . Final validation checks were performed using MOLPROBITY. Surface Plasmon Resonance (SPR) of compounds 2a, 2d and 2j:
SPR measurements were performed on a Biacore 4000 instrument from GE Healthcare as previously described in refs. 7,20. Briefly, recombinant N-HSP90 with 17-desmethoxy-17-N,N-dimethylaminoethylamino-geldanamycin (17-DMAG, Merck Millipore) was immobilized on a Biacore CM5 chip at 25°C at a flow rate of 10 µL/min using amine coupling at pH 4.50 according to Biacore’s standard protocol. HBS-N (10 mM Hepes pH 7.40, 0.15 M NaCl) served as the running buffer during immobilization and all SPR binding kinetics measurements assays were performed in 20 mM HEPES pH 7.50, 150 mM NaCl, 0.05% Tween 20, 1 mM DTT, 0.1 mM EDTA, 2% DMSO. Data sets were processed and analyzed using the Biacore 4000 Evaluation software, version 1.1. Solvent corrected and double-referenced association and dissociation phase data were fitted to a simple 1:1 interaction model with mass transport limitations.
Simulation setup:
TMD calculations were performed with Gromacs v4.6.5 (ref. 38) using the AMBER99SB forcefield for protein and ions, and the TIP3P water model. Crystal structures for compounds and were taken from PDB IDs 1OSF and 3TUH, respectively. Structures of compounds and were taken from PDB IDs 2YKC, 2YKI and 2YKJ. Due to their high similarity, the structure of compound was modeled based on the protein-ligand complex by removing a single terminal methyl group of the respective butenyl side chain. The initial structure of compound was taken from the structure published herein. Crystal structures of all other compounds were determined within the Kinetics for Drug Discovery consortium and are published in refs. 7,20,21 (see Table 1). Ligand parameters were created with antechamber and acpype using GAFF parameters and AM1-BCC charges. Protein/ligand crystal structures together with present crystal water molecules were centered in a cubic box with 7 nm side length, missing protons added, protonated, solvated, and sodium ions added to ensure a charge neutral simulation box. Protonation states of amino acids were determined by propka. Protonation states of other compounds were determined based on pK a values predicted using Chemicalize TMD calculations:
Simulations were carried out with PME for electrostatics (minimal real space cut-off of 1 nm) and a van der Waals cut-off of 1 nm. Hydrogen atom bonds were constrained via the LINCS algorithm. The prepared systems were first minimized with the conjugate gradient method, and subjected to a short equilibration runs in the NPT ensemble at 300 K and 1 bar, using the Berendsen thermostat and barostat, with an integration step size of 2 fs and a trajectory length of 100 ps. For each ligand, 30 statistically independent equilibration runs were performed, in which differed velocity distributions were attributed to the minimized systems to generate an initial equilibrium state ensemble. Non-equilibrium TMD calculations using the Gromacs PULL code in constraint mode were then carried out by continuing the 30 independent equilibration runs for 200 ps in the NPT ensemble at 300 K and 1 bar, using the Nosé-Hoover thermostat and Parrinello-Rahman barostat, with a fixed constraint velocity of 0.01 nm/ps and an integration step size of 1 fs. Constraint pseudoforces were written out for each time step. The 1 st reference group for COM pulling along path 1 consisted of all C(alpha) atoms of the beta-sheet forming the ligand binding site (see Fig. S1) and of all C(alpha) atoms of helix 1 for path 2, the 2 nd group was formed by the ligand heavy atoms. Integrating f c along the pathway as 𝑊(𝑥) = ∫ 𝑓 , (𝑥 - ) 𝑑𝑥′ (1) yields the non-equilibrium work W performed to remove the ligand. In our simulations, we obtain a resolution of 10 fm along x, and calculate W via trapezoidal numeric integration. Trajectory evaluation was then carried out with Gromacs tools, and data evaluation in Python using numpy and scipy libraries and Jupyter notebooks. Stationary thermodynamic integration simulations were performed by extracting 21 equidistant snapshots from a random non-equilibrium simulations, and carrying out equilibration simulations of 10 ns trajectory length with them, setting the constraint velocity to zero (for a detailed explanation see ref. 28). Mean constraint pseudoforces 〈𝑓 , 〉 were calculated from the last 2.5 ns of these simulations. Free energy profiles as given in Fig. S2A were then calculated by integrating the mean forces along the distance from the binding site x as Δ𝐺(𝑥) = ∫ 〈𝑓 , (𝑥′)〉𝑑𝑥′ ≈ ∑ 〈𝑓 ,,9 〉 Δ𝑥 : ; (2) with N x being the number of steps between x and x, and D x the distance between snapshots. Error calculation:
For experimental results, we calculated errors as standard error of the mean (SEM)
Δx = σ/√𝑁 with standard deviation s and number of experiments N . To calculate the errors of theoretical results, we used random sampling bootstrapping with replacement as implemented in the Python scikit-learn library, using 10000 iterations. To keep comparability with the experimental SEM, the displayed error bars represent the 1 s confidence level. Table 2 . Computational results for mean work
1j 1ja
2a 2aa
2b 2ba
2c 2ca
2d 2da
2e 2ea
2f 2fa
2g 2ga
2h 2ha
2i 2ia
2j 2ja
2k 2ka Results & Discussion A linear non-equilibrium energy relationship for unbinding kinetics.
At the beginning of our investigations, we attempted to characterize ligand unbinding kinetics of Hsp90 ligands by determining their free energy profile along the unbinding pathway via standard stationary Thermodynamic Integration (TI) calculations. The most probable unbinding pathway for the ligand appeared to be the passage through an opening between helices 1 and 3 (pathway 1 in Figure S1). Figure S2A displays the resulting free energy surface for compounds , , and , which is in good general agreement with free energy curves for other Hsp90 binding ligands obtained by umbrella sampling. The three investigated compounds exhibit 1-2 free energy barriers between the ligand bound and unbound state. Interpreting the shape and peak height by means of the Eyring equation for rate constants, 𝑘 BCC = 𝜅 exp(−𝛽Δ𝐺 I ) , (3) with the friction-dependent prefactor k , the inverse temperature b = and the free energy difference between bound state and unbinding transition state D G ≠ , we find that effectively does not exhibit a barrier, but a slope between bound and unbound state, and consequently should be the fastest unbinding of the three test compounds. and exhibit a comparable transition barrier of ca. 65 kJ/mol, and both a 2 nd small barrier at 1.8 nm. and should therefore unbind equally fast, which does not agree with the experimental results (see Table 1). Apparently, D G ≠ is not sufficient as descriptor for predicting unbinding kinetics, and would require taking into account the prefactor k in Equation (3), which according to Kramers includes friction effects. A further problem we faced when applying stationary TI calculations was the large number of necessary equilibration points along the unbinding pathway that need several nanoseconds of equilibration for reliable determination of the free energy surface, significantly raising the computational cost for investigating a large set of compounds. Furthermore, in our two investigated compound groups, about half of all compounds exhibit two possible protonation states ( , , and the full series ). As an example, the morpholine side chain in can exist in a protonated state with a charge of +1 e (see Figure 1A), or in a deprotonated state with a charge of 0 e . All ligands in compound group are bound to the protein by a hydrogen bond between nitrogen atoms in aromatic rings (pK a range of ca. 1.5-5) and Asp93 (see Figure S3), or via highly polarized water molecules mediating this contact. Assigning the correct protonation state for such protein-ligand-water complexes is a challenging task, as the protein environment can significantly alter pK a values. To avoid a bias from wrongly chosen charge states, we needed a method that allowed us to carry out simulations of multiple compounds in 2-3 possible protonation states, with TI calculations simply being too inefficient for this task. Surprisingly, when we looked at the mean non-equilibrium work profiles
K9LL (6) with dissipative work W diss , which contains the second and all higher moments. After an increase of
𝐶 = 𝛽 MJ ln 𝜅 + 𝑊 K9LL , which serves as a basis of understanding the apparent linear non-equilibrium energy relationship. C effectively is a function of b , but in the following is treated as an independent fit factor, as we otherwise encountered instabilities in non-linear curve fitting. In the following, we approximate C to be constant, which is only valid in the case that the friction during unbinding is the same for all ligands. Figure 2.
Comparison of experimentally derived k off constants and calculated TMD work
Influence of protein conformation and electrostatics.
As a starting point for investigating molecular effects influencing unbinding rates, we focused on a dependence on helix/loop 3 conformation as implied by our analysis in Fig 2. In our simulations, helix binding compounds with decreasing k off display an increasing unbinding
Influence of ligand radius on unbinding kinetics of the resorcinol series compounds. Group helix binders in red, loop binders in blue, group in cyan. Linear regression as full lines. Line colors match the color of data points used for fits. A: difference of radii of gyration in respect to the radius in the unbound state integrated over pulling distance. B: mean radius of gyration during unbinding. C: absolute integral of radii of gyration difference to unbound state over the full course of simulation in the N-heterocycle compound series . Protonation states were chosen according to pK a predictions. Fit excluding outliers as full line, fit with outliers as dashed line. Based on the predicted pK a values in Table 2, , and should be the relevant structures at pH = 7.5. However, considering pH values at which the alternative protonation states claim 10% of the full ligand population ( : 7.1; : 8.8; : 5.2), we realize that a small population of alternatively protonated forms could be present within an ensemble of protein-ligand complexes. Following our line of argumentation on path selection, these deprotonated ligand forms with smaller
Electrostatic facilitation in compound
A: Ligand dihedral potential energies (referenced to the mean value of the last 0.5 nm and smoothed with a Gauss filter with s = 0.2 nm) and radius of gyration as average of N = 30 simulations. Trajectory mean as lines, 1 s error level from bootstrap analysis (see Methods) as shaded area. Uncharged ligand in red, protonated form in yellow. B: molecular details of ligand unbinding. A favorable charge interaction of the azaindoline moiety with Asp102 facilitates the contraction of the ligand from (1) to (2), and guides the ligand out of the binding pocket into the unbound conformation (3). Electrostatic locking vs. facilitation.
Focusing on the effect of electrostatics on
In recent years, several other novel methods have been established for fast and efficient computation of binding kinetics (see refs. 73,74 for reviews), and our approach presented here shares similarities with methods based on metadynamics and steered MD. While the calculation of a constraint numerically is more complex than the addition of an additional bias potential, the major advantage of our method is that employing a constraint allows to scan the preset pulling coordinate with a linear velocity that is accurate down to machine precision. Furthermore, we experienced in other works that employing biasing potentials causes problems at steep potential gradients, causing long simulation times or even simulation crashes. Such problems are overcome by employing constraints, as the biasing force f c exactly cancels out such gradients. As a prerequisite, we need to have initial information on unbinding pathways to create a suitable reaction coordinate to apply the target bias: as stated in the Introduction, while the ligand is free to diffuse perpendicular to the pulling vector, this dynamics is restricted by the ratio between diffusion rate and constraint velocity. At the start of the simulation, the initial pulling vector thus needs to roughly point into the direction of the presumed binding pathway. Besides taking educated guesses, this information can be provided by other methods that employ more general reaction coordinates. It was recently shown that TMD simulations can be used to effectively push a molecular system of choice along a reaction coordinate that correctly mimics the pathway taken under equilibrium conditions. The first major strength of our non-equilibrium method is the significant reduction of necessary computational power: unbinding can be enforced within 0.1 to 0.3 ns simulated time, which allows us to reduce the necessary calculation time by a factor of 30 in comparison to stationary TI (180 CPUhs for one non-equilibrium unbinding event sampled with 30 trajectories vs. 6000 CPUhs for one equilibrium free energy pathway analyzed by 20 intermediate steps on a recent octacore CPU workstation). Secondly, the non-equilibrum work rapidly converges (see Figure S4), and each simulation by definition results in an unbinding event, which reduces the necessary number of simulations to a number well below that for Markov State Model creation. Thirdly, we do not change the full system Hamiltonian, but merely add a perturbation, avoiding artifacts such as protein unfolding that appear in smoothed/scaled MD.
Summarizing the results of our investigation, it becomes clear that estimation of binding / unbinding rates and elucidating the relevant underlying molecular mechanisms is far more complex than free energy calculations of binding poses, as it not only require to correctly characterize a single binding pose, but to assess the full dynamics along binding/unbinding pathways. While binding simulations apparently are more easy to carry out due to the shorter inherent timescales, special care needs to be taken that the simulations result in the correct binding pose. In this respect, biased unbinding simulations are simpler, as they only require enforced unbinding, and it is always “simpler to break things than to make things”. However, to gain meaningful results that can be compared to experimental data, the biased simulations need to reproduce the correct protein-ligand dynamics, i.e., ligands need to leave the protein along the correct pathway, which needs to be identified first. Another challenge is found in the theoretical basis of the Jarzynski equality (4) itself: biased simulations need to start from a representative equilibrium ensemble of the employed protein. In our case, we were interested in developing a “fast” method for prediction of unbinding kinetics, and thus performed only short equilibration runs of 0.1 ns before initiating unbinding. This requires that the employed protein crystal structure is close to an equilibrium structure. Furthermore, if unbinding kinetics are coupled to protein conformational dynamics, e.g., by only shortly opening a binding site to allow unbinding, our approach will not succeed. Besides the results presented here, we found that flexible proteins with challenging unbinding pathways pose a problem for our non-equilibrium TMD method, as a similar investigation with ligands bound to the b adrenergic receptor performed by us did not succeed in obtaining a successful linear non-equilibrium energy relationship (data not shown). The reason for this appears to be the presence of a second intermediate ligand binding site that increases the complexity of the underlying ligand diffusion pathways, and protein conformational dynamics. Last, as the details of ligand conformational dynamics and protein-ligand charge interactions during unbinding seem to be crucial for correct predictions, small errors in ligand parameterization can cause significant problems, which require elaborate parameterization schemes. As we want to be able to perform calculations with a larger number of ligands, we here use a quick standard parameterization scheme via antechamber using semi-empirical charges. It thus may be that a part of the spread of data points around linear fits in Fig. 2B is a result from errors in dihedral angle potentials and atomic charge assignment. In RAMD simulations on Hsp90 it was indeed shown that including charge details like halogen s -holes improves the prediction of unbinding kinetics. As the problems listed above are principal effects coming from the underlying Hamiltonian dynamics of the protein-ligand complexes, other biased simulations approaches will face similar challenges. Conclusion and future perspective
To elucidate the molecular determinants for unbinding kinetics, we here combined preexisting and novel data from SPR binding kinetics measurements and X-ray crystallography with non-equilibrium targeted MD simulations on the N-terminal domain of Hsp90 for two compound series. The non-equilibrium work
Three supporting tables, seven supporting figures (PDF) SMILES annotations (CSV)
Accession Codes
The crystallographic coordinates of novel compound are deposited in the Protein Data Bank under the accession code 5LRL. Authors will release the atomic coordinates and experimental data upon article publication. Author Information Corresponding Author Information *phone: +49-761-203-5913 eMail: [email protected]
Author contributions
S.W. and M.A. contributed equally to this work. S.W. designed and performed all computational studies. M.A. designed and performed SPR and X-ray studies. M.L., F.V. and D.M. solved selected crystal structures. J.G. assisted in data analysis. J.B., M.K.D. M.F., J.S. and KG supervised the studies. All authors contributed to writing of the manuscript.
Acknowledgements This work was supported by the EU/EFPIA Innovative Medicines Initiative (IMI) Joint Undertaking K4DD (grant No. 115366). 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. We acknowledge the Partnership for Advanced Computing in Europe for awarding us access at CINECA Italy (project No. 2015133089). We are grateful to M. Bianciotto, D. Kohk and R. Wade for helpful discussions.
Abbrevations Used
Hsp90: heat shock protein 90; MD: molecular dynamics; SPR: surface plasmon resonance; TMD: targeted molecular dynamics References (1) Swinney, D. C. Opinion: Biochemical Mechanisms of Drug Action: What Does It Take for Success?
Nat. Rev. Drug Disc. , , 801–808. (2) Copeland, R. A.; Pompliano, D. L.; Meek, T. D. Drug-Target Residence Time and Its Implications for Lead Optimization. Nat. Rev. Drug Disc. , (9), 730–739. (3) Lu, H.; Tonge, P. J. Drug-Target Residence Time: Critical Information for Lead Optimization. Curr. Opin. Chem. Biol. , , 467–474. (4) Pan, A. C.; Borhani, D. W.; Dror, R. O.; Shaw, D. E. Molecular Determinants of Drug–Receptor Binding Kinetics. Drug Discov. Today , (13-14), 667–673. (5) Romanowska, J.; Kokh, D. B.; Fuller, J. C.; Wade, R. C. Computational Approaches for Studying Drug Binding Kinetics. In onlinelibrary.wiley.com ; Methods and Principles in Medicinal Chemistry; Wiley-VCH Verlag GmbH & Co. KGaA: Weinheim, Germany, 2015; Vol. 12, pp 211–235. (6) 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.; Müller-Fahrnow, A.; Ecker, G. F. Kinetics for Drug Discovery: an Industry-Driven Effort to Target Drug Residence Time. Drug Discov. Today , , 896–911. (7) Amaral, M.; Kokh, D. B.; Bomke, J.; Wegener, A.; Buchstaller, H. P.; Eggenweiler, H. M.; Matias, P.; Sirrenberg, C.; Wade, R. C.; Frech, M. Protein Conformational Flexibility Modulates Kinetics and Thermodynamics of Drug Binding. Nat. Commun. , , 2276. (8) Swinney, D. C. Applications of Binding Kinetics to Drug Discovery. Pharm. Med. , , 23–34. (9) Sukkar, E. Bound to Work Better: Binding Kinetics Can Be Used to Better Inform the Design and Development of Drugs. Pharmaceut. J. , , 7818. (10) Klebe, G. Applying Thermodynamic Profiling in Lead Finding and Optimization. Nat. Rev. Drug Disc. , , 95–110. (11) Klebe, G. The Use of Thermodynamic and Kinetic Data in Drug Discovery: Decisive Insight or Increasing the Puzzlement? ChemMedChem , , 229–231. (12) Bortolato, A.; Deflorian, F.; Weiss, D. R.; Mason, J. S. Decoding the Role of Water Dynamics in Ligand-Protein Unbinding: CRF1R as a Test Case. J. Chem. Inf. Model. , , 1857–1866. (13) Segala, E.; Guo, D.; Cheng, R. K. Y.; Bortolato, A.; Deflorian, F.; Doré, A. S.; Errey, J. C.; Heitman, L. H.; IJzerman, A. P.; Marshall, F. H.; Cooke, R. M. Controlling the Dissociation of Ligands From the Adenosine A2A Receptor Through Modulation of Salt Bridge Strength. J. Med. Chem. , , 6470–6479. (14) Taipale, M.; Jarosz, D. F.; Lindquist, S. HSP90 at the Hub of Protein Homeostasis: Emerging Mechanistic Insights. Nat. Rev. Mol. Cell Biol. , , 515–528. (15) Ratzke, C.; Hellenkamp, B.; Hugel, T. Four-Colour FRET Reveals Directionality in the Hsp90 Multicomponent Machinery. Nat. Commun. , , 4192. (16) Verba, K. A.; Wang, R. Y.-R.; Arakawa, A.; Liu, Y.; Shirouzu, M.; Yokoyama, S.; Agard, D. A. Atomic Structure of Hsp90-Cdc37-Cdk4 Reveals That Hsp90 Traps and Stabilizes an Unfolded Kinase. Science , , 1542–1547. (17) Whitesell, L.; Lindquist, S. L. HSP90 and the Chaperoning of Cancer. Nat. Rev. Cancer , , 761–772. (18) Trepel, J.; Mollapour, M.; Giaccone, G.; Neckers, L. Targeting the Dynamic HSP90 Complex in Cancer. Nat. Rev. Cancer , , 537–549. (19) Vallée, F.; Carrez, C.; Pilorge, F.; Dupuy, A.; Parent, A.; Bertin, L.; Thompson, F.; Ferrari, P.; Fassy, F.; Lamberton, A.; Thomas, A.; Arrebola, R.; Guerif, S.; Rohaut, A.; Certal, V.; Ruxer, J.-M.; Delorme, C.; Jouanen, A.; Dumas, J.; Grépin, C.; Combeau, C.; Goulaouic, H.; Dereu, N.; Mikol, V.; Mailliet, P.; Minoux, H. Tricyclic Series of Heat Shock Protein 90 (Hsp90) Inhibitors Part I: Discovery of Tricyclic Imidazo[4,5- C]Pyridines as Potent Inhibitors of the Hsp90 Molecular Chaperone. J. Med. Chem. , , 7206–7219. (20) Kokh, D. B.; Amaral, M.; Bomke, J.; Grädler, U.; Musil, D.; Buchstaller, H.-P.; Dreyer, M. K.; Frech, M.; Lowinski, M.; Vallée, F.; Bianciotto, M.; Rak, A.; Wade, R. C. Estimation of Drug-Target Residence Times by Τ-Random Acceleration Molecular Dynamics Simulations. J. Chem. Theory Comput. , , 3859–3869. (21) Güldenhaupt, J.; Amaral, M.; Kötting, C.; Schartner, J.; Musil, D.; Frech, M.; Gerwert, K. Ligand-Induced Conformational Changes in HSP90 Monitored Time Resolved and Label Free-Towards a Conformational Activity Screening for Drug Discovery. Angew. Chem. Int. Ed. Engl. , , 9955–9960. (22) Jez, J. M.; Chen, J. C. H.; Rastelli, G.; Stroud, R. M.; Santi, D. V. Crystal Structure and Molecular Modeling of 17-DMAG in Complex with Human Hsp90. Chem. Biol. , , 361–368. (23) Ying, W.; Du, Z.; Sun, L.; Foley, K. P.; Proia, D. A.; Blackman, R. K.; Zhou, D.; Inoue, T.; Tatsuta, N.; Sang, J.; Ye, S.; Acquaviva, J.; Ogawa, L. S.; Wada, Y.; Barsoum, J.; Koya, K. Ganetespib, a Unique Triazolone-Containing Hsp90 Inhibitor, Exhibits Potent Antitumor Activity and a Superior Safety Profile for Cancer Therapy. Mol. Cancer Therapeut. , , 475–484. (24) Proia, D. A.; Bates, R. C. Ganetespib and HSP90: Translating Preclinical Hypotheses Into Clinical Promise. Cancer Res. , , 1294–1300. (25) Schuetz, D. A.; Richter, L.; Amaral, M.; Grandits, M.; Grä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. , , 4397–4411. (26) Schlitter, J.; Engels, M.; Krüger, P. Targeted Molecular Dynamics: a New Approach for Searching Pathways of Conformational Transitions. J. Mol. Graph. , , 84–89. (27) Schlitter, J.; Swegat, W.; Mülders, T. Distance-Type Reaction Coordinates for Modelling Activated Processes. J. Mol. Model. , , 171–177. (28) Ernst, M.; Wolf, S.; Stock, G. Identification and Validation of Reaction Coordinates Describing Protein Functional Motion: Hierarchical Dynamics of T4 Lysozyme. J. Chem. Theory Comput. , , 5076–5088. (29) Wolf, S.; Stock, G. Targeted Molecular Dynamics Calculations of Free Energy Profiles Using a Nonequilibrium Friction Correction. J. Chem. Theory Comput. , , 6175–6182. (30) Betz, R. M.; Dror, R. O. How Effectively Can Adaptive Sampling Methods Capture Spontaneous Ligand Binding? J. Chem. Theory Comput. , , 2053–2063. (31) Kabsch, W. Xds. Acta Crystallogr. , , 125–132. (32) Battye, T. G. G.; Kontogiannis, L.; Johnson, O.; Powell, H. R.; Leslie, A. G. W. iMOSFLM: a New Graphical Interface for Diffraction-Image Processing with MOSFLM. Acta Crystallogr. D , , 271–281. (33) Brünger, A. T.; Adams, P. D.; Clore, G. M.; DeLano, W. L.; Gros, P.; Grosse-Kunstleve, R. W.; Jiang, J. S.; Kuszewski, J.; Nilges, M.; Pannu, N. S.; Read, R. J.; Rice, L. M.; Simonson, T.; Warren, G. L. Crystallography & NMR System: a New Software Suite for Macromolecular Structure Determination. Acta Crystallogr. D , , 905–921. (34) Murshudov, G. N.; Vagin, A. A.; Dodson, E. J. Refinement of Macromolecular Structures by the Maximum-Likelihood Method. Acta Crystallogr. D , , 240–255. (35) Bricogne, G.; Blanc, E.; Brandl, M.; Flensburg, C.; Keller, P.; Paciorek, W.; Roversi, P.; Sharff, A.; Smart, S. O.; Vonrhein, C.; Womack, O. T. Buster Version 2.11.6. . (36) Emsley, P.; Lohkamp, B.; Scott, W. G.; Cowtan, K. Features and Development of Coot. Acta Crystallogr. D , , 486–501. (37) Chen, V. B.; Arendall, W. B.; Headd, J. J.; Keedy, D. A.; Immormino, R. M.; Kapral, G. J.; Murray, L. W.; Richardson, J. S.; Richardson, D. C. MolProbity: All-Atom Structure Validation for Macromolecular Crystallography. Acta Crystallogr. D , , 12–21. (38) Pronk, S.; Páll, S.; Schulz, R.; Larsson, P.; Bjelkmar, P.; Apostolov, R.; Shirts, M. R.; Smith, J. C.; Kasson, P. M.; van der Spoel, D.; Hess, B.; Lindahl, E. GROMACS 4.5: a High-Throughput and Highly Parallel Open Source Molecular Simulation Toolkit. Bioinformatics , , 845–854. (39) Hornak, V.; Abel, R.; Okur, A.; Strockbine, B.; Roitberg, A.; Simmerling, C. Comparison of Multiple Amber Force Fields and Development of Improved Protein Backbone Parameters. Proteins , , 712–725. (40) Best, R. B.; Hummer, G. Optimized Molecular Dynamics Force Fields Applied to the Helix−Coil Transition of Polypeptides. J. Phys. Chem. B , , 9004–9015. (41) 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. , , 926–935. (42) Wang, J.; Wang, W.; Kollman, P. A.; Case, D. A. Automatic Atom Type and Bond Type Perception in Molecular Mechanical Calculations. J. Mol. Graph. Model. , , 247–260. (43) Sousa da Silva, A. W.; Vranken, W. F. ACPYPE - AnteChamber PYthon Parser interfacE. BMC Res. Notes , , 367. (44) Wang, J. M.; Wolf, R. M.; Caldwell, J. W.; Kollman, P. A.; Case, D. A. Development and Testing of a General Amber Force Field. J. Comput. Chem. , , 1157–1174. (45) Jakalian, A.; Bush, B. L.; Jack, D. B.; Bayly, C. I. Fast, Efficient Generation of High-Quality Atomic Charges. AM1-BCC Model: I. Method. J. Comput. Chem. , , 132–146. (46) Jakalian, A.; Jack, D. B.; Bayly, C. I. Fast, Efficient Generation of High-Quality Atomic Charges. AM1-BCC Model - II. Parameterization and Validation. J. Comput. Chem. , , 1623–1641. (47) Olsson, M. H.; Søndergaard, C. R.; Rostkowski, M.; Jensen, J. H. PROPKA3: Consistent Treatment of Internal and Surface Residues in Empirical P K a Predictions. J. Chem. Theory Comput. , (49) Darden, T.; York, D.; Pedersen, L. Particle Mesh Ewald: an N ⋅ Log( N) Method for Ewald Sums in Large Systems.
J. Chem. Phys. , , 10089–10092. (50) Hess, B.; Bekker, H.; Berendsen, H. J. C.; Fraaije, J. G. E. M. LINCS: a Linear Constraint Solver for Molecular Simulations. J. Comput. Chem. , , 1463–1472. (51) Berendsen, H. J. C.; Postma, J. P. M.; van gunsteren, W. F.; DiNola, A.; Haak, J. R. Molecular Dynamics with Coupling to an External Bath. J. Chem. Phys. , , 3684–3690. (52) Nosé, S. A Molecular Dynamics Method for Simulations in the Canonical Ensemble. Mol. Phys. , , 255–268. (53) Hoover, W. G. Canonical Dynamics: Equilibrium Phase-Space Distributions. Phys. Rev. A , , 1695–1697. (54) Parrinello, M.; Rahman, A. Polymorphic Transitions in Single Crystals: a New Molecular Dynamics Method. J. Appl. Phys. , , 7182–7190. (55) Oliphant, T. E. Python for Scientific Computing. Comput. Sci. Eng. , 10–20. (56) Kluyver, T.; Ragan-Kelley, B.; Pérez, F.; Granger, B.; Bussonier, M.; Frederic, J.; Kelley, K.; Hamrick, J.; Grout, J.; Corlay, S.; Ivanov, P.; Avila, D.; Abdalla, S.; Willing, C.; Team, J. D. Jupyter Notebooks-a Publishing Format for Reproducible Computational Workflows. In Positioning and Power in Academic Publishing ; Loizides, F., Schmidt, B., Eds.; books.google.com, 2016; pp 87–90. (57) Berendsen, H.
Simulating the Physical World ; Cambridge University Press, 2007. (58) Pedregosa, F.; Varoquaux, G.; Gramfort, A.; Michel, V.; Thirion, B.; Grisel, O.; Blondel, M.; Prettenhofer, P.; Weiss, R.; Dubourg, V.; Vanderplas, J.; Passos, A.; Cournapeau, D.; Brucher, M.; Perrot, M.; Duchesnay, E. Scikit-Learn: Machine Learning in Python.
J. Mach. Learn. Res. , , 2825–2830. (59) Hughes, I.; Hase, T. Measurements and Their Uncertainties ; Oxford University, Press, 2010. (60) Eyring, H. The Activated Complex in Chemical Reactions.
J. Chem. Phys. , , 107–115. (61) Kramers, H. A. Brownian Motion in a Field of Force and the Diffusion Model of Chemical Reactions. Physica , , 284–304. (62) Wolf, S.; Freier, E.; Cui, Q.; Gerwert, K. Infrared Spectral Marker Bands Characterizing a Transient Water Wire Inside a Hydrophobic Membrane Protein. J. Chem. Phys. , , 22D524. (63) Wolf, S.; Freier, E.; Gerwert, K. A Delocalized Proton-Binding Site Within a Membrane Protein. Biophys. J. , , 174–184. (64) Gohlke, H.; Klebe, G. Approaches to the Description and Prediction of the Binding Affinity of Small ‐ Molecule Ligands to Macromolecular Receptors.
Angew. Chem. Int. Ed. Engl. , , 2644–2676. (65) Jarzynski, C. Nonequilibrium Equality for Free Energy Differences. Phys. Rev. Lett. , , 2690–2693. (66) Hummer, G.; Szabo, A. Free Energy Surfaces From Single-Molecule Force Spectroscopy. Acc. Chem. Res. , , 504–513. (67) Post, M.; Wolf, S.; Stock, G. Principal Component Analysis of Nonequilibrium Molecular Dynamics Simulations. , , 204110. (68) Hendrix, D. A.; Jarzynski, C. A “Fast Growth” Method of Computing Free Energy Differences. J. Chem. Phys. , , 5974. (69) Park, S.; Khalili-Araghi, F.; Tajkhorshid, E.; Schulten, K. Free Energy Calculation From Steered Molecular Dynamics Simulations Using Jarzynski's Equality. J. Chem. Phys. , , 3559–3566. (70) Schuetz, D. A.; Bernetti, M.; Bertazzo, M.; Musil, D.; Eggenweiler, H.-M.; Recanatini, M.; Masetti, M.; Ecker, G. F.; Cavalli, A. Predicting Residence Time and Drug Unbinding Pathway Through Scaled Molecular Dynamics. J. Chem. Inf. Model. , , 535–549. (71) Zou, K. H.; O’Malley, A. J.; Mauri, L. Receiver-Operating Characteristic Analysis for Evaluating Diagnostic Tests and Predictive Models. Circulation . (72) Schneider, M.; Wolf, S.; Schlitter, J.; Gerwert, K. The Structure of Active Opsin as a Basis for Identification of GPCR Agonists by Dynamic Homology Modelling and Virtual Screening Assays.
FEBS Lett. , , 3587–3592. (73) Chen, Y.-C. Beware of Docking! Trends Pharmacol. Sci. , , 78–95. (74) Potterton, A.; Husseini, F. S.; Southey, M. W. Y.; Bodkin, M. J.; Heifetz, A.; Coveney, P. V.; Townsend-Nicholson, A. Ensemble-Based Steered Molecular Dynamics Predicts Relative Residence Time of a 2AReceptor Binders. J. Chem. Theory Comput. , , 3316–3330. (75) Oostenbrink, C.; van Gunsteren, W. F. Calculating Zeros: Non-Equilibrium Free Energy Calculations. Chem. Phys. , , 102–108. (76) Plattner, N.; Doerr, S.; De Fabritiis, G.; Noé, F. Complete Protein–Protein Association Kinetics in Atomic Detail Revealed by Molecular Dynamics Simulations and Markov Modelling. Nat. Chem. . (77) 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. , , 11539. (78) 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. , , 1–10. (79) Rasmussen, S. G. F.; Choi, H.-J.; Rosenbaum, D. M.; Kobilka, T. S.; Thian, F. S.; Edwards, P. C.; Burghammer, M.; Ratnala, V. R. P.; Sanishvili, R.; Fischetti, R. F.; Schertler, G. F. X.; Weis, W. I.; Kobilka, B. K. Crystal Structure of the Human Β2 Adrenergic G-Protein-Coupled Receptor. Nature , , 383–387. (80) Dror, R. O.; Pan, A. C.; Arlow, D. H.; Borhani, D. W.; Maragakis, P.; Shan, Y.; Xu, H.; Shaw, D. E. Pathway and Mechanism of Drug Binding to G-Protein-Coupled Receptors. Proc. Natl. Acad. Sci. USA , , 13118–13123. (81) Schaudinnus, N.; Bastian, B.; Hegger, R.; Stock, G. Multidimensional Langevin Modeling of Nonoverdamped Dynamics. Phys. Rev. Lett. ,115