The starting structure of epothilone B was taken from the crystallographic coordinates that are deposited (1) in the Cambridge Structural Database (Refcode: TIPFON). The crystal structure was optimized with the AM1 and PM3 semiempirical Hamiltonians using the MNDO94 program. Of these two Hamiltonians, the PM3 model was found to converge closest to the crystal structure geometry (RMS deviation = 0.087 Å for PM3, and 0.575 Å for AM1, for all heavy atoms in the 16 atom membered epothilone ring). Molecular mechanics potentials from the CFF91 force field (2) were selected to produce an optimized geometry (until the norm of the gradient was 1.0 x 10¯³ kcal/mole) as close as possible to the crystal structure and the PM3 model. The resulting overlap of all molecular models are shown in Figure 1. The RMS deviation of the CFF91model to the crystal structure was 0.115 Å. Figure 2. displays the potentials that yielded the close geometry from molecular mechanics optimization.
The epothilones are molecules that are inherently flexible, have hydrogen bond donors and acceptors, as well as other interesting features. Our objective was to sample a sufficient amount of conformational space in order to determine the possibility that epothilone B utilized the preorganization of certain key molecular substituents to bind to the taxane-tubulin site. We chose high temperature molecular dynamics to conformationally sample epothilone B.
An initial pilot study found that 1000° Kelvin would not cause inversions of the chiral centers of the epothilone scaffold using the CFF91 potentials that were described above. Therefore, a molecular dynamics simulation temperature of 1000° Kelvin was chosen for the conformational sampling phase. For all molecular dynamics and minimization calculations, the entire functional form of the CFF91 potential energy expression was enabled. Each molecular dynamics simulation started with the minimized crystal structure of epothilone B, and a new random number seed. Temperatures were increased slowly, from 0° by 25° increments for 1 ps, until the target simulation temperature was reached (1000° K). All simulations used a time step of 0.5 femtoseconds, and were initiated with an equilibration phase using direct velocity scaling for 10 picoseconds. This phase was followed by 100 picoseconds of molecular dynamics using Berendsen method of temperature-bath coupling. During this phase of dynamics, a sample structure was collected every 100 femtoseconds to be minimized with conjugate gradients until the norm of gradient was 1.0 x 10¯³ kcal/mole.
Since the epothilone ring structure is capable of forming a variety of conformations that are stabilized by intramolecular hydrogen bonds, we employed a variety of dielectric conditions in order to screen the electrostatics thereby enabling less bias towards intramolecular hydrogen bonding. The general molecular dynamics simulation described above was repeated using the following dielectric conditions: in vacuo (1.0), distance dependent (1.0 x r), water (80.0), dimethyl sulfoxide (46.7), methanol (32.8), ethanol (24.3), methyl ethyl ketone (15.8), chlorobenzene (5.708), iso propyl ether (3.38) and n-decane (1.991). This resulted in a total of 10,000 energy minimized structures (1000 minimized conformations per dielectric condition). Generally, from each dielectric condition (N=1000), 100 of the lowest energy conformations usually fell within 15 kcal/mole of energy, wheras 50 of the lowest energy conformations fell within 5 kcal/mole of energy between conformers. On average, 35 of these 50 isomers were conformationally unique.
Computed distances between the epoxide O and other atoms of epothilone B revealed a consistent pattern. It was observed that the centroid of the 4-gem-dimethyl group, the 7-OH, and 3-OH, fell at a mean distance of about 6.5 Å (± 1.25 Å) from the epoxide O within a dielectric condition (for 35 unique, lowest energy conformers, results from the distance-dependent dielectric condition given as an example). When the same measurement was made from 10 unique, lowest energy structures across all dielectric conditions (N = 100), the mean distance between the pseudoatom and the epoxide oxygen was 6.93 Å with a standard deviation of ±1.18Å.
This pharmacophoric similarity between the baccatin ring system's 16-gem-dimethyl group, 9-carbonyl, 1-OH, oxetane O and epothilone's 4-gem-dimethyl group, 7-OH, 3-OH, epoxide O led to the formation of a more complete overlap model. This was accomplished by template forcing the remaining hydrocarbon portions of the epothilone B's 16 atom ring system onto the corresponding groups of the baccatin core of docetaxel. Eleven, template-forced minimizations, employing a force of 100 Kcal/mole/Ų , were performed before the most complete overlap was found. The final model of the overlapped minima of docetaxel, epothilone B as well as paclitaxel are presented in Figure 3.
The depostited tubulin structure (3) (Refcode: 1TUB) was used as a starting model for the formation of the epothilone B/tubulin molecular model. First, using tethered minimizations, an energy-refined, full atomic structure model of 1TUB was created. The objective was to find the nearest local minimum to the empirically determined coordinates (1TUB). Our strategy entailed applying a 2000 Kcal/mole/Ų of force that was stepped off the structure in 100 Kcal/mole decrements by minimizing with conjugate gradients, until the norm of the gradient was 1.0 x 10¯³. This process was repeated until all applied external force was removed. The resulting model had a total RMS deviation of 1.97 Å across all heavy atoms from the reported tubulin crystal structure coordinates. The conformation of the epothilone B model from the overlap model described above, was docked into the binding site formerly occupied by docetaxel. Once Van der Waals violations of 0.25 Å were removed, the ligand/protein structure coordinates were tether minimized in the same manner as described above. The RMS deviation between the docked epothilone B and the conformation extracted from the overlap model was only 0.71 Å. Optimization of epothilone B conformation from the template-forced model with the semiempirical PM3 method also resulted in close geometries (RMS deviation between the two 16 atom rings is 0.331 Å when superimposed). In the docetaxel/tubulin binding site, epothilone B was found to dock to the protein in two different binding modes. In binding mode I, epothilone B's thiazole overlaps with the C2 benzoate of the baccatin ring system. The proposed overlap with paclitaxel for specific atoms of epothilone B for binding model I is presented in Figure 4. In binding mode II, the thiazole of epothilone B overlaps with a phenyl group of the taxane side chain subsituted at C13 of the baccatin ring system. The proposed overlap with paclitaxel for specific atoms of epothilone B for binding mode II is presented in Figure 5. The overlap of the docked epothilone conformations from binding mode I and II with docetaxel in the tubulin binding site demonstrate complete agreement with our pharmacophore model. Of the two, binding mode II is presently believed to be the most preferable according to the latest structure activity relationship data in the mutant cancer cells.
Conformational analysis revealed two major ring puckers for the sarcodictyin A analog. Of these, only conformation I fits the pharmacophore and the refined docetaxel/tubulin receptor model. The methyl imidazole side chain overlaps in a similar manner to binding mode II of the epothilones. Docking of this conformer in the model revealed an unfavorable but tolerable hydrophobic-polar interaction (dashed line in light blue). However, when threonine-274 is mutated to an isoleucine, this same interaction becomes a highly favorable hydrophobic-hydrophobic interaction. Similar to the epothilones, sarcodictyins are less bulky than docetaxel in the binding site.
The docetaxel/tubulin coordinates indicate that some of the polar groups of docetaxel are not directly bound to the tubulin protein via hydrogen bonding. Using the optimized model described above, single water molecules were built into the model and re-minimized, until docetaxel was completely ligated by water to the tubulin structure. Waters were also included into the region occupied by ß-R282 to form a contiguously hydrogen-bonded solvent network. These results of this model are depicted in Figure 6.
For additional information, answers to questions, or comments, please feel free to contact Rick Gussio.
1. G.Hofle, N.Bedorf, H.Steinmetz, D.Schomburg, K.Gerth, H.Reichenbach. Angew. Chem., Int. Ed. Engl., 35, 1567,1996
2. Maple, J.R., Hwang, M.-J., Jalkanen, K.J., Stockfisch, T.P., Hagler, A.T. Derivation of Class II Force Fields:V. Quantum Force Field for Amides, Peptides, and Related Compunds. J. Comp. Chem., 19, 430-458, 1998. (cff91 force field)
3. E. Nogales, S.G.Wolf, K.H.Downing. Structure of the Alpha Beta Tubulin Dimer by Electron Crystallography. Nature, 391, 199-203, 1998.
