|Home | About | Journals | Submit | Contact Us | Français|
We determine the binding mode of a macrocyclic radicicol-like oxime to yeast HSP90 by combining computer simulations and experimental measurements. We sample the macrocyclic scaffold of the unbound ligand by parallel tempering simulations and dock the most populated conformations to yeast HSP90. Docking poses are then evaluated by the use of binding free energy estimations with the linear interaction energy method. Comparison of QM/MM-calculated NMR chemical shifts with experimental shift data for a selective subset of back-bone 15N provides an additional evaluation criteria. As a last test we check the binding modes against available structure-activity-relationships. We find that the most likely binding mode of the oxime to yeast HSP90 is very similar to the known structure of the radicicol-HSP90 complex.
The heat shock protein 90 (HSP90) is an ATP-dependent chaperone whose activity is required for the functioning of a number of proteins that promote the growth (or survival) of cells. [1, 2, 3] In the absence of the chaperoning activity by HSP90, the oncoproteins are targeted for degradation. The prospect of stalling several oncoproteins by inhibiting a single protein has made HSP90 an attractive target for chemotherapy. [4, 5] The natural product radicicol (1, Fig. 1), a 14-membered macrolide, is known to be a potent inhibitor of HSP90’s chaperone activity. Crystal structures and NMR studies of HSP90 bound with radicicol show that radicicol binds to the ATP-binding pocket even though it lacks structural or topological similarity with ATP. [6, 7, 8] While radicicol is a potent inhibitor of HSP90, it is not active in vivo presumably due to metabolic instability.
We aim to design new inhibitors for HSP90 that are based on the scaffold of radicicol but have metabolic stability. The oximes 2 and 3 fulfill this requirement and indeed show in vivo activity.  To optimize the binding specificity of these inhibitors, it is crucial to know where and how they bind to the protein; i.e., what the interactions of an inhibitor with the protein are and what is its conformation. This is termed the binding mode of the inhibitor. Once the binding mode is established, rational design is feasible, i.e., modifications of the inhibitor that would increase the binding affinity can be designed in silico. This is a common and often successful procedure in structure-based drug design. 
The structural characterization of the binding mode is therefore a key step in the process of optimizing inhibitors 2 and 3. From an experimental point of view, crystallography and nuclear magnetic resonance (NMR) spectroscopy are the ideal methods for addressing this question. Experimental measurements of the 15N and 1H chemical shifts of the protein backbone are available for yeast HSP90 in the presence of the ligands 1 and 3 (the configuration of the C=N double bond is E). In principle a complete NMR structure determination for a protein of this size is possible, but it is considerably more time consuming (than measurement of the chemical shifts) and so far has not been performed.
To determine the possible binding mode of a ligand, we present here an approach that combines experimental data with computational methods. The calculations involve enhanced sampling of ligand conformations, docking, molecular dynamics simulations and quantum mechanics calculations; the experiments include chemical shift measurements and structure-activity-relationships (SAR).
Previous computational and NMR studies [11, 12] recognized the flexibility of the macrocyclic scaffold as an important aspect for target binding. The macrocycle can adopt several conformations in solution, including the bioactive conformation. The population of these conformations influences the activity of the ligand. We therefore sample the conformations accessible to the ligand 3 in solution and test their binding interaction with HSP90; i.e., we dock the solution conformations to the protein and obtain a set of feasible binding modes. Several criteria are then applied to evaluate those modes. We estimate the relative free energy of binding from the docking scoring function and from molecular dynamics simulations using the linear interaction energy (LIE) approach.  Furthermore we calculate for a subset of protein residues the 15N NMR shifts with a QM/MM-Hamiltonian  and compare the results with experimental data. The possible binding modes are also evaluated based on their consistency with SAR data.
We apply the approach to inhibitor 3 bound to yeast HSP90. Since both X-ray structure and chemical shift data are available for radicicol (1) bound to HSP90, the known binding mode of 1 serves as test case, in particular for the NMR calculations.
Most of the studies on HSP90, including this one, have been performed on the yeast protein. It is clear, however, that from a therapeutic point of view the human protein is of primary interest. We consider therefore the question of transferability of our results to the human protein.
Because the available experiments (chemical shift measurements and SAR) do not permit a determination of the binding mode of E-oxime 3 to HSP90, we supplement the experiments with computational techniques to obtain a unique answer. The characteristics of the approach are summarized in this section; computational details in the application are given in the following sections.
We dock 3 to HSP90 and evaluate the resulting binding modes by four criteria:
Based on these criteria we determine the most probable binding mode(s) of E-oxime 3, i.e., the mode(s) which satisfy the largest fraction of criteria 1)- 4).
We modeled the Hamiltonian of the ligands by the Merck Molecular Force Field (MMFF94) [20, 21, 22, 23, 24] as implemented in the program CHARMM . A dielectric constant of 80 was used to simulate the effect of the solvent;  since no charged groups are involved such a simple treatment of solvation is expected to be sufficiently accurate.
The conformational space of radicicol and E-oxime 3 was sampled by a MD-based parallel tempering (PT) [26, 27, 28] approach with 12 heat baths (T1, T2, …, T11, T12 [in K]: 285, 300, 335, 370, 425, 480, 555, 630, 725, 820, 920, 1020). The temperature of the baths was controlled via Langevin dynamics with a friction coefficient of 40 ps−1. The integration time step was 1 fs; the mass of the hydrogen atoms was set to 10 amu to avoid integration problems at very high temperatures. Swapping between heat baths was attempted every 10 ps; the total length of the PT simulation was 200 ns (20’000 attempts). We saved snapshots every 1 ps for the last 150 ns (50 ns served as equilibration) for the system at 300 K. This yielded equilibrated ensembles of 150’000 structures.
The structures of the ensembles were split into conformational clusters with the leader algorithm . A root-mean-square (RMS) deviation of 0.5 Å (when comparing the non-hydrogen atoms of the macrocycle) was used as threshold criterion for the clustering. For the E-oxime simulation a small percentage (< 8 %) of structures showed a Z-configuration of the oxime C=N double bond, i.e., isomerization of the oxime occurred due to the high temperature heat baths used in the PT simulation. Those structures were deleted from the clusters.
The clusters were ranked according to their population with the first cluster having the highest population. The relative free energy of the clusters was estimated using , with k the Boltzmann constant, T the temperature (300 K), Pi the population of cluster i, and P1 the population of the first cluster; i.e., ΔGrel 0 for the first cluster. In total we obtained 23 clusters with a free energy range of 0–5.8 kcal mol−1 for radicicol and 31 clusters ranging from 0 to 5.7 kcal mol−1 for E-oxime 3, respectively.
From the low free energy clusters (i.e., clusters that are not higher than 3 kcal mol−1 in free energy than the highest populated cluster; there are six such clusters for radicicol and four clusters for E-oxime 3 respectively), the conformation with the lowest potential energy was picked as a representative and minimized by 2000 steps of steepest-descent followed by 5000 steps of adapted-basis Newton-Raphson minimization.
The docking was performed with the program suite AutoDock 4 [15, 16]. Each conformation of the ligands (see below) was docked to each HSP90 structure (section “Preparation of the protein”) with default AutoDock settings; except that docking was repeated 100 times (instead of 10 times). The reported value in Tab. 1 (denoted with “AutoDock”) is the mean docked energy.
The docking of the ligands 1 and 3 needed to be done with fixed macrocycle conformations (see below, section “Preparation of the ligands”) and therefore the contribution of the relative stability of the macrocycle conformation to the binding affinity was missing in the scoring function. To correct for this, we summed up the mean docked energy and the conformational free energy of the macrocycle in solution (ΔGrel) to obtain an estimate of the binding affinity (which is given relative to the most favorable mode in Tab. 1).
The ligands were prepared with the script “prepare_ligand4.py” provided with AutoDock. As a consequence the internal degrees of freedom of all ring structures were fixed (as required by the AutoDock program); i.e., the dihedral angles of the macrocycle and the cyclohexyl group (oxime side-chain of 3) remained unchanged. Also, the oxime C=N bond and the amide C-N bond of 3 were kept fixed since we are only interested in the E-configuration of the oxime and the cis/trans-isomers of the amide are chemically identical. The only flexible degrees of freedom were the dihedral angles for the following rotatable bonds: Ar—OH (hydroxy groups on the aromatic ring), =N—O- (oxime), -O—C(H2)- & -(H2)C—C(O)- (both of the oxime side-chain).
In order to reflect the flexibility of the macrocycle, multiple conformations of radicicol and the E-oxime 3 were used for docking; these are the representative structures of the conformational clusters as obtained from parallel tempering, clustering, and minimization (Fig. 3). The cyclohexyl group of the oxime side-chain of 3 shows a chair-like conformation in all cluster representatives; this conformation was kept for the docking.
To account for protein flexibility in approximative way, docking was performed to a number of X-ray structures of HSP90 obtained in the presence of inhibitors (PDB entries 1BGQ, 2BRC, 2CGF, 2IWS, 2IWU, 1ZWH, 2IWX). Target flexibility is an important aspect in drug design,  and cross-docking  is a common strategy to account for the different conformations accessible to the protein.
The protein structures were prepared with the script “prepare_receptor4.py” provided with AutoDock, i.e., the protein structures were kept fixed. Autogrid4 was used to generate the docking grid, centered on the radicicol binding pocket, with default settings except that 60 points (instead of 40) were used in each direction. This was done because of the large size of the ligand (e.g., E-oxime 3).
We used the CHARMM potential energy function  with parameter set 22 . Since the ligands 1 and 3 have not been parametrized for the CHARMM force field so far, we chose the bonding parameters and the partial charges of the ligands from the Merck Molecular Force Field MMFF94 [20, 21], and the van der Waals (VDW) parameters were taken from the CHARMM potential by analogy. All parameters are listed in the “Supplementary Material” Fig. S4 & S5, Tab. S3-S7). A cutoff radius of 12 Å was used for atom based force shifting of electrostatic interactions and potential shifting of VDW interactions.
The docking pose for PDB entry 1BGQ served as the starting point for the setup of the structural model for each binding mode. This structure contained only the protein and the ligand. The crystal waters were added and a cubic box (78 Å) of pre-equilibrated waters was overlaid, the shape of the water box was trimmed to a truncated octahedron. Those overlaid waters that were overlapping with the crystal waters or with the protein were deleted. This resulted in a total of 10856 water molecules; the water shell around the protein measured about 12 Å. We kept the same number of water molecules for all simulations and therefore the crystal waters and overlaid box waters overlapping with the ligand were not deleted but shifted away from the ligand into the bulk water. The shifting magnitude and direction was chosen randomly so that the shifted water molecule was at least 5 Å apart from the ligand and the protein.
Protein and ligand atoms were fixed during the equilibration of the water molecules. We first minimized the waters by 1000 steps of steepest descent (SD) and 1000 steps of adopted-basis Newton-Raphson (ABNR) minimization, and then let them relax by a 500 ps MD simulation at 300 K. The temperature was maintained by the Hoover algorithm, the volume was kept constant at this early stage. The integration time step was 2 fs and SHAKE constraints  were applied to covalent hydrogen bonds of the waters. At the end of the MD run the waters were again minimized by 1000 steps of SD and 1000 steps of ABNR minimization. The constraints on the protein and ligand atoms were replaced by harmonic position restraints with a force constant of 200 kcal mol−1 Å−1. No restraints were set on the hydrogen atoms of the protein or the ligand. The system was deeply minimized by 10000 steps of ABNR minimization, and the force constant of the harmonic restraints was slowly scaled to zero during this minimization.
The ligand and the water molecules were heated from 5 K to 300 K within 60 ps of MD with the the protein kept fixed. We switched from constant volume to constant pressure (Langevin-Piston algorithm) and relax the system for 1 ns. Then the constraints on the protein were replaced by harmonic position restraints with a force constant of 0.05 kcal mol−1 Å−1. No SHAKE constraints nor harmonic position restraints were set on the hydrogen atoms of the protein and of the ligand. The integration time-step was therefore lowered to 1 fs. Within the next 2 ns the force constant of the harmonic restraints was slowly decreased to 0.0025 kcal mol−1 Å−1. It was kept at this value during the following production phase of 500 ps. We saved snapshots every 500 fs during the production phase which yielded ensembles of 1000 structures.
The free energy of binding for a given binding mode relative to the most favorable binding mode, ΔΔGbind, can be estimated by the linear interaction energy approach :
with the ensemble average of the non-polar van der Waals (VDW) interaction energy between the ligand and its surrounding when bound to the protein; the Δ indicates that the difference with respect to the most favorable mode (i.e., the mode with the smallest calculated value of ). The expression stands for the corresponding term of the polar electrostatic interaction energy.
In the LIE approach the ensemble averages of the interaction energies are usually calculated once for the ligand bound to the protein and once for the free ligand in solution. Since we calculate here relative binding free energies between binding modes of the same ligand and same protein, we can neglect the solution contributions because the terms would cancel when calculating ΔΔGbind.
In early papers on the LIE method by Åqvist et al., , Paulsen & Ornstein, , Jones-Hertzog & Jorgensen, , and Wang et al. [36, 37] the coefficient β was fixed at 0.5 since the electrostatic contribution to the hydration energy of a single ion is equal to half the corresponding ion-water interaction energy.  The suggested values for α then ranged from 0.16 [13, 39] to 1.043  based on comparisons of calculated and experimental values for different protein/ligand sets. In more recent works the value of β has been a variable as well; values in the range 0.17 to 1.15 have been used. [40, 41] It has been shown that the values of α and β do only slightly depend on the applied force field;  for the Amber95 force field  (a similar force field to the one used in this work, CHARMM) the optimal agreement was found with α=0.30 and β=0.17. These values are very similar to those that were determined by Ganguly & Mukhopadhyay with the CHARMM force field for the protein ricin B lectin and sugar-type ligands (α=0.36,β=0.16).  In the most recent works by Åqvist and coworkers,  a constant value of 0.18 for α and a range of 0.33-0.50 for the value of β is suggested; β depends on the nature of the ligand .
The optimal values of α and β are not known for E-oxime 3-like ligands bound to HSP90. We tested therefore three different protocols, i.e., three different tuples of α,β that have been recently suggested in the Literature: I) α=0.18,β=0.50 (the coefficients of Åqvist and coworkers with an upper bound of β); [17, 39] II) α=0.18,β=0.33 (same as I but with a lower bound for β); and III) α=0.36,β=0.16 (coefficients obtained for the CHARMM force field) . The VDW contribution increases with respect to the electrostatic contribution from protocol I) to III) so that I) and III) represent two extremes for the relative weighting of the VDW and electrostatic contributions.
The ensemble averages and were calculated from the snapshots of the MD trajectory using the same non-bonded parameters as used in the MD simulation. To test the robustness of the estimation, we calculated the interaction energies also with a much larger cutoff value of 22 Å.
From the MD production phase we selected snapshots every 50 ps along the trajectory. The snapshot structures were minimized by 2000 steps of SD and 5000 steps of ABNR minimization. This yielded 10 structures of every binding mode (see Fig. 2 for the 10 structures of the radicicol binding mode) for which the backbone 15N NMR shifts were calculated.
The NMR calculations are very CPU intensive. One QM/MM-calculation (see below) takes 6-12 hours on 4 CPUs of an IBM SP4 machine; for a single backbone nitrogen nucleus the determination of seven chemical shifts (one for the binding mode of radicicol and six for the binding modes of E-oxime 3) requires more than 1500 hours of CPU time. We limited therefore the calculations to a subset of ten protein residues whose backbone 15N chemical shift shows a significant influence as a function of the ligand, i.e., residues with a large experimental shift difference between radicicol-bound and E-oxime 3-bound protein. These are the six residues with the largest negative shift difference (Δδ = δE−oxime − δradicicol < −1.3: Q9, L89, I90, N92, T95, and I167), the three residues with the largest positive shift difference (Δδ > 0.8: G83, G94, and F156), and residue F124 that is in direct contact with the ligand and that also shows a considerable shift difference for the backbone nitrogen (Δδ = -0.7 ppm).
The molecular system, i.e., protein, ligand, and solvent, is too large to be described entirely by quantum mechanics. Instead a small set of atoms was cut out from the entire molecular systems and the chemical shift for the nucleus in the center of the set was calculated with a QM/MM-Hamiltonian.  We focused on the ten backbone 15N of interest (see above) and for each of them we cut out an individual set of atoms. This means that for a given MD snapshot we performed an independent QM/MM-calculation for each of the ten sets (= ten backbone nitrogen nuclei). In total there were 700 QM/MM-calculations (7 binding modes × 10 snapshots/binding mode × 10 sets/snapshot).
For each residue of interest a separate QM/MM-partitioning was made. The midpoint between the nitrogen and hydrogen served as the center of the quantum-mechanics (QM) region. The inclusion of atoms into the QM region was made based on the following criteria: protein residues, water molecules, and the ligand were completely included if they have at least one atom within 5 Å of the center. The selected atoms that had an open valence due to this selection criterion (i.e., N- and C-terminal atoms of protein residues connecting to protein residues outside of the QM region) were capped by an acetyl group (N-terminus) or a methyl-amino group (C-terminus). To form the cappings we additionally included the following atoms in the QM region. From residues connecting to N-terminal QM atoms: C,O,Cα/Hα atoms; from residues connecting to C-terminal QM atoms: the N,HN,Cα/Hα atoms, respectively. Missing hydrogen atoms in the cappings groups were generated with the HBUILD facility of CHARMM. Remaining residues (or water molecules and the ligand) were included in the molecular-mechanics (MM) region if they had at least one atom within 20 Å of the center of the QM region. From MM residues that connect to QM residues, only those atoms were included in the MM region that were more than 3 Å away from the QM capping atoms. This was done to avoid unrealistic interactions between those atoms. The sum of the charges of these neglected atoms was distributed equally among the included MM atoms of the corresponding residue.
The NMR shielding tensors were calculated by the GIAO method as implemented in the program GAUSSIAN03 . We used the density functional theory (DFT) method B3LYP/3-21G* for the QM region and represented the MM region as point charges applying the command CHARGE of GAUSSIAN03. The small split-valence basis set 3-21G* has been shown to produce meaningful results for small peptide fragments and calculated shift differences are comparable in quality to larger basis sets like 6-311+G(2d,p).
The chemical shifts were calculated by
with σj the isotropic shielding of nuclei j on an absolute scale in ppm and < … > indicates an average of the 10 ensemble structures. The value σref is needed to obtain chemical shifts relative to the indirect experimental reference (see below). We used the backbone nitrogen of residue 90 as an internal reference, which means that we added a constant value so that the calculated value of the backbone nitrogen of residue 90 became identical with its experimentally measured chemical shift. The offset value σref differs for each binding mode, the range is 277.0-278.5 ppm. The value of the offset has no influence on the calculation of the correlation coefficient.
Uniformly 13C- and 15N-labeled yeast HSP90 N-terminal domain (residues 1-210) was expressed and purified as described in Ref. . NMR measurements were performed on a Bruker DMX 750 spectrometer (Bruker, Rheinstetten, Germany) operating at 750 MHz proton frequency. Chemical shifts were measured by 1H,15N-HSQC experiments. Shifts were measured first for the unbound protein (500 μM in 40 mM potassium phosphate, pH = 7.5 at 303K) as a reference and then in presence of a 2:1 excess of the ligands radicicol (1) or E-oxime 3.
Referencing of chemical shifts: The proton chemical shift were referenced directly to DSS (2,2-dimethyl-2-silapentane-5-sulfonic acid); the proton shifts of the methyl group of DSS were set to 0 ppm. The heteronuclear chemical shifts were referenced indirectly using the following equation: [47, 48]
where ν0X is the absolute frequency of 0 ppm for the X spin (with X = 13C or 15N), is the absolute frequency of 0 ppm for the 1H spin, and Ξ is the relative frequency for the X spin compared to the 1H spin (0.251449530 for 13C, 0.101329118 for 15N) .
Assignment of the ligand-bound HSP90 N-terminal domain: The published chemical shift assignments of the free domain [49, 7] were used to identify the chemical shifts of unbound HSP90. The shifts of the ligand-bound protein were assigned by comparison with the chemical shifts of unbound yeast HSP90 and using a set of standard HNCA and HN(CO)CA triple resonance NMR experiments. The assignment procedure was facilitated with the program PASTA. 
By analyzing the equilibrium ensembles of radicicol in solution, we find six low energy clusters within 3 kcal mol−1 (see Tab. 1). The first cluster is more than 50 times more populated than the other clusters taken together. In this global minimum basin the ligand adopts a conformation where the planes of the aromatic ring and the macrocycle span an angle of about 90 degrees (Fig. 3). This conformation corresponds to the L-shape of Ref.  and to the conformation NAMFIS-2 of Ref. ,respectively; it is the bioactive conformation in the cocrystal structure , i.e., the conformation of the ligand when it is bound to HSP90. Beside this L-shape basin, there are two other main conformational basins: clusters where the macrocycle is positioned on the opposite site of the aromatic ring (L’-shape clusters 1c, 1e) and planar-like (i.e., less bent) conformations (P-shape clusters 1b, 1d, 1f). Structures 1e and 1c differ by their orientation of the methyl group at carbon position 17. In 1c the methyl group is equatorial relative to the aromatic ring, in 1e it is aligned axially. The same difference can be noted between the P-shape conformations 1d and 1f. The third P-shape conformation, 1b, features also an equatorial orientation of this methyl group but the carbonyl bond of the ketone points “above” the aromatic ring (“above” indicates the same side as of the macrocycle in the L-shape conformation), whereas in all other structures this bond is equatorial relative to the aromatic ring (1a, 1c, 1e) or it points below the ring (1d, 1f). Structures 1e and 1f of the highest energy basins correspond to the L’- and P-shape structures of Ref.  and the values of the relative stability given there (derived from potential energies) agree reasonably with the free energy values of this work. The highest stability of the L-shape conformation in solution is also confirmed by density functional theory calculations. 
For the E-oxime 3 we find only four clusters within 3 kcal mol−1 of the highest populated cluster (Fig. 3). The two highest populated clusters contain L’-and L-shape conformations. Due to the lack of a chiral center (which contrasts this ligand with respect to radicicol), the L’- and L-shape conformations are mirror images and the population of these enantiomeric conformations should be identical in an achiral environment such as the aqueous solution; the parallel tempering simulation yields a population difference of less than 10 %. The two other clusters are much less populated; they contain each one of the two enantiomeric planar-like (P- and P’-shape) conformations. The difference in stability between the L and P conformation is very similar to that reported for E-oxime 11 of Ref. , which is the same oxime as 3 except that it has the macrocyclic scaffold of radicicol (instead of C17-demethylated pochonin D).
The L-shape conformation 1a docks for all employed X-ray protein structures in the same way as in the X-ray structure of radicicol-bound HSP90, i.e., the aromatic ring sticks into the binding pocket formed by Leu34, Val36, Asn37, Ile77, Asp79, Phe124, Thr171, Leu173; and the macrocycle is bent “upwards” (see Fig. 4) so that the epoxide oxygen interacts with the side-chain of Lys44. This binding mode was also found in two other docking studies. [52, 12]
The binding energy calculated with the AutoDock scoring function for 1a is significantly stronger than for the conformations 1b–f (see Tab. 1). If the free energy of binding is combined with the conformational free energy of the macrocycle (in solution) the binding mode of the L-shape conformation is more than 2.5 kcal mol−1 favored over the other modes. We use this value as a threshold criterion for the evaluation of the binding modes of the E-oxime, i.e., all modes within 2.5 kcal mol−1 of the most favorable mode will be regarded as “likely” modes.
The L-shape conformation 3b yields a binding mode very similar to that of radicicol, this mode is therefore termed L.usual (Fig. 5). The cyclohexyl group of the oxime side-chain partially fills an aliphatic pocket in close vicinity to the ATP-binding pocket, as it was predicted by Moulin et al. for a similar E-oxime compound (see Fig. 4 of Ref. ). For the L’-shape conformation 3a we find two binding modes, L’.outside and L’.upside-down. The latter closely resembles the L.usual binding mode but the L’ conformation is rotated by 180 degrees to obtain an L-like orientation of the ligand. Both P-shape conformations, 3c and 3d, dock outside of the ATP-binding pocket (P.outside, P’.outside1, P’.outside2). The value of the scoring function is very similar for all binding modes; but in combination with the conformational free energy of the macrocycle only three binding modes are found within a margin of 0.3 kcal mol−1, the other modes are above the threshold criterion of 2.5 kcal mol−1. Based on this criterion the modes L.usual, L’.outside and L’.upside-down are “likely” modes.
We test three protocols of the LIE method (I-III, Tab. 2) to estimate the relative binding free energy, ΔΔGbind, for the six modes of the E-oxime 3. With the recommended coefficients by Åqvist and coworkers for the LIE method, [17, 39] i.e., α is 0.18 and β ranges from 0.50 (protocol I) to 0.33 (protocol II), the binding mode L.usual is largely favored over the others; the modes L’.outside, L’.upside-down and P’.outside2 are extremely unfavorable. We obtain the same result for the larger cutoff value of 22 Å. By increasing the weight of the VDW interaction term, α=0.36,β=0.16 (protocol III), the modes P’.outside1 and L’.upside-down become, however, more favorable, even more favorable than the mode L.usual. But modes L’outside and P’.outside2 are still very unfavorable.
For binding free energies estimated by the LIE protocols I)-III) standard errors of 1.0 kcal mol−1 have been reported.  We use twice this value as a threshold criterion, i.e., if a mode is within 2.0 kcal mol−1 of the most favorable mode in at least one of the three protocols we regard it as “likely”; these are the modes L.usual, P.outside, P’.outside1.
The NMR analysis focuses on the backbone 15N chemical shifts for a subset of HSP90 residues (Fig. 2); the nuclei of this set feature significantly different experimental chemical shifts for radicicol-bound and E-oxime 3-bound HSP90. (The complete lists of measured chemical shifts and the 1H,15N-HSQC spectra are given in the “Supplementary Material”, Tab. S1 and Fig. S2)
The calculated and experimental chemical shifts for radicicol-bound HSP90 are in good agreement (see Fig. 6; Tab. S2 of the “Supplementary Material” lists all calculated shifts). The sample correlation coefficient r is 0.81, and linear regression yields a slope of 0.99. The standard error of the correlation coefficient can be approximated by  (with r the correlation coefficient and n the number of data points); it is 0.11. This allows us to set a threshold for the evaluation of the E-oxime binding modes, i.e., all modes with a correlation coefficient higher than 0.7 (= the correlation coefficient of radicicol-bound HSP90, 0.81, minus its standard error, 0.11) are regarded as “likely” binding modes.
A particularly broad distribution of NMR shifts is calculated for residues F156 and I167 (root-mean-square fluctuation > 3.0 ppm, see also horizontal bars in Fig. 6); it is caused by a change of the three-dimensional structure during the molecular dynamics simulation. We note, for example, that the distance d between the center of the F156 phenyl ring and the backbone nitrogen fluctuates during the simulation (3.43 - 3.71 Å) and that the value of the chemical shift correlates with the distance d due to the anisotropy effect of the aromatic moiety (see Fig. 7, A1-A3). This distance fluctuation results in a broad distribution of calculated NMR shifts for the backbone nitrogen of residue F156 on the pico- to nanosecond time-scale of the MD simulation.
The case of I167 is more complex, i.e., we do not find a single geometric criterion that correlates with the chemical shift. We see, however, that there is a reversible interchange between three different hydrogen bond configurations (Fig. 7, B2-B4): 1) the backbone amide proton of I167 forms a hydrogen bridge with a water molecule that then forms a hydrogen bridge with the amide oxygen of the Q145 side-chain (B2); 2) the backbone oxygen of I167 is bridged to the amide protons of the Q145 side-chain by a water molecule (B3; as is the case in the X-ray structure of radicicol-bound HSP90); and 3) the backbone oxygen of I167 forms directly a hydrogen bridge with the side-chain of Q145 (B4). Since the partner of the hydrogen bond and its length influences the chemical shift,  the calculated shielding of the backbone 15N of I167 fluctuates in this dynamic, interchanging hydrogen-bond network.
Interestingly, the HSQC experiment shows a significant line broadening for residues 167-169 (see Fig. 8) indicative of conformational changes on the micro-to millisecond time-scale. The high flexibility on the pico- to nanosecond time-scale, as seen by the molecular dynamics simulation, could be responsible for the initiation of such larger time-scale conformational changes.
The correlation between experiment and calculation for the six binding modes of E-oxime 3 varies significantly (Fig. 9). The binding modes L.usual and L’.outside yield correlation coefficients similar to that of radicicol (0.83 and 0.78), binding modes P’.outside1 and P’.outside2 have higher correlations (0.93 and 0.89) whereas the remaining two modes show correlation coefficients of only 0.46 (L’.upside-down) and 0.54 (P outside). The linear regression yields slopes close to 1 (0.9 to 1.2) for the first four modes and 0.6 for the latter two.
The experimental chemical shifts are consistent with a single binding mode of HSP90-bound ligand 3, i.e., only a single set of chemical shifts is measured for the ligand-bound protein. Rapid exchange between different binding modes can be excluded due to the extremely slow process of ligand unbinding (several hours).  We can therefore exclude binding modes that lack exp.-vs.-calc. correlation. With the threshold value established for radicicol-bound HSP90 we can exclude binding modes L’.upside-down and P’.outside. In other words, based on the correlation between calculated and experimental NMR shifts only the modes L.usual, L’.outside, P’.outside1, and P’.outside2 are “likely” modes.
The activity of a number of ligands that differ from E-oxime 3 by a slight modification has been measured for human HSP90. If we assume that a given ligand has the same binding mode in the yeast and human protein (see below, section “Concluding Discussion”) we can use the activity data to provide a further criterion for evaluating the binding modes of 3.
The size and character of nitrogen substituent of the amide group in E-oxime 3 influences the activity of the human protein HSP90 significantly.  The exchange of the cyclohexyl group by the smaller N-methyl group decreases the affinity for human HSP90 by a factor of 50. On the other hand, a less hydrophobic group like morpholine amide leads to a ten-fold decrease in affinity for HSP90. This suggest that the cyclohexyl group of 3 is involved in important lipophilic interactions.
For modes L.usual, P.outside, P’.outside1, and P’.outside2 the cyclic amino group is in a shallow protein pocket (see also Fig. S3 of the “Supplementary Material”; yellow spheres) with lipophilic character, whereas in the modes L’.upside-down and L’.outside this is not the case. For the former four modes the exchange of the amide group should have a stronger effect on the binding affinity than for the latter two modes.
The E-oximes 2 and 3 differ by the methyl group at the carbon center C17. The affinity of the two ligands for human HSP90 is on the same order of magnitude for both ligands. 
Replacing the hydrogen atom by a methyl group results in a strong steric overlap of the ligand with the protein for modes L’.upside-down, L’.outside, P.outside, and P’.outside, but not for modes L.usual and P’.outside2 (Fig. S3 of the “Supplementary Material”; red spheres). For the latter two modes the binding affinity is expected to change less strongly with an additional methyl group at C17 than for the former four modes.
Based on their agreement with the SAR data, the modes L.usual and P’.outside2 are “likely” modes.
Our combined method of a) searching ligand conformations in solution, and b) docking these conformations to the protein, is capable to find the correct binding mode of radicicol that is bound to yeast HSP90. The bioactive conformation of radicicol is by far the most stable conformation in solution at 300 K and docking of this conformation yields the same, correct binding mode for all studied crystal structures. Docking of other ligand conformations yields much lower values for the scoring functions. The calculated NMR chemical shifts for this binding mode agree well with the experimental values.
When applying the same method to the E-oxime 3 the result is more heterogeneous, e.g., we obtain six modes that have a scoring function within 0.5 kcal mol−1. We need to evaluate the likeliness of those modes by additional methods. We use four different evaluation criteria. Two criteria are based purely on calculations (relative binding affinities from AutoDock, AD, and from linear interaction energy approaches, LIE). A third test compares calculated NMR shifts with experimental data, and the last one is a qualitative criterion that checks the modes against available structure-activity (SAR) data. The results of these evaluations are summarized in Tab. 3. Accordingly, the most likely mode is L.usual; it is the only mode that has a positive evaluation for all four criteria. All other modes have at least two negative evaluations. The L.usual mode is very similar to the one of radicicol.
Interestingly the L.usual binding mode is largely favored over the other modes when the relative free energy of binding is calculated with Åqvist’s parameters for the linear interaction energy method. As noted earlier, a recent study has shown that these parameters are capable of correctly distinguishing possible binding modes obtained from docking in the case of HIV-RT with the force field OPLS-AA.  We conclude that the parameters work also for our system with the CHARMM force field.
Due to the apparent different flexibility of human HSP90 relative to yeast HSP90 (see “Supplementary Material”, Fig. S1 and text), our conclusions regarding the yeast binding mode of E-oxime 3 have to be transferred with caution to the human binding mode. Long-time molecular dynamics simulations indicate that a region close to the ATP-binding site (i.e., amino acids from helix α5 to helix α6)1 shows a considerably different flexibility pattern in yeast than in human HSP90.
For four ligands (ADP, geldanamycin, CT5 2, H64 3) both the human and yeast structure has been solved by X-ray spectroscopy; for H64 no coordinates are available of the region between helices α5 and α6 for the yeast protein. In all four cases the binding mode of the ligand, i.e., conformation of the ligand and its orientation/position with respect to the protein, does not change when comparing the yeast protein with its human analogue. But the protein conformation matches only for the natural ligand ADP (PDB entries: 1AM1, 1AMW for yeast, 1BYQ for human). With geldanamycin (3C11 for yeast, 1YET for human) and CT5 (2BRC for yeast, 2CDD for human) the region between α5 and α6 adopts a different conformation in the human protein with respect to the yeast protein. It is therefore unclear if human HSP90 adopts the exact same conformation as yeast HSP90 when the E-oxime 3 binds to the protein. The binding mode, however, should be similar for both proteins, i.e., L.usual.
We thank Vincent Zoete for helping us in determining force-field parameters for the ligands. M.S. thanks the University of Strasbourg for hosting him as guest assistant professor. A.T. acknowledges support from the Fondation pour la Recherche Medicale (FRM). F.H. and H.K. were supported by the Deutsche Forschungsgemeinschaft (SFB 594 and the Center of Integrated Protein Science Munich, CIPSM). Funding from the Agence Nationale de la Recherche (ANR) is also gratefully acknowledged. Work done at Harvard University was supported by the National Institute of Health (NIH).
1The numbering of the helices corresponds to PDB entries 1BGQ (yeast) and 1OSF (human).
24-[4-(2,3-dihydro-1,4-benzodioxin-6-yl)-3-methyl-1H-pyrazol-5-yl]-6-ethylbenzene-1,3-diol; PDB ligand ID: CT5.
Publisher's Disclaimer: This is a PDF file of an unedited manuscript that has been accepted for publication. As a service to our customers we are providing this early version of the manuscript. The manuscript will undergo copyediting, typesetting, and review of the resulting proof before it is published in its final citable form. Please note that during the production process errors may be discovered which could affect the content, and all legal disclaimers that apply to the journal pertain.