Calculated conformer energies for organic molecules with multiple polar functionalities are method dependent: Taxol (case study).

BACKGROUND: Molecular mechanics (MM) and quantum chemical (QM) calculations are widely applied and powerful tools for the stereochemical and conformational investigations of molecules. The same methods have been extensively used to probe the conformational profile of Taxol (Figure 1) both in solution and at the beta-tubulin protein binding site.
RESULTS: In the present work, the relative energies of seven conformations of Taxol derived from NMR and X-ray analyses were compared with a set of widely used force fields and semiempirical MO methods coupled to a continuum solvent treatment. The procedures not only diverge significantly in their assessment of relative conformational energies, but none of them provide satisfactory agreement with experiment.
CONCLUSIONS: For Taxol, molecular mechanics and semiempirical QM methods are unable to provide a consistent energetic ranking of side-chain conformations. For similar highly polar organic structures, "energy-free" conformational search methods are advised.

Topological structure of Taxol (paclitaxel).

Background

Conformational and structural analysis of complex organic molecules has been significantly advanced by the development of molecular mechanics schemes parameterized for a wide variety of organic functionalities. For small organic molecules, the Allinger family of programs has served the community very well for many years.[1] One widely used package that incorporates a range of force fields and features of the Allinger protocols, solvation continuum models and conformational searching options is MacroModel.[2] Two studies by Liljefors and colleagues devoted to an evaluation of quantitative aspects of conformational analysis using the MacroModel force fields demonstrate them to perform rather well for a wide range of organic structures with few polar substitutents.[3] A more recent investigation by Halgren on similar structures points out that there is still work to be done to accurately and completely map conformational energy profiles for organic molecules. [4] An area in which molecular mechanics and conformational analysis are critical is in the evaluation of drug candidates and in the molecular design of novel analogs. A case in point is the intense activity around Taxol (paclitaxel), one of the more clinically effective chemotherapeutic agents against a range of otherwise intractable cancers. During the past decade, numerous NMR studies in solution coupled to conformational analysis have led to suggestions that either the polar[5] or nonpolar[6] conformations represent the bioactive one. The recent electron crystallographic 3.7 Å resolution structure of αβ tubulin and Taxol[7] has stimulated a range of suggestions for both the binding mode and the bioactive conformation[8] including a T-Taxol (butterfly) conformer.[9] Many of the proposals arise from molecular mechanics or constrained molecular dynamics conformational searching focused on low energy conformers of the molecule. Our long-standing interest in the binding forms of Taxol, epothilone and analogs[10] has led to a concern for the ability of current computational methodologies to accurately treat the conformational energy manifolds of such molecules either in solution or at a protein binding center.

The seven-conformer dataset

The Taxol molecule is a complex diterpenoid with a conformationally immobile core consisting of the fused A-D rings. Side chains critical for bioactivity of the molecule and analogs are those emanating from the core at C2, C4 and C13. To evaluate the energetic performance of various computational schemes, seven Taxol conformations derived from experimental data with differing torsional angles in the C13 fragment were examined. The diversity of C3' side chain orientations is illustrated in Figure 2, in which the diterpenoid core A-C rings have been superimposed, and amplified with reference to specific torsions in Table 1. The first five are those with the highest estimated populations (4 to 35%) with a ΔG range of 0.0–1.3 kcal/mol from an NMR/NAMFIS analysis in CDC13 solution. [11,12] Structure 1 corresponds to the non-polar conformation observed as the predominant species in CDC13 or CD2Cl2 as depicted in Figure 3. Characteristic of the form is the "hydrophobic collapse" [13] of the benzamido phenyl at C3' and the benzoyl phenyl at C2. The point is illustrated by the short 5.4 Å distance between the centers of the corresponding phenyl rings. Conformer 2 experiences a similar collapse with a somewhat shorter ring-to-ring distance derived from an alternative set of torsions along the C13 side chain (Table 1). With respect to the close approach between pendant hydrophobic centers, 1 and 2 resemble the semisynthetic analog of Taxol, Taxotere (docitaxel). The latter compound differs constitutionally from Taxol in that the C10 acetate becomes an OH and the NHC(=O)Ph benzamido group is replaced by NHC(=O)O-t-Bu. The single crystal X-ray structure[14] and the 2-D NMR (CDCl3) [6d] of Taxotere likewise demonstrate hydrophobic association, in this case between the tert-butyl and the C2 phenyl group. Structure 5 is the polar form frequently found in DMSO-d6/D2O. It too exhibits hydrophobic collapse, but between the Taxol phenyl ring attached directly to C3' and the C2-benzoyl ring (Figure 2). Isomers 3 and 4 are extended rotamers in which the C2 benzoyl phenyl ring is distant from both C3' terminal rings. Rotamer 4 corresponds to the recently proposed bioactive conformation of Taxol bound to β-tubulin.[9] The perspective given as 4' illustrates the "T" relationship between the three phenyl rings of the molecule (Figure 2). The final two Taxol conformers included in our dataset, 6 and 7, appear together in the unit cell of an X-ray crystal structure determination.[15] The latter is a C13 side chain extended structure, while hydrophobically collapsed 6 is very similar to 5. Table 1 complements Figure 3 in providing a selected set of dihedral angles to illustrate explicitly the conformational variation among the seven Taxol torsional isomers.

Figure 2

Seven conformations of Taxol superimposed within the diterpenoid core shown at the right. The diversity of C13 side-chain orientations with respect to the torsionally rigid A-C ring core is illustrated at the left.

Table 1

C13 side chain dihedral angles for Taxol conformations 1–7 used as starting points for the optimization results recorded in Tables 2-5, deg.

	C12-C13-O-C	C13-O-C1'-C2'	O-C1'-C2'-C3'	C1'-C2'-C3'-N	C2'-C3'-N-C	C1-C2-O-C(O)
1	-155	175	76	73	-89	-98
2	-100	-160	105	-48	-65	-88
3	-122	154	78	93	-149	-88
4	-100	-170	73	79	-89	-89
5	-99	-167	94	164	-168	-89
6a	-101	-177	103	179	-155	-86
7a	-104	180	159	176	-117	-86

a Taxol conformations determined in the solid state;15 optimized using AMBER* with all non-terpenoid core dihedral angles frozen; 6 polar; 7 extended.

Figure 3

Seven conformations of Taxol showing distances (Å) between the centroids of the C3' phenyl rings and the C2 terminus. The T-Taxol or butterfly conformation (4 and 4') is pictured in two views; 4' illustrates the "T" relationship among the three aromatic rings.

Results

Five conclusions are immediately evident from Tables 2,3,4. 1) For a given medium (gas, CHCl3, H2O), there is no consistency in energy ranking within the molecular mechanics methods. MM3(96) and MM3 (2000), performing similar to MM3*, offer no exceptions. 2) None of the force-field/ solvation models posit non-polar 1 to be favored in CHCl3, while only MM3* predicts polar 5 and 6 to be most stable in the water continuum model (Table 2). 3) About half of the force field protocols in Table 2 suggest the extended (uncollapsed) conformer 4 to be lowest in energy. 4) Use of 6-31G*-quality ESP charges as recommended for the continuum solvation models [2b,16] leads to conformer 2, a non-polar type conformer distinct from 1, as the uniform global minimum for both gas phase and solution models (Table 3). The next lowest conformer is predicted to be 15–25 kcal/mol higher in energy, while overall conformer ranking between the different force fields generally inharmonious. Thus, incorporation of ESP charges amplifies and alters the energy rankings. 5) AM1 and PM3 calculations in the gas phase vary in suggesting 2, 6 and 7 as most stable (Table 4). While H2O solvation correctly identifies polar 6, the CHCl3 model also favors either 6 or extended 7 rather than 1 or 2.

Table 2

NMR/NAMFIS and X-ray structure determined conformations of Taxol evaluated energetically by six force fields in the gas phase and two solvation continuum models; Relative energies, kcal/mol.a

	MMFF			AMBER*			MM2*			MM3*			MM3(96)	MM3(2000)
	Gas	CHC13	H2O	Gas	CHC13	H2O	Gas	CHC13	H2O	Gas	CHC13	H2O	Gas	Gas
1	0.0	4.2	1.1	4.1	4.4	4.5	4.8	3.1	0.0	5.9	2.0	1.6	4.2	3.5
2	5.8	8.9	7.5	4.0	8.1	8.0	6.6	2.4	5.4	0.0	0.0	2.5	0.0	0.0
3	1.5	6.6	6.6	8.1	2.0	0.0	8.8	6.0	4.7	10.4	9.1	6.1	16.0	15.2
4	0.0	0.0	0.0	2.5	0.0	3.9	0.0	0.0	5.8	2.1	0.6	2.4	4.9	3.6
5	3.6	6.0	1.7	1.2	0.1	0.5	3.5	0.6	3.9	3.3	1.4	0.0	1.5	0.9
6b	2.5	6.9	6.0	0.0	1.4	1.5	4.8	1.9	6.2	5.6	2.7	1.4	1.8	1.2
7b	2.1	5.6	6.3	5.4	5.6	7.3	11.3	5.0	11.9	12.8	8.5	8.3	4.3	6.8

a Each structure was optimized with the indicated force field and the accompanying GBSA solvation model. [2b]b Taxol conformations determined in the solid state;15 optimized using AMBER* with all non-terpenoid core dihedral angles frozen; 6 polar; 7 extended.

Table 3

NMR/NAMFIS and X-ray structure determined conformations of Taxol evaluated energetically by four force fields in the gas phase, two solvation continuum models with the use of scaled ESP atomic charges; Relative energies, kcal/mol.a

	MMFF			AMBER*			MM2*			MM3*
	Gas	CHCl3	H2O	Gas	CHCl3	H2O	Gas	CHCl3	H2O	Gas	CHCl3	H2O
1	60.0	57.8	57.5	76.0	74.8	88.5	47.8	43.2	40.1	46.3	41.9	41.9
2	0.0	0.0	0.0	0.0	0.0	0.0	0.0	0.0	0.0	0.0	0.0	0.0
3	66.1	62.6	61.5	81.6	80.9	86.6	39.2	29.7	25.8	33.2	30.1	30.1
4	22.6	18.0	24.9	21.7	21.5	24.5	24.3	15.2	21.7	28.1	16.4	16.4
5	36.5	31.6	35.3	33.4	30.9	37.8	33.0	28.9	28.5	32.0	27.4	27.4
6b	23.7	23.0	28.1	20.0	24.9	27.2	30.0	22.6	24.6	28.6	22.9	22.9
7b	50.2	49.9	60.6	38.8	44.2	53.6	72.4	63.2	72.3	81.3	68.0	68.0

a Each structure was optimized with the indicated force field and the accompanying GBSA solvation model. [2b]b Taxol conformations determined in the solid state;15 optimized using AMBER* with all non-terpenoid core dihedral angles frozen; 6 polar; 7 extended.

Table 4

NMR/NAMFIS and X-ray structure determined conformations of Taxol evaluated energetically by semiempirical methods; Relative energies, kcal/mol

	AM1a						PM3a
	AM1//AMBER*			AM1//AM1			PM3//AMBER*			PM3//PM3
	gas	CHCl3	H2O	gas	CHCl3	H2O	gas	CHCl3	H2O	gas	CHCl3	H2O
1	3.7	4.0	5.3	4.0	2.0	2.5	4.6	5.8	7.6	0.2	3.8	2.1
2	1.8	4.7	5.1	0.0	0.8	0.2	2.4	6.6	7.7	0.3	6.9	4.3
3	9.3	9.4	11.5	6.7	4.7	5.9	8.5	9.2	11.5	2.6	5.7	4.4
4	1.7	1.4	4.7	0.4	2.1	0.3	2.8	2.8	6.1	2.0	4.4	4.1
5	1.4	1.9	2.2	4.1	2.5	1.8	1.3	2.1	2.4	1.4	4.6	1.3
6	0.0	0.0	0.0	3.1	0.9	0.0	0.0	0.0	0.0	1.2	3.7	0.0
7	2.8	0.3	4.3	4.7	0.0	3.2	3.8	1.4	5.4	0.0	0.0	0.4

a AMSOL, PM3/SM5.4a [reference 18] Solvation energies calculated at AMBER* geometries.

Discussion

Electrostatic interactions dominate and differ across methods

For Taxol, application of commonly applied force fields as well as the semiempirical methods AM1 and PM3 results in an ordering of conformational energies that is method dependent and mostly inconsistent with experimental data. In an attempt to understand this behavior within the force field framework, we examined the various energy contributions for each conformation and method. In most cases the overwhelming factor is the electrostatic term, a component that is seriously amplified by incorporating ESP charges. The point is likewise illustrated by damping rather than magnifying intramolecular electrostatics. Taxol includes ten polar functionalities: five 3-atom units (four esters, one amide) and five 1- or 2-atom units (three OHs, one ether, one C9 carbonyl). In one set of calculations, the latter five groups were converted to the hydrocarbon analogs (CH3, CH2 and C=CH2, respectively). The optimized structures (Table 5) are very similar to those cited in Tables 2 and 3, but the individual global minima are shifted to different conformers. When the five C(=O)-X units are converted to trans butene moieties, now removing all heteroatoms and most of the electrostatics in the molecule, the energy minima shift once again and the average energy spread between high and low unconstrained energy conformations (i.e. 1-5) diminishes to an average of 2.8 kcal/mol (Table 5). The latter is to be compared with an average of 7.5 and 57.6 kcal/mol for the fully optimized structures of Tables 2 and 3, respectively.

Table 5

MM2* energetics of Taxol conformations denuded of polar groups; Relative energies; kcal/mol.a

	Taxol b			Taxol-HClc			Taxol-HC2d			Taxol-HC3e
Conf	Gas	CHCl3	H2O	Gas	CHCl3	H2O	Gas	CHCl3	H2O	Gas	CHCl3	H2O
1	4.8	3.1	0.0	0.1	2.1	0.0	0.5	0.2	0.0	1.4	0.0	1.3
2	6.6	2.4	5.4	0.0	1.2	1.6	1.3	2.6	3.4	2.6	3.3	2.4
3	8.8	6.0	4.7	7.4	6.1	6.7	1.8	2.4	2.5	1.7	0.1	2.0
4	0.0	0.0	5.8	1.8	2.0	0.5	0.0	0.0	0.4	0.0	0.0	0.0
5	3.5	0.6	3.9	0.6	0.0	0.3	1.3	1.0	1.7	0.6	0.0	1.0
6f	4.8	1.9	6.2	3.6	3.2	3.4	6.2	5.4	7.2	9.1	8.2	9.5
7f	11.3	5.0	11.9	7.5	4.2	6.5	9.1	6.7	11.8	8.6	7.2	9.9
ΔEg	8.8	6.0	5.8	7.4	6.1	6.7	1.8	2.6	3.4	2.6	3.3	2.4

a The number of low quality MM2* parameters for the hydrocarbon (HC) analogs are very few (HC1 1.0, HC2 1.0, HC3 0.0%). Thus, the energetic changes are primarily electrostatic in origin. b MM2* relative energies as presented in Table 1. c Cl, C7 and C2' OHs in 1 were replaced with CH3; C5 ether with CH, and C9=O with C=CH2. d C2, C4, C10 and C13 esters and C3' amide were replaced with trans-CH=CH. e All ten polar groups replaced with hydrocarbon as in b and c. f Taxol conformations determined in the solid state13; optimized using AMBER* with all non-terpenoid core dihedral angles frozen; 6 polar; 7 extended. g The energy spread (kcal/mol) between the highest and lowest unconstrained energy conformations; i.e. 1-5.

The molecular mechanics results might be ascribed to the presence of "low quality"[17] parameters in the force fields employed. Indeed MM3* and MMFF involve 74 and 62 such interactions out of a total of 439 and 555 (17 and 11%), respectively. However, MM2* and AMBER* incorporate only 2 and 6% such parameters, respectively, yet their ability to accurately predict the purported experimental forms is not improved. The two semiempirical methods each locate 2–6 conformations with energies no more than 2.5 kcal/mol above the lowest, but none of the calculations selects the same set of low energy conformers.

Modern force fields with stretch, bend, torsional, hydrogen bond and cross-terms are parameterized primarily within three-bond units (A-B-X-Y). Of the two additional key terms, the van der Waals component operates with a very short range force and essentially prevents atom-atom interpenetration. In the absence of damping, the single long-range interaction that persists from one end of the molecule to the other is the electrostatic contribution. While individual bond dipoles are characteristic of the Allinger MM family of force fields, other molecular mechanics methods employ atomic charges. In either case, no general scheme for parameterizing the electrostatic interaction among multiple polar groups in three-space in terms of conformational energies has yet emerged.

Similarly, within a quantum chemical semiempirical framework, charge polarization notwithstanding, conformational energies are sensitive to method and charge distribution. Inspection of Table 4 reveals that, unlike the force-field results described in Tables 2 and 3, all four semiempirical recipes predict polar conformers 5 and 6 to be the low and "global" minimum energy forms, respectively, when the aqueous solvation continuum model is employed. The one exception, AM1//AM1 suggests both polar and nonpolar forms 6 and 2 to be equally populated. This limited success is negated, however, by the incorrect prediction that the same polar species are also the low energy forms in the gas phase and in the CHCl3 solvation model. In no instance do the semiempirical calculations predict either nonpolar forms 1 or 2 to be dominant in CHCl3. As mentioned above, the AM1//AM1 calculations posit equal populations for 2 and 6 within the latter regime, but they simultaneously relegate nonpolar 1 to a higher energy. The remaining semiempirical methods predict the empirically verified nonplar conformers 1 and 2 to be unobservable in the chlorocarbon solvent with energies ranging from 3.8 to 6.9 kcal/mol above the polar collapsed conformation. It might be concluded that the water model is robust, while the chloroform model is ill-parameterized, although the AMSOL literature on solvation[18] gives little basis for this contention. Alternatively, the few inerrant predictions for 6, combining both structure and continuum aqueous solvation energies, might well be fortuitous. The likelihood that this is correct is accentuated by the fact that Taxol conformational analysis by NMR has never been performed in pure water, but always as a mixture of DMSO-d6 and D2O.[5] In addition, in anhydrous methanol, a solvent with a high dielectric and the capacity for explicit hydrogen bonding as in water, the polar form goes undetected. [5b]

T-Taxol (the butterfly conformation)

For the unsuspecting organic or medicinal chemist, however, either interpretation is equally unfortunate. Faced with a large, poly-polar organic structure and a solvation-equipped semiempirical package, the user is limited by the proposition that the aqueous solvation model might provide the correct global minimum, but that other conformers may or may not be assigned an appropriate relative energy. In less polar chlorocarbon solvents, the Taxol structure and other molecules with similar complexity do not appear to be capable of even qualitative ranking by the various methods with respect to energy. An interesting but fortuitous outcome in this respect is that conformer 4, the T or butterfly conformer of Taxol, is predicted to be the dominant conformation by MMFF, AMBER* and MM2* under various protocols (Table 2). Semiempirical models likewise predict it to be of low energy, though not the most stable (Table 4). This conformer was assigned a low population in the multi-conformational analysis of Taxol in chloroform[12], used in the fitting of the electron crystallographic density of the β-tubulin/Taxol complex in zinc-stabilized sheets, and ultimately assessed as the binding conformation for Taxol in this structure. [9]. A unique, but difficulty observed conformer, computational methodology overestimates its stability and might have led to its discovery by misrepresentation. While the operation of serendipity is much to be desired in the drug seeking process, errant methodology is the least desired path to discovery.

Conclusions

As a consequence of the above considerations, for highly polar molecules like Taxol it is overoptimisitic to expect that conformational searching based on a computational energy criterion will yield results faithful to experiment. The use of solvent continuum models clearly does not ameliorate the situation. Our observations and those of Halgren[4] suggest that a given force field or quantum chemical method selected for potential application to a highly-polar molecular system be carefully validated to assure that conformational energy comparisons are not spurious. Until this problem can be addressed at a fundamental theoretical level, it would seem prudent to employ an "energy-free" conformational search protocol. We find the NAMFIS approach valuable in this respect;[11,12,19] however there are a number of equally attractive alternatives[20] that, in principle, can also bypass the energy catastrophe. Future work will be devoted to examining a range of functionalized organic molecules in order to define the boundary of molecular polarity within which standard conformational search methodology can be applied with confidence.

Materials and Methods

Force fields and semiempirical methods applied to Taxol

The study was conducted by employing four popular force fields in MacroModel6.5 (MM2*, MM3*, AMBER* and MMFF).[2] Each was used in the "gas phase" with two continuum solvation models (CHCl3, H2O) [2b] with default MacroModel atomic charges (Table 2). To provide a comparison with MacroModel, we also tested solvated MM3(96) and polarization enhanced MM3(2000).[21] In addition, the seven AMBER* conformations were fitted with MNDO electrostatic potential (ESP) charges scaled to the 6-31G* level[22,23] and the calculations repeated (Table 3). The structures were also fully optimized with the AM1 and PM3 semiempirical methods[22] and supplemented by AMSOL SM5.4a solvation energies (Table 4). [24]