on-line at: www.actabp.pl Computational study of binding of epothilone A to β-tubulin

Understanding the interactions of epothilones with β-tubulin is crucial for computer aided rational design of macrocyclic drugs based on epothilones and epothilone derivatives. Despite numerous structure-activity relationship investigations we still lack substantial knowledge about the binding mode of epothilones and their derivatives to β-tubulin. In this work, we reevaluated the electron crystallography structure of epothilone A/β-tubulin complex (PDB entry 1TVK) and proposed an alternative binding mode of epothilone A to β-tubulin that explains more experimental facts.


InTRoDUcTIon
Epothilones (Fig. 1) are 16-membered macrocyclic lactones that were discovered and first described over 15 years ago.They are extracted from the myxobacteria Sorangium cellulosum as products of their metabolism (Höfle et al., 1993;Gerth et al., 1996).Epothilones demonstrate microtubule-stabilizing properties (Bollag et al., 1995).Their good water solubility, potent anticancer activity exceeding that of taxane family members, and retained activity against MDR cell lines suggest that epothilones could soon replace taxanes in cancer treatment (Wartmann & Altmann, 1996).Several natural and synthetic epothilones are currently undergoing clinical trials (De Jonge & Verweij, 2005).Because of their highly promising features, numerous attempts have been undertaken to unravel the mechanism of action of epothilones.Understanding the binding mode of epothilones to microtubules is crucial in this quest.Many studies have been dedicated to the determination of epothilone conforma-tions in a range of chemical environments to identify the bioactive conformation (Heinz et al., 2005) and to facilitate design of more potent epothilone analogues.One of the first experiments dedicated to the study of epothilone binding to β-tubulin was an NMR experiment conducted in solution (Carlomagno et al., 2003) where the interactions of epothilone molecule with tubulin were unraveled.
The hypothesis of a pharmacophore common to taxanes and epothilones has recently been discarded by Heinz and coworkers (2005) and Nettles and coworkers (2004) because of the different modes of interactions between taxol and epothilone A with β-tubulin.Due to the differences in size, the two compounds exploit the binding cavity in a unique and qualitatively different manner.The epothilone molecule, being substantially smaller than taxol, occupies only part of the taxol binding cavity.While both the taxol oxetane ring and the epothilone molecule interact through hydrogen bonding with a few common residues (Thr274, Arg276), the much larger taxol molecule interacts mainly through van der Waals forces with residues such as Leu215, Leu217, Leu228, Ala231, Ser234, Phe270 and Pro358 (Löwe et al., 2001).
In this work the Nettles' conformer of epothilone A, available from the complex structure of epothilone A/α,β-tubulin dimer (see Fig. 2.), in which Zn 2+ -stabilized tubulin layers form a system mimicking microtubules, despite certain discrepancies, such as antiparallel orientation of protofilaments (Makowski, 1995) was chosen for docking.This conformation (I) has the epoxide ring in the endoorientation, similarly to the one proposed by Wang et al. (1999).Epothilone and β-tubulin form a hydrogen bond between the amine moiety of Arg276 and =O1.The functional groups O3−H, =O5 and O7−H assume an almost parallel orientation, which enables the interaction with Thr274 and Arg282.The nitrogen atom of the thiazole moiety of the epothilone side chain forms an H-bond with His227 (see Fig. 3.).The purpose of our work was to probe computationally the possibility of epothilone A binding to β-tubulin in a different pose and at different locations of the binding pocket than that observed in the Nettles' complex.

MATeRIAls AnD MeThoDs
The structures of the ligand and receptor were built on the basis of the 1TVK complex structure from the Protein Data Bank.The structures were parameterized with Gasteiger atomic partial charges using Chimera software (Pettersen et al., 2004).The docking procedure was performed using Auto Dock software (Morris et al., 1998) with a cubic grid box with a spacing of 0.375 Å and dimensions of 22 Å.Lamarckian genetic algorithm was used for 2 × 10 7 energy evaluations for each of 200 docking trials.To obtain the docked conformation of epothilone A, preliminary clustering at 2 Å, and then final clustering at 0.5 Å were performed.
First, an attempt was made to reconstruct the 1TVK structure by docking the ligand with either restrained or unrestrained torsion angles to a rigid β-tubulin receptor.During the flexible docking procedure the following torsion angles H−O3−C3−C2, H−O7−C7−C6, O16− C15−C17−C18 and C17−C18−C19−N were allowed to change.
Next, we tried to determine the position of the binding using both a flexible ligand and flexible amino-acid residues of His227, Thr274, Arg276 and Arg282.These residues were chosen for two reasons.First, they are known to interact with epothilone, so they had to be rotated in order to eliminate the interactions at the cavity entrance to allow epothilone to penetrate deeper into the binding pocket.Second, these residues have a rather high rotational freedom, so making them rotatable reproduces their behavior in the real biological system.
The docking procedure consisted of ten docking runs, each employing 200 docking trials, with parameters set as in the rigid docking procedure.Seven of the docking runs produced a structure which was considered as an example of a plausible, novel binding mode of epothilone to β-tubulin (conformation IV).
Table 1 contains the binding energies, relative binding energies and hydrogen bond parameters for the obtained structures (conformations I-IV), with conformation I being the crystallographic structure.
In the second part of our calculations, two models of mini receptor derived from 1TVK were built.These models consisted of amino acids common to the two binding modes of epothilone A that were located less than of 5 Å from the epothilone molecule: Glu22, Val23,  Asp26, Glu27, Leu215, Gly223, Asp224, His227, Leu228,  Ala231, Phe270, Ala271, Pro272, Leu273, Thr274,  Arg276, Arg282, Pro358, Arg359, Gly360, Leu361 and  Ser364.In the first mini receptor, all of the amino acids were taken directly from the 1TVK structure, and in the second one the procedure was similar to the previous one with exception of His227, Thr274, Arg276 and Arg282, whose conformations were taken from the Auto Dock output.For evaluation of the binding energy of both models, the rescoring procedure of Fanfrlík et al. (2010) was employed, as the results of their calculations for HIV-1 protease inhibitors yielded results consistent with known experimental facts.
Both mini receptors with bound ligands, without ligands, and free ligands were geometrically optimized in continuum solvent model COSMO (Klamt & Schuurman, 1993) and in vacuo using semiempirical PM6 method (Stewart, 2007), with dispersion and hydrogen bonding corrections (PM6-DH2) (Korth et al., 2010), implemented in MOPAC2009 (Stewart, 2008).Here, our aim was to compare the binding energy of the ligands in both models.
We adopted the following approach based on the thermodynamic cycle of Raha and Merz (2005): (1) ΔG w bind is the estimated binding free energy in water environment, ΔH f w (X) is the heat of formation in aqueous solution, where X denotes the protein-ligand complex, free protein or free ligand.Similarly, ΔH f (X) are the heats of formation in vacuo of the protein-ligand complex, free ligand or free protein.ΔH f complex (X) corresponds to the enthalpy of the protein or ligand molecule in the complex conformation.ΔG int is the in vacuo interaction energy and is calculated as:   where Deformation and desolvation enthalpies are denoted as: where X means protein or ligand molecule.All heats of formation were calculated at 298 K. Entropy changes (ΔS) were calculated in gas phase using AMBER force-field (Cornell et al., 1995).The binding free energy (ΔG w bind ) is thus the sum of the gas phase interaction enthalpy, gas phase interaction entropy and the sum of desolvation and deformation enthalpies of protein and ligand, as all these terms play important roles in determining the binding energy of a ligand.

ResUlTs AnD DIscUssIon
Docking of the rigid ligand to the rigid protein receptor yielded complexes (Fig. 3) fairly similar to the crystallographic structure of the 1TVK complex (see Table 1).The H-bond between His227 and the nitrogen of the epothilone thiazole moiety was reconstructed, as were the hydrogen bonds between O7−H and Arg282 to which both amine moieties of the arginine contributed.Also the H-bonds between Arg276 amine groups and =O1, and O3−H and Thr274 carbonyl group were reproduced.An interaction between the amide moiety of Thr274 and =O5 (conformation II) replaced hydrogen bonds between hydroxyl group of Thr274 and =O5 and O7−H.The mean value of binding energy was -4.5 kcal/mol, and the overall RMS between the original ligand structure and the computed structures ranged between 1.60 and 1.70 Å.
In the flexible structure of the docked ligand (but still the rigid protein, conformation III, see Fig. 3), all the interactions that were reproduced in the fully rigid docking procedure were present, but an additional intramolecular H-bond between O7−H and =O5 was formed.The mean value of the binding energy was lowered to -5.1 kcal/mol, and the RMS between the Nettles' conformer and the computed structures increased slightly to 1.75-1.85Å.The quite accurate reproduction of the experimental procedure by Auto Dock validates the docking procedure used in this work.
The flexible docking (i.e., both the ligand and the amino acid side chains of tubulin flexible) yielded surprising results.The resulting complexes differed qualitatively from the ones generated in the previous dockings, enabling the epothilone molecule to penetrate the binding pocket more deeply and to establish contacts with additional amino acids, thereby lowering substantially the binding free energy of the complex.The mean binding Four different criteria were taken into account during the evaluation of the computed complexes.The first criterion was the binding of the epothilone molecule to the appropriate fragment of the M loop, as hydrogen bonding interaction of epothilone to the Thr274 and Arg282 residues seems to be crucial (Carlomagno et al., 2003).Next, the models with the highest number of hydrogen bonds were considered, as the hydrogen bonding contributes the most to the binding energy, and thus determines the overall stability of the complex.Third, the depth of epothilone molecule penetration was evaluated, since a deeper and tighter fit into the cavity decreases the possibility of solvation of the ligand.Fourth, the consistency of the computed model with the results of numerous SAR experiments was assessed.It should be noted, however, that due to the inaccuracy of the scoring functions implemented in the docking programs, the Auto Dock score was not taken into consideration during the evaluation of the docked complexes.
In the proposed model of the epothilone A/β-tubulin complex (conformation IV, Fig. 4) the position of epothilone differs substantially from that observed in the 1TVK structure, as indicated by the relatively high RMS distance between the superimposed structures.The hydrogen bond with His227 is absent in IV because the distance between the thiazole moiety in the Nettles' structure and in conformation IV exceeds 8 Å, however, the thiazole group of conformation IV is involved in a hydrogen bond with the amide of Gly360.Unlike in the crystallographic structure, the proposed model has the thiazole moiety very close to the protein surface, which, because of the hydrophobic properties of sulfur, could provide better protection against solvation of epothilone.The oxygen-rich region of epothilone (C3-C7 backbone fragment) provides main hydrogen bond interactions and thus forms the "anchor" which contributes the most to the binding energy of the epothilone/β-tubulin complex.The hydrogen bonding pattern in this fragment also differs from that observed in the 1TVK structure.In IV all contacts with Arg276 are lost and those with Thr274 are substantially altered.The O7−H hydroxyl forms hydrogen bonds with Thr274 amide and the carbonyl group of Pro272, =O5 forms a hydrogen bond with the amine moiety of Arg282 and additionally an intramolecular H-bond with O7−H.The O3−H group interacts with both Arg282 amine groups and with the carbonyl group of Gly360.The hydrogen bonds between the hydroxyl group of Thr274 and the carbonyl group of Pro272, and between the carbonyl group of Pro272 and the amine of Arg282 have a stabilizing effect on the interactions in the C3-C7 region (see Fig. 4).Thus, in this disposition the vulnerable hydrophilic fragment of epothilone is buried within the binding cavity and is less solventaccessible than in the original structure, which indicates higher stability of the complex.
The structure of the complex of epothilone conformation IV/β-tubulin also explains a large amount of structure-activity relationship data.
First, our model conforms well with various modifications of the epoxide moiety in the C12-C13 region.The substitution of epoxide by cyclopropanes or cyclobutanes yields epothilone analogues that demonstrate significant antiproliferative activity (Johnson et al., 2000;Nicolaou et al., 2001).The substituted aziridine epothilone analogues also exhibit rather high antiproliferative activity (Regue-iro-Ren et al., 2001).Here, the deoxyepothilones possessing large groups such as CH 2 OC(O)Ph, as reported by Nettles et al. (2004) also act as considerable tubulin polymerization promoters.The flexibly docked model can easily explain these findings.The C12-C13 fragment lies in a location that can accommodate various modifications of this fragment.It is obvious, however, that the presence of a very large substituent in the C12-C13 region will decrease epothilone analogue activity due to the steric clashes with the tubulin residues or with other fragments of the analogue macrocycle, thus preventing it from adopting the binding conformation.
Second, various modifications of the side chain yielded results that are compatible with our model.Nicolaou et al. (2000) found that pyridines possessing a nitrogen atom at position 2 relative to the side chain spacer are one to two orders of magnitude more cytotoxic than 3and 4-analogues.The reason of such activity can clearly be seen in Fig. 4. which shows hydrogen bonding to the Gly360 amide moiety.Other placements of the nitrogen atom prevents the molecule from "anchoring" the glycine.It should be noted, however, that in methylpyridine analogues the position of the methyl substituent is also important.The sutructure-activity relationship data provided by Nicolaou and coworkers clearly shows that the analogues with a methyl group at positions 3-and 4-to the side chain spacer demonstrate extraordinary antiproliferative activity, whereas the analogue with a methyl at the 2-position with respect to the nitrogen atom proved to have moderate potency and the compound with the 2-placement of the methyl group relative to the side chain spacer was virtually inactive (Nicolaou et al., 2000).These results are well explained by our model.Taking the pyridine analogues into consideration, we found that compounds with a methyl group at 3-or 4-positions to the side chain spacer can easily interact with Gly360, and do not produce any steric clashes with the tubulin residues.The epothilone analogue with a methyl group placed at the 2-position to the nitrogen atom exhibits over an order of magnitude lower cytotoxic activity, mainly due to the repulsive interaction with amides of Gly360 and Arg359.In the last considered analogue, the methyl group at the 2-position to the spacer falls in the region of the Asp26 side chain carboxyl group, which produces quite a severe steric clash and, moreover, the very strong repulsive effect of the electron-rich group also contributes to the analogues inability of proper binding to β-tubulin.One of the most prominent analogues of epothilone B is 20-desmethyl-20-methylsulfanyl-Epo B (Nicolaou et al., 2002) that demonstrates outstanding antiproliferative activity.The only difference between the parent molecule and the analogue is the methylsulfanyl instead of methyl group connected to the thiazole moiety.The authors' explanation for such astonishing activity is the higher binding affinity to β-tubulin.This was achieved by introducing the sulfur atom between the thiazole and the methyl group.The sulfur atom lies in favorable position to the side chain of Arg359, therefore enabling strong hydrogen bonding between the arginine and the sulfur atom.The quinoline analogues synthesized by Altmann's group (Altmann et al., 2000) are also cytotoxic agents with great potency.The incorporation of an olefinic spacer into the aromatic ring causes higher rigidity of the system.Such rigidity significantly reduces the side chain conformational space, and thus the analogue adopts the proper conformation more readily.One can also suspect that, due to the high rotational freedom of glutamic acid, the benzene fragment of the quinoline Binding of epothilone A to β-tubulin moiety can fall into the region of interaction of Glu27, and the bound conformation could be additionally stabilized by π-electron interaction between benzene and the deprotonated carboxyl group of Glu27.
Third, the effects of modifications of the C3-C7 fragment can also be explained in the context of our model.The substitution of the O3-H group with a cyano group did not alter substantially the cytotoxicity of epothilone B (Regueiro-Ren et al., 2002).The H-bond with Gly360 backbone is lost, but the -CN group comes into favorable hydrogen bonding contact with the Arg282 side chain.The unnatural (3R-cyano) epimer is about two orders of magnitude less potent.Our model clearly shows that the (3R) configuration severely disrupts the hydrogen bonding network and therefore significantly reduces activity.
Fourth, the structure of the complex of epothilone conformation IV/β-tubulin also explains the very low efficiency of a C6-C8 bridged epothilone analogue recently synthesized (Zhan et al., 2008).Steric clashes between the C6-C8 cyclohexyl bridge, the side chain of Leu273 and the backbone of Thr274 force epothilone to form a high-energy complex with, presumably, weaker hydrogen bonding interactions.
The most important feature of all the described models, including ours, is the great importance of binding to Arg282, and the hydrogen bonding network that includes Pro272, Thr274 and Arg282 which provide an excellent binding place for electro-donor groups.All the three models also occupy virtually the same region of the binding site, and all of them presumably alter the M loop conformation and thus enable the very strong lateral contacts between tubulin dimers and promote polymerization, although they exploit the binding cavity in different manners.
In the rescoring procedure adapted from Fanfrlík and coworkers (2010) the conformation of epothilone A molecule was taken from the Nettles' structure, and the conformation of the mini receptor was taken from the computational model, as those structures had lower heats of formation in aqueous solution and thus were considered "global" minima.While the rescoring function is very reliable in reproducing the relative binding affinity between ligands and can easily discriminate the binding ones from non-binders, a direct comparison of binding free energy obtained by isothermal titration calorimetry and by this method is strictly speaking impossible, due to the difference in construction of the scoring function (Fanfrlík et al., 2010).Nevertheless, such rescoring should provide a plausible ranking of computational models.Therefore, we followed the example of the authors of the method and use the term "score" instead of binding free energy.
The results of the rescoring procedure for the Nettles' model and for ours are shown in Figs. 5 and 6, respectively.
A lower value of the score means a higher binding affinity of the ligand.It can clearly be seen from Fig. 5 and 6 that the enthalpy of interaction in gas phase plays a crucial role as it is the main contributor to the free energy of binding.The other two terms, that is entropy change and deformation and desolvation enthalpies work as opposing force to the binding and contribute slightly less to the total score, but are necessary for proper reproduction of interactions and thus the final ranking of the model.

conclUsIons
The overall location of our modeled epothilone A binding site is in the same part of tubulin as the one determined by Nettles et al. (2004), and both conformers share two out of four interacting amino-acid residues.The similarity to the Nettles' and the Carlomagno's models goes even further, because all the discussed models require electron-donor groups at C3, such as hydroxyl, cyano, or even C2-C3 double bond, as the C3 carbon is biologically an extremely important site (Carlomagno et al., 2003).However, the model of epothilone A binding to β-tubulin proposed here shows several advantages over the experimentally derived structures.
First, the strong hydrogen bonding interactions increase the stability of the complex, lowering the total energy of the system substantially.
Second, the oxygen atoms in the C3-C7 region lie in a more favorable location, less exposed than those in the Nettles' structure, and are protected from solvent accessibility by the C4, C6 and C8 methyl groups, which greatly impedes the possibility of solvation.
Third, the thiazole moiety lies very close to the surface of β-tubulin and the sulfur atom is exposed and lies in the direction of the macrocycle, which should reduce the number of water molecules able to penetrate the cavity because of the strong hydrophobic properties of sulfur.Fourth, the model explains significant amount of structure-activity relationship data up to date, including the most recent ones (Zhan et al., 2008).
Fifth, the results of the calculations using Fanfrlík et al. (2010) approach are in favor of our computational model and due to reliability of the method, they should be considered as rather meaningful.
Due to the above features of the computed epothilone A/β-tubulin complex structure we believe it should be considered as an alternative to the electron crystallography-derived model and used as a reference in computer-aided design of novel epothilone analogues.

Figure 3 .
Figure 3. superposition of rigid and flexible conformations of epothilone A docked to the rigid crystallographic structure of β-tubulin Interactions of epothilone A with β-tubulin: conformation II (rigid, blue) and conformation III (flexible, magenta) with hydrogen bonding pattern.

Figure 4 .
Figure 4. Predicted binding conformation of epothilone A in β-tubulin binding pocket with hydrogen bonding pattern Interactions of predicted binding conformation (conformation IV) with β-tubulin sub-site.

Figure 5 .
Figure 5.Total score and contributing physical terms for the nettles' complex ΔH int -enthalpy of interaction in gas phase, -TΔS -entropy term, ΔG c -sum of deformation and desolvation enthalpies.

Figure 6 .
Figure6.Total score and contributing physical terms for the computed complex ΔH int -enthalpy of interaction in gas phase, -TΔS -entropy term, ΔG c -sum of deformation and desolvation enthalpies for the most stable computed complex (IV).