- Email: [email protected]
Molecular mechanics conformational analysis of tylosin
J o u r n a l of
MOLECULAR STRUCTURE ELSEVIER
Journal of Molecular Structure 440 (1998) 121- 130
Molecular mechanics conformational analysis of tylosin P e t k o M. I v a n o v ~a,b ~lnstitute of Organic Chemistry. with Centre of Phytochemistry, Bulgarian Academy of Sciences, ul. Acad. G. Bonchev, bloc 9, 1040 Sofia, Bulgaria bDepartament de Quimica, Facultat de Cikncies, Universitat Autbnoma de Barcelona, 08193 Bellaterra, Barcelona, Spain
Received 17 March 1997; revised 16 May 1997; accepted 16 May 1997
Abstract
The conformations of the 16-membered macrolide antibiotic tylosin were studied with molecular mechanics (AMBER* force field) including modelling of the effect of the solvent on the conformational preferences (GB/SA). A Monte Carlo conformational search procedure was used for finding the most probable low-energy conformations. The present study provides complementary data to recently reported analysis of the conformations of tylosin based on NMR techniques. A search for the low-energy conformations of protynolide, a 16-membered lactone containing the same aglycone as tylosin, was also carried out, and the results were compared with the observed conformation in the crystal as well as with the most probable conformations of the macrocyclic ring of tylosin. The dependence of the results on force field was also studied by utilizing the MM3 force field. Some particular conformations were computed with the semiempirical molecular orbital methods AMI and PM3. © 1998 Elsevier Science B.V. Keywords: tylosin; protynolide; conformation; molecular mechanics; AMBER*; MM3; MC search
1. I n t r o d u c t i o n
An assignment of the IH and 13C NMR spectra of tylosin (1, Fig. 1) in CDCI 3 was recently reported [1], and conclusions had been also drawn regarding the solution conformation of the macrocyclic ring of 1 in aprotic solvents, based on vicinal proton-proton coupling constants and NOE enhancements. The authors further supported their experimental data with results from molecular mechanics calculations of protynolide (2), a 16-membered lactone containing the same aglycone as tylosin [1]. Two low-energy conformations of 2 were found in that study and i E-mail: [email protected] 2 Recent references to papers related to techniquesfor the search of low-energyconformationsare listed in [5].
discussed as possible candidates for the preferred conformations of 1 in solution. The predominance of the conformation corresponding to the conformation of 2 in the crystal [2] is proposed (conformation 2t.1 in Fig. 4). Not always, however, does the congruence between averaged computed quantities and measured observables provide a guarantee for assessing the conformational ordering, and different computational methods can be at variance in the details of the computed conformational distribution [3]. The results from the computational conformational analysis of a macrocyclic system of the size of tylosin could happen to depend significantly on the molecular mechanics model used [4], as well as on the way the search for the most probable conformations is carried out. 2 The author recently stressed the importance of
0022-2860/98/$19.00 © 1998 Elsevier Science B.V. All rights reserved PII S0022-2860(97)00233-0
122
P.M. lvanov/Journal of Molecular Structure 440 (1998) 121-130
IIt° 6
"1
lvle ~ . . . . . ~ ' 1 1
12 1112
"' Me ~5'~...---O
HO
23
-4 /]J 13
t.,,
,9
7 ~,.
18
661
~le
~
0
\
20
t CH2CHO ........
,
1'
a
0
5'
/X---o
-
6'
...
4'
..
a o ~ ' ~ OMe
16 [ / ~ ~ " . . . . . . , ,
OH
171
/ N Me 7'
2,, ~ Me 8'
~ u
~T~ 7" ~ ) H 5"\ ~1" Me 6 t*
Fig. 1. Numbering of atoms in tylosin.
s e l e c t i n g the p r o p e r m o l e c u l a r m e c h a n i c s m o d e l w h e n i n c l u s i o n c o m p l e x e s are studied [4]. T h e s a m e r e m a r k s are v a l i d also for the case o f a large m o l e c u l e w h e r e d i s t a n t f r a g m e n t s o f o n e a n d the s a m e m o l e c u l e h a v e the f r e e d o m to i n t e r a c t in a w a y as if t h e y belonged to different molecules. 3 Molecular m e c h a n i c s m o d e l s h a v e b e e n d e v e l o p e d m o s t l y to m a t c h specific p u r p o s e s ; 4 thus, it s h o u l d n o t b e s u r p r i s i n g if t h e y d i s p l a y a d i v e r s e b a l a n c e o f intera c t i o n s in s o m e cases. T h e m o l e c u l a r m e c h a n i c s conformational analysis of tylosin presents such an example. T h e p r e s e n t c o m p u t a t i o n a l s t u d y is a i m e d at finding low-energy conformations of 1 and 2 u s i n g the o n l y a v a i l a b l e (to us) s t a n d a r d s e a r c h t e c h n i q u e [15] i m p l e m e n t e d in the M a c r o M o d e l and BatchMin V4.0 molecular modelling p r o g r a m s [16]. T h e f o l l o w i n g q u e s t i o n s c a n be formulated that deserve elucidation in this c o m p u t a t i o n a l study: 1. D o 1 a n d 2 h a v e the s a m e o r d e r i n g o f c o n f o r m a t i o n s d i s p l a y e d b y the g e o m e t r i e s o f the 16-membered ring? 2. A r e t h e r e o n l y o n e or t w o l o w - e n e r g y c o n f o r m a t i o n s o f 2, as r e p o r t e d b y o t h e r a u t h o r s [1], a n d w h a t are the f o r c e s w h i c h r e s t r i c t the p o t e n t i a l for flexibility o f the 1 6 - m e m b e r e d ring to o n l y o n e or t w o s t r o n g l y f a v o u r e d forms? 3. Is the crystal c o n f o r m a t i o n o f 2 a m o n g the m o s t s t a b l e f o r m s o f the 1 6 - m e m b e r e d r i n g s o f 1 a n d 2 in s o l u t i o n ?
3Crystal calculations using molecular mechanics are an appropriate test for the ability of a force field to handle intermolecular and distant intramolecular van der Waals interactions. Even tiny modifications of the van der Waals parameters can affect significantly the balance of nonbonded interactions in large molecules or crystals. A comparison between three generic parametrizations, MM2 [6], MM2' [7] and MM3 [8], carried out with the KESSHOU program [9] (see below), can serve to illustrate the problem [10]. Calculations of the n-hexane crystal using the MM3 van der Waals potentials for atom-atom interactions closely reproduce observed unit cell constants. This is not unexpected, because crystal structural data had been extensively used for the parametrization of the MM3 force field [8]. MM2 produces about 15.0 ~3 smaller unit cell volume than the one found by experiment [11]. The MM2' parametrization [7], which was specifically developed with the intention of improving the performance of MM2 in reproducing rotational barriers, gives even larger deviations from experimental crystal data; i.e., the balance of nonbonded interactions has been sacrificed in MM2' for the purpose of improving performance in reproducing rotational barriers. It is well known that MM2 enhances the stabilization of the conformations when van der Waals interactions form the prevailing contribution to the energy balance of intramolecular interactions [3,12]; thus, using MM2 for studying the conformations of protonolide 2 was not the best choice the authors could have made [1]. Besides, modelling of hydrogen bonding was also not included [1]. KESSHOU [9] was recently further developed to handle also electrostatic interactions [9c]. MM3 crystal calculations with full molecular geometry optimization, using CHELPG electrostatic potential derived atomic charges [13], gave perfect reproduction of the packing of the 1.4dicyanobenzene molecules in the crystal ttriclinic crystal lattice, space group PL Z = 1, crystal radius 20.0 A, 373 unit cells considered; full geometry optimization of the reference molecule). The packing of the C60 molecules in the low-temperature phase of crystalline C~ is also exactly reproduced (cubic, Pa3, Z = 4, crystal radius 50.0 A, 332 unit cells). These results certify the improved quality of the MM3 van der Waals potential functions and parameters for modelling nonbonded atom-atom interactions. 4The present status of molecular mechanics force field development was recently reviewed by the author [14].
123
P.M. lvanov/Journal of Molecular Structure 440 (1998) 121-130
2. Computational details The Monte Carlo conformational searcher [15,16] in BatchMin was used (MCMM command) with the AMBER* force field as a model, and the parameters for CHC13 were adopted for modelling the solvent with Still's GB/SA method [17]. Extended cutoff distances for the nonbonded (van der Waals and electrostatic) interactions were set to 20.0 ,~. For protynolide, 2, the maximum number of MC steps was set to 6000 using the global search mode as a parameter for the MCMM command. The minimum and maximum number of degrees of freedom changed in an MC step were set to 2 and 11, respectively. The MC starting structure selection command (MCSS) was executed with use-directed mode (argl = 2) and an energy window of 25.0kJ tool <. Optimized structures with relative energies less than 25.0 kJ mol -j were kept during the search (the DEMX command), while structures with relative energy (after 200 minimization iterations) higher than 50.0 kJ mol -~ were excluded from further energy minimization. All 28 non-hydrogen atoms were used for comparison when eliminating duplicate minima. Twelve endocyclic torsion angles were varied. The four torsional angles of the C = C - C = C - C O fragment were not included in the list of variable torsion angles. The default values for a minimum and maximum dihedral angle variation were adopted. Conjugate gradients, followed by Newton-Raphson minimization, were applied in the final calculations of the 198 unique structures stored after the MCMM search was completed. We followed the same computational scheme for 1. Twenty-one torsional angles were assumed to vary in the case of tylosin. The number of heavy atoms used for comparison when eliminating duplicate minima is 52. These comprise the atoms of the 16-membered ring and the most important atoms from the side chains. The list of 21 variable torsional angles includes the same dihedrals as in the case of the protynolide (2), plus all dihedrals defining the rotational isomerism about side chains bonds. Endocyclic dihedrals of the six-membered rings were not included in the list of variable torsional angles; thus ring closure conditions were imposed only for the 16-membered ring. The starting conformation had the observed conformation of the
macrocycle of 2 in the crystal. The maximum number of MC steps was set to 10 000 using the global search mode as a parameter for the MCMM command. The search ended with 494 unique structures stored. Only three of them satisfied the criteria for good convergence. Thus, all stored structures were further subjected to consecutive conjugate gradients and Newton-Raphson optimizations. Restrictions were not imposed to keep the starting trans configurations around the C ( 1 0 ) = C ( l l ) and C(12)=C(13) double bonds during the MCMM search. Optimized low-energy structures were obtained for the 12-trans and 12-cis isomers of 2 and are designated as 2t and 2c, respectively. All low-energy structures obtained for 1 have cis configuration around the C(12)=C(13) bond. The global minimum conformation of 12-trans-1 is with 19.0 kJ mol -~ higher energy than the global minimum conformation of 12-cis-1. Only the conformations of the 12-cis isomer of 1 will be discussed. These conformations are designated lc.
3. Results and discussion Computed conformational data for lc are summarized in Table 1. Structural and energy data of 2t and 2c are presented in Table 2. These are relative energies of the low-energy conformations obtained with the AMBER* force field, as well as a letter-code representation of important dihedral angles. The numbering of the atoms of tylosin is displayed in Fig. 1. For the sake of compactness of the presentation we used a letter-code designations of dihedral angles, as given in Fig. 2. A letter G designates a dihedral angle in the range from 20 ° to 100 °, G from 100° to 150 °, A from 150° to - 150 °, etc. .20°
-1~0
20°
--
~50°
Fig. 2. Designation of dihedral angles.
124
P.M. lvanov/Journal of Molecular Structure 440 (1998) 121-130
Table 1 Computed conformational data of the 12-cis isomer (le) of tylosin (1)" 1-2-3-4-5-6-7-8-9-
I0-11-12-13-14-15-0-1
G' G' A G ' G' A G ' G A A ~ ' G' A G' G' A G' G A A ~ ' G' A G' G' A G' G A A ~ ' G' A G' G' A G' G A A ~ ' G' A G' G' A G' G A A ~ ' G' g G' G' g G' G g g ~G'AG'G'GG'AAAG'E~A ~ ' G' A G' G' A G' G A A ~ ' G' A G' G' A G' G A A
I
a
n
A A A A A g
E G G' EGG' EGG' EGG' EGG' E~G'
G' G' G' G' G' G' G A E ~ G ' G' A E ~ G ' G' rr
A A A A A A A A A
Xi
X2
X3
X4
X5
G
G G G G G G G G G'
G' ~' ~' ~ G' ~' G' ~' ~'
A A A G A G g A A
G G' ~ A G G ~' ~ G
~' G ~' G' m
AE (kJ mol )) 0.0 4.9 7.7 7.8 8.9 10.0 10.7 I1.9 12.0
Conformation
le.1 le.2 lc.3 lc.4 lc.5 lc.6 le.7 lc.8 to.9
.
X~=3-4-O-1 X2=3-4-O-H X3=3 - 2 - O - M e X 4 = 1 5 - 1 4 - 2 3 - O Xs = 1 4 - 2 3 - O - l " . Other side chain dihedral angles have the same values in all conformations: 2 - 3 - O - H (G'), 4 - 5 - O - 1' (G'), 5 - 6 - 1 9 - 2 0 (G'), 1 4 - 1 5 - 1 6 - 1 7 (A), 6 - 1 9 - C ( 2 0 ) = O (E), 5 - O - 1' - 2 ' (A), l ' - 2 ' - O - H (G), 2 ' - 3 ' - N - M e ( 7 ' ) (G), 4"-3"Y/-O-H (G'), 2 3 - O - 1 ' " - O (G'L 4 " - 3 " - O - M e (G'), 3 " - 4 " - O - H (G').
Table 3 contains dihedral angles and calculated vicinal coupling constants for some pairs of protons in the macrocyclic ring of low-energy conformations of lc, 2t and 2c. Relative energies of the most probable conformations estimated with different molecular mechanics force fields and semiempirical molecular orbital methods are given in Table 4. Stereographic views [18] of some conformations of lc are displayed in Fig. 3. Analogous presentations of conformations of 2t and 2e are given in Fig. 4. 3.1. Protynolide (2)
The computed lowest-energy structures of 2t and 2c are practically of the same energy. Except for the configuration around the C(12)~C(13) double bond, these two structures differ mostly in the dihedral angle C(7)-C(8)-C(9)-C(10), -59.0 ° and 104.0 ° in 2t.1 and 2c.1, respectively. Six-membered hydrogen bonded rings are formed in 2t.l and 2c.1 and these include the C(3) hydroxyl group and the C(1)=O oxygen. The more favourable torsional energy contribution (AEtors = 7,1 kJ mol -I) for the trans isomer, 2t.1, is balanced by advantageous van der Waals and electrostatic interactions in the cis form, 2c.1. The hydrogen bond distances O - H . . . O = C in 2t.1 and 2e.1 are about 2.1 A. The computed lowest-energy conformation of 2t, conformation 2t.1, indeed corresponds to the conformation of 2 found in the crystal [2]. The other six local minima have relative energies ranging from 4.6
to 8.8 kJ tool -l. The total number of conformations with relative energies less than 12.5 kJ mol -I is 10. Conformation 2t.2 has the same type of intramolecular hydrogen bonding as 2t.1. The preference of conformation 2t.1 is due to more advantageous angle bending and van der Waals interactions. Next in energy is a conformation with the C(3) hydroxyl group hydrogen bonded to the oxygen atom of the 16-ring (2t.3: AE = 5.1 kJ mol-l). Conformation 2t.4 is stabilized by the formation of a six-membered hydrogen bonded ring with the participation of the C(3) and C(5) hydroxyl groups of 2 AE = 5.9kJmol-l). The lowest-energy local minimum without hydrogen bond, 2t.6, is 8.5 kJ mo1-1 above the global minimum conformation. Hydrogen bonding through the ring (2t.10: O-H...O--C(9); 8-membered hydrogen bonded ring) is accompanied by a 12.5 kJ mol -l increase of the energy of the corresponding conformation. The solvation energy contribution does not affect markedly the computed ordering of conformations. Twenty-four conformations with relative energies less than 12.5 kJ moV I were found for 2c. Conformations 2c.3 and 2c.4 have exactly the same geometry of the macrocyclic ring and the same type of hydrogen bonding as conformation 2c.1. The lowest-energy local minimum without hydrogen bond, 2c.5, is 6.3 kJ mol -~ above the global minimum conformation 2c.1. Conformation 2c.2 (Fig. 4) deserves a brief comment: this conformation displays a completely different geometry for the 16-membered ring, as well as a different pattern for intramolecular hydrogen
125
P.M. lvanov/Journal of Molecular Structure 440 (1998) 121-130
Table 2 Computed conformational data of the 12-trans (2t) and 12-cis (2c) isomers of protynolide (2) 1 - 2 - 3 - 4 - 5 - 6 - 7 - 8 - 9 - 1 0 - 1 | - 12-13-14-15-0-1 12-trans-2 A A A G' G' A A A A G' G' A G' A A G' G' A G G' A G' G' A ~ A A ~ ' G' A G G' A G' G' G ~ G' A ~ G' ~ ~ ' G' G A A G' G A A A G G ~ A A G' G' A 12-cis-2 G A A G' G' G A G' G A A G G A A G' G' G ~ A A G' G' ~ G G' A G' G' G ~ ' A A G' G' ~ G G' A G' G' G ~ G' A G' G' ~ ~ G' A G' G' ~ G'G'A G'G'~ ~ G' A G' G' A G' A A G' G' G A G' A G' G' ~ G' G' A G' G' G ~ G' A G' G' ~ G G' A A A G' G A A a G G ~ G' A G' G' G A G' A G' G' G A G' A G' G' A G' G' A G' G' A A G' G A A G A G' A G' G' G G G' A A G G
G' G' G' G' G' G' G' G' G' G'
AE a (kJ mol-I)
Conformation
G' G' G' G' G' G' G' G ~' G'
A A G ~ ~ A G ~ A A
A A A A A A A A A A
A A G ~ A A A A A A
A A A A A A A A A A
G ~ G ~ ~ ~ ~ ~ ~ ~
G' G' G' G' G' G' G' G' G' G'
G G G' G ~ G ~ ~ G G
A A A A A A A A A A
0.0 4.6 5.1 5.9 6.8 8.5 8.8 12.3 12.3
2t.1 2t.2 2t.3 2t.4 2t.5 2t.6 2t.7 2t.8 2t.9
12.5
2t.lO
G' G A G G' G G' G A G G' G A G A G A G G'G G' G ~' G G' G' G' G G' G G' G G' G G' ~ A G G' G ~ ' G' A G A G A G
A E A A E A E E E A A A E A A G' a a E A E E E E
A A A A A A A A A A A A A A A A A A A A A A A A
A A A A A A A A A A G' ~ A A A A A A A A A A A A
E E E E E E E E E E E E E E E E E E E E E E E E
G G ~ G G G ~ G G G ~' ~ ~ G ~ G ~ G G G ~ G G G
G' G A G G' G G' G A G G' G' A G A G A G G'G'A G' G' G' G' G' G' G' G' G' G' G' G G' G G' ~ ' A G G' G' G' G A G A G A G
A A A A A A A A A
O. 1 2.4 3.0 4.2 6.4 6.7
2c.1 2c.2 2c.3 2c.4 2c.5 2c.6 2c.7 2c.8 2c.9 2c.10 2c.ll 2c.12 2c.13 2c.14 2c.15 2c.16 2c.17 2c.18 2c.19 2c.20 2c.21 2c.22 2c.23 2c.24
m
A A A A A A A A A A A A A A
7.1 7.5 8.4 8.7 8.8 9.4 9.6 9.8 10.0 10.1 10.4 11.3 11.4 11.4 11.7 11.8 12.1 12.3
aThe energies are given relative to the energy of the global minimum conformation of the 2t isomer (conformation 2t.l). b o n d i n g . A n e i g h t - m e m b e r e d h y d r o g e n b o n d e d ring is f o r m e d in this c a s e w i t h the p a r t i c i p a t i o n o f C ( 1 ) = O and the h y d r o x y l g r o u p at C(5). T h e c o n f o r m a t i o n has the m o s t d i s a d v a n t a g e o u s a n g l e b e n d i n g and van der W a a l s e n e r g y c o n t r i b u t i o n s , and the m o s t f a v o u r a b l e torsional and e l e c t r o s t a t i c terms. The O - H . . . O = C d i s t a n c e in 2c.2 is o n l y 1.82 A. A n a n a l o g o u s l o w - e n e r g y c o n f o r m a t i o n w a s not f o u n d for the 2t i s o m e r d u r i n g the M C M M search. Such a conformation of 2t may, however, be found if o t h e r s e a r c h p r o c e d u r e s or o t h e r m o l e c u l a r m e c h a n i c s m o d e l s are used.
3.2. T y l o s i n (1)
T h e c o n f o r m a t i o n s d i s c u s s e d in this s e c t i o n r e f e r to the i s o m e r w i t h c i s c o n f i g u r a t i o n a r o u n d the C ( 1 2 ) = C ( 1 3 ) d o u b l e b o n d , i.e. i s o m e r l c . O n l y nine c o n f o r m a t i o n s with relative e n e r g i e s less than 12.0 kJ m o l -j w e r e f o u n d for l c with the A M B E R * f o r c e field (Table 1). The d i h e d r a l 2 - 3 - 4 - 5 has the o b s e r v e d e x t e n d e d a n t i c o n f o r m a t i o n o f 2 [2] for this f r a g m e n t , and this c o m p l i e s with the four n e i g h b o u r ing all-equatorial s u b s t i t u t i o n s o f a n i n e - m e m b e r e d h y d r o g e n b o n d e d ring f o r m e d b y the C(3) h y d r o x y l
126
P.M. Ivanov/Journal of Molecular Structure 440 (1998) 121-130
lc.1
lc.mm3
Fig. 3. Stereographic views [ 18] of the lowest-energy conformations of lc computed with the AMBER* (lc.1) and MM3 (le.mm3) force fields.
group and the aldehyde oxygen (C(3)-OH . . . O = C H C H z C ( 6 ) , d i s t a n c e 1.80 A in the g l o b a l m i n i m u m c o n f o r m a t i o n ) . All n i n e c o n f o r m a t i o n s have the same pattern of intramolecular hydrogen b o n d i n g . E i g h t o f t h e m h a v e the s a m e c o n f o r m a t i o n o f the 1 6 - m e m b e r e d ring a n d d i f f e r f r o m e a c h o t h e r b y r o t a t i o n s a b o u t side c h a i n b o n d s . A local m i n i m u m w i t h r e l a t i v e e n e r g y A E = 10.7 kJ tool -l
(conformation lc,7) has different dihedrals about the C ( 6 ) - C ( 7 ) a n d C ( 1 4 ) - C ( 1 5 ) - O bonds. Conf o r m a t i o n l c . 7 h a s the m o s t f a v o u r a b l e n o n b o n d e d (van der Waals and electrostatic) interactions, but disadvantageous torsional and solvation energy c o n t r i b u t i o n s . E x c e p t for the d i f f e r e n c e in the c o n f i g u r a t i o n a b o u t the C ( 1 2 ) = C ( 1 3 ) d o u b l e b o n d , n o o n e a m o n g the A M B E R * c o m p u t e d c o n f o r m a t i o n s
Table 3 Dihedral angles (degrees) and calculated vicinal coupling constants (Hz) for protons in the macrocyclic ring of low-energy conformations of le, 2t and 2c Vicinal pair
H 2a./H3
H 2b/H3 H3/H4
Ha/H5 Hs/H6 H6]H 7a
H 6/I-I7b H7a/H 8 HTb/H8 ~Ref. [ll.
lc.1
2t.l
2c.l
Jaxp
2c.2
~b
Jcatc
~b
Jc.tc
~
Jcalc
~b
J~l~
179 65 50 -173 87 -166 79 92 -153
11.5 3.4 2.8 10.8 0.7 11.7 1.2 1.0 10.1
169 -70 58 178 67 -171 73 65 180
10.7 1.1 1.7 10.6 1.0 12.1 1.6 2.4 12.4
166 -78 71 -168 72 -140 105 70 -174
10.5 1.0 0.7 10.8 0.7 7.7 1.9 1.9 12.3
-163 82 --175 -58 -166 -177 --68 -65 179
11.3 1.6 10.7 1.7 9.5 12.3 2. I 2.4 12.4
10.3 0.8 < 1 9.1 1.3 10.1 3.4 5.0 11.0
P.M. lvanov/Journal of Molecular Structure 440 (1998) 121-130
127
Table 4 Relative energies of some conformations of le, 2t and 2e estimated with molecular mechanics force fields and semiempirical molecular orbital methods Conformation
AE~ (kJ mol -I) AMBER*
MM3
AM 1
PM3
0.0 4.9 15.2
14.6 10.0 0.0
2t.1
0.0
0.0
2c.1 2c.2 2c.3 2c.4
0.1 2.4 3.0 4.2
8.4 15.1 7.1 c
8.4 6.7 0.0 5.0 1.3 11.7 0.0 2.5
9.2 6.7 0.0 0.0 4.2 17.6 2.5 6.3
lc.l 1c.2 le.mm3 h
~Total energies of global minimum conformations (kJ mol-I): lc ( - 186.3 (AMBER*), 375.3 (MM3), - 1216.7918 x 105 (AM1), -1138.2315 x 1 0 3 (PM3)); 2t ( - 34.3 (AMBER*), 165.7 (MM3), - 461.4403 x 103 (PM3)); 2e ( - 489.5350 × 103 (AM1)). ~I'his conformation is not given in Table I. ~Converged to conformation 2c.3. of l c with relative energies less than 20.0 kJ mo1-1 has the same geometry o f the macrocycle as the observed ring conformation of 2t in the crystal [2].
3.3. Comparison of the computed data with geometries proposed from measured coupling constants and NOE enhancements [1] Table 3 presents dihedral angles and computed vicinal coupling constants [19] for some pairs of protons from the macrocyclic ring of the global m i n i m u m conformations l e A , 2t.1 and 2 e . l , as well as conformation 2e.2. In the rightmost column of the table are the measured couplings. Conformations lc.1 and 2t.1 have calculated vicinal coupling constants in accordance with the observed ones [ 1]. The macrocyclic ring conformation of the computed lowest-energy form of I t is very similar to structure lc.1 (a letter-code designation G ' G ' A G ' G ' A G ' G ' A A A A G G ' G ' A ) . Thus, the measured coupling constants cannot be used unequivocally to make an analogy for the preferred geometries of 1 and 2. The computed distances for pairs of protons in 2t.1 most closely match the NOE data [1]: H(3)/H(6) (2.2~,), and H(10)/H(21) (2.4 A). The global minimum conformation of the cis isomer, 2e.1, has the distances between protons H(3)/H(6) (2.1 A), H(3)/ H ( l l ) (3.3 A), while conformation 2c.2 is the only one with a short contact o f 2.2 ,~ for the pair H(8)/ H(11). The A M B E R * global minimum conformation
lc.1 has only one pair of protons in close proximity, H(3)/H(6) (2.1 A).
3.4. Force field dependence of the results. Semiempirical molecular orbital calculations of some low-energy conformations of lc, 2t and 2c All A M B E R * local minima of l c with relative energies less than 16.0 kJ mo1-1 were subsequently minimized also with the MM3* [16] force field parameters. The results confirmed the expectation for possible significant dependence of the relative energies of the conformations on the molecular mechanics force field used. A local minimum with A M B E R * relative energy 1 5 . 2 k J m o l -~ now has the lowest energy with MM3*. This conformation (designated l c . m m 3 in Table 4) and the conformations l e a and lc.2 were further calculated with the original MM3 version of Allinger's force field [8]. The same calculations were carried out also for 2t.1 and the four lowest energy conformations of 2c obtained with AMBER*. These results are compared in Table 4 with the relative energies of the same structures estimated with the semiempirical molecular orbital parametrizations AM1 [20] and PM3 [21], using M O P A C 93.0 [22]. MM3 and the semiempirical molecular orbital methods are consistent in predicting the conformation l c . m m 3 to have the lowest energy (Table 4). The three methods reverse the order of preference for
128
P.M. lvanov/Journal of Molecular Structure 440 (1998) 121-130
2t.1
2c.1
2c.2 Fig. 4. Stereographic views 118] of low-energy conformations of 2t and 2e.
P.M. lvanov/Journal of Molecular Structure 440 (1998) 121-130
conformations lc.1 and lc.mm3, computed with the AMBER* force field. The preference of lc.1 with AMBER* is determined by more favourable torsional and electrostatic contributions. A completely different balance of interactions is displayed by the MM3 force field. The angle bending and torsional terms have the decisive role in lowering the energy of conformation lc.mm3 when the MM3 force field is used. The computed vicinal coupling constants for pairs of protons in the macrocyclic ring of lc.mm3 are also in accordance with the corresponding measured couplings [1], and this conformation has the same short distances for pairs of protons as structure 2c.1. In accordance with the AMBER* results, MM3 and PM3 predicted the conformation 2t.1 to have lower energy than the four conformations of 2c considered. Only AM 1 gives a significantly different order of stability for the set of structures of 2t and 2c given in Table 4. In summary, the analysis of the conformational data obtained from the search for low-energy conformations of lc, 2t and 2c with the AMBER* force field, and the comparative studies of some local minima also with MM3, AM 1 and PM3, suggest that macrocyclic ring conformations having the hydrogen bonding patterns of lc.1 and lc.mm3 also have to be considered among the probable conformations of tylosin in solution.
4. Conclusions A Monte Carlo conformational search technique was used with the AMBER* force field as a model to search for low-energy conformations of tylosin, 1, and protynolide, 2. Selected optimized structures obtained from the search were also further studied with the MM3 force field [8] and semiempirical molecular orbital methods [20-22]. Structural and energy data were obtained for the 12-cis isomer of 1 (lc), and the 12-trans (2t) and 12-cis (2c) isomers of 2. The diene system at C(10), C(ll), C(12) and C(13) of 2 has the t r a n s - t r a n s configuration [2]. The NOE experiments show that the two C=C double bonds of 1 have the s - t r a n s arrangement [1 ]. The computed data for lc, 2t and 2c provide grounds for conclusions to be drawn about probable conformations of 1 and 2 in solution. Irrespective of the configuration around the C(12)=C(13) double bond, the 12-atoms fragment
129
of the 16-membered ring, not containing the diene system, has enough flexibility to adopt different low-energy hydrogen bonded conformations. The answers to the questions formulated at the beginning can be briefly summarized as follows. The geometry of the 16-membered ring of the AMBER* computed low-energy conformations of le, as well as the lowest-energy conformation of lt, is not the same as the geometry of the macrocycle in the crystal of 2 [2], nor does it resemble the geometry of the low-energy conformations of 2t and 2c that were computed with the same force field. Intramolecular hydrogen bonding is among the important factors which determine the preferred geometry of the macrocycle. The pattern of intramolecular hydrogen bonding is also different in the most probable conformations of 1 and 2 found with the AMBER* force field. The lowest energy conformation of 2t, conformation 2t.1, indeed corresponds to the conformation found for 2 in the crystal [2]. The computed relative energies of the conformations of the isomers of 1 and 2 strongly depend on the molecular mechanics force field used. Thus, in addition to the most probable conformation lc.1 of lc found from the Monte Carlo search with the AMBER* force field, we propose a conformation with the 16-ring conformation and hydrogen bonding pattern of lc.mm3 (Fig. 3), also to be considered as a possible candidate for a low-energy conformation of 1. Other conformations may also be found if other search procedures are applied [5], or if more advanced molecular mechanics models are used (see [23] for some recent developments). Testing different values for the arguments of the MCMM (Monte Carlo Multiple Minimum), TORS (variable TORSion selection) and RCA4 (Ring Closure Atoms(4)) commands [15,16], and/or adding also the TORC (TORsional Constraint) command [15,16] of the search procedure used in our study, could also influence the results. Such a systematic approach, requiring significant computational resources, is, however, beyond our reach at present.
Acknowledgements Financial support from the Direcci6 General de Universitats from Generalitat de Catalunya, through
130
P.M. Ivanov/Journal of Molecular Structure 440 (1998) 121-130
a grant for a visiting professorship at the Universitat Aut6noma de Barcelona, Spain (1995), is gratefully acknowledged. Thanks are due also to Prof. Carlos Jaime for providing computational resources, and to Prof. E. Osawa for facilities for preparing the manuscript.
References
[10] [11] [12] [13]
[1] S. Simova, G. Ivanova, Magn. Reson. Chem. 34 (1996) 255. [2] S. Omura, H. Matsubara, A. Nakagawa, A. Fukusaki, T. Matsumoto, J. Antibiot. 33 (1980) 915. [3] P.M. Ivanov, J. Mol. Struct., in press. [4] (a) P.M. Ivanov, C. Jaime, Anal. Quim. Int. Ed. 92 (1996) 13. (b) P.M. Ivanov, C. Jaime, J. Mol. Struct. 377 (1996) 137. [5] (a) H. Goto, E. Osawa, J. Mol. Struct. (Theochem) 285 (1993) 157. (b) H. Goto, E. Osawa, J. Chem. Soc., Perkin Trans. 2, (1993) 187. (c) F. Villamagna, M.A. Whitehead, J. Chem. Soc., Faraday Trans. 96 (1994) 47. (d) J. Koca, S. Perez, A. Imberty, J. Comput. Chem. 16 (1995) 296. (e) L. DosenMicovic, Tetrahedron 51 (1995) 6789. (f) S.B. Engelsen, J. Koca, I. Braccini, C.H. du Penhoat, S. Perez, Carbohydr. Res. 276 (1995) 1. (g) A. Smellie, S.D. Kahn, S.L. Teig, 25 (1995) 285. (h) J.C. Meza, R.S. Judson, T.R. Faulkner, A.M. Treasurywala, J. Comput. Chem. 17 (1996) 1142. (i) I. Kolossvary, W.C. Guida, J. Am. Chem. Soc. 118 (1996) 5011. (j) C.S. Wang, J. Comput. Chem. 18 (1997) 277. [6] (a) N.L. Allinger, J. Am. Chem. Soc. 99 (1977) 8127. (b) U. Burkert, N.L. Allinger, Molecular Mechanics, American Chemical Society, Washington, DC, 1982. [7] C. Jaime, E. Osawa, Tetrahedron 39 (1983) 2769. [8] (a) N.L. Allinger, Y.H. Yuh, J.H. Lii, J. Am. Chem. Soc. 111 (1989) 8551; 8566; 8576. (b) MM3(92): Technical Utility Corporation, 235 Glen Village Court, Powell, OH 43065, USA. [9] (a) P.M. Ivanov, J.M. Rudzinski, E. Osawa, results presented at the 13th Symposium on Chemical Information and Computer Sciences, 28-29 November 1990, Toyohashi,
[14]
[15]
[16]
[ 17] [18]
[19] [20] [21] [22]
Japan. (b) E. Osawa, P.M. Ivanov, J.M. Rudzinski, KESSHOU: Molecular mechanics program for molecular and crystal structure optimization, presented in part in: E. Osawa, Potential energy calculations of crystals, Ch. 1. I of Reactivity in Molecular Crystals, Y. Ohashi (Ed.), Kodansha, Tokyo and VCH, Weinheim, 1993, pp. 1-9. (c) P.M. Ivanov, J. Mol. Struct., in press. P.M. Ivanov, unpublished results, 1991. N. Norman, H. Mathisen, Acta Chem. Scand. 15 (1961) 1755. P.M. Ivanov, T.G. Momchilova, I.G. Pojarlieff, J. Comput. Chem. 12 (1991) 427. C.M. Breneman, K.B. Wiberg, J. Comput. Chem. 11 (1990) 361. (a) P.M. Ivanov, Jpn Chem. Program Exch. Newslett. 7(3) (1995) 3. (b) P.M. Ivanov, Jpn Chem. Program Exch. Newslett. 7(4) (1996) 13. (a) G. Chang, W.C. Guida, W.C. Still, J. Am. Chem. Soc. 111 (1989) 4379. (b) M. Saunders, K.N. Houk, Y.D. Wu, W.C. Still, M. Lipton, G. Chang, W.C. Guida, J. Am. Chem. Soc. 112 (1990) 1419. F. Mohamadi, N.G.J. Richards, W.C. Guida, R. Liskamp, C. Caufield, G. Chang, T. Hendrickson, W.C. Still, J. Comput. Chem. 11 (1990) 440. W.C. Still, A. Tempczyk, R.C. Hawley, T. Hendrickson, J. Am. Chem. Soc. 112 (1990) 6127. M.N. Burnett, C.K. Johnson, ORTEP-III: Oak Ridge thermal ellipsoid plot program for crystal structure illustrations, Oak Ridge National Laboratory Report ORNL-6895, 1996. C.A.G. Haasnoot, F.A.A.M. De Leeuw, C. Altona, Tetrahedron 36 (1980) 2783. M.J.S. Dewar, E.G. Zoebisch, E.F. Healy, J.J.P. Stewart, J. Am. Chem. Soc. 107 (1985) 3902. J.J.P. Stewart, J. Comput. Chem. 10 (1989) 221. MOPAC 93.00: J.J.P. Stewart, Fujitsu Ltd, Tokyo, Japan, 1993.
[23] (a) AMBER: W.D. Cornell, P. Cieplak, Ch.I. Bayly, I.R. Gould, K.M. Merz, Jr, D.M. Ferguson, D.C. Spellmeyer, Th. Fox, J.W. Caldwell, P.A. Kollman, J. Am. Chem. Soc. 117 (1995) 5179. (b) MM4: N.L. Allinger, K. Chen, J.H. Lii, J. Comput. Chem. 17 (1996) 642. (c) Merck Molecular Force Field: T.A. Halgren, J. Comput. Chem. 17 (1996) 616.
MOLECULAR STRUCTURE ELSEVIER
Journal of Molecular Structure 440 (1998) 121- 130
Molecular mechanics conformational analysis of tylosin P e t k o M. I v a n o v ~a,b ~lnstitute of Organic Chemistry. with Centre of Phytochemistry, Bulgarian Academy of Sciences, ul. Acad. G. Bonchev, bloc 9, 1040 Sofia, Bulgaria bDepartament de Quimica, Facultat de Cikncies, Universitat Autbnoma de Barcelona, 08193 Bellaterra, Barcelona, Spain
Received 17 March 1997; revised 16 May 1997; accepted 16 May 1997
Abstract
The conformations of the 16-membered macrolide antibiotic tylosin were studied with molecular mechanics (AMBER* force field) including modelling of the effect of the solvent on the conformational preferences (GB/SA). A Monte Carlo conformational search procedure was used for finding the most probable low-energy conformations. The present study provides complementary data to recently reported analysis of the conformations of tylosin based on NMR techniques. A search for the low-energy conformations of protynolide, a 16-membered lactone containing the same aglycone as tylosin, was also carried out, and the results were compared with the observed conformation in the crystal as well as with the most probable conformations of the macrocyclic ring of tylosin. The dependence of the results on force field was also studied by utilizing the MM3 force field. Some particular conformations were computed with the semiempirical molecular orbital methods AMI and PM3. © 1998 Elsevier Science B.V. Keywords: tylosin; protynolide; conformation; molecular mechanics; AMBER*; MM3; MC search
1. I n t r o d u c t i o n
An assignment of the IH and 13C NMR spectra of tylosin (1, Fig. 1) in CDCI 3 was recently reported [1], and conclusions had been also drawn regarding the solution conformation of the macrocyclic ring of 1 in aprotic solvents, based on vicinal proton-proton coupling constants and NOE enhancements. The authors further supported their experimental data with results from molecular mechanics calculations of protynolide (2), a 16-membered lactone containing the same aglycone as tylosin [1]. Two low-energy conformations of 2 were found in that study and i E-mail: [email protected] 2 Recent references to papers related to techniquesfor the search of low-energyconformationsare listed in [5].
discussed as possible candidates for the preferred conformations of 1 in solution. The predominance of the conformation corresponding to the conformation of 2 in the crystal [2] is proposed (conformation 2t.1 in Fig. 4). Not always, however, does the congruence between averaged computed quantities and measured observables provide a guarantee for assessing the conformational ordering, and different computational methods can be at variance in the details of the computed conformational distribution [3]. The results from the computational conformational analysis of a macrocyclic system of the size of tylosin could happen to depend significantly on the molecular mechanics model used [4], as well as on the way the search for the most probable conformations is carried out. 2 The author recently stressed the importance of
0022-2860/98/$19.00 © 1998 Elsevier Science B.V. All rights reserved PII S0022-2860(97)00233-0
122
P.M. lvanov/Journal of Molecular Structure 440 (1998) 121-130
IIt° 6
"1
lvle ~ . . . . . ~ ' 1 1
12 1112
"' Me ~5'~...---O
HO
23
-4 /]J 13
t.,,
,9
7 ~,.
18
661
~le
~
0
\
20
t CH2CHO ........
,
1'
a
0
5'
/X---o
-
6'
...
4'
..
a o ~ ' ~ OMe
16 [ / ~ ~ " . . . . . . , ,
OH
171
/ N Me 7'
2,, ~ Me 8'
~ u
~T~ 7" ~ ) H 5"\ ~1" Me 6 t*
Fig. 1. Numbering of atoms in tylosin.
s e l e c t i n g the p r o p e r m o l e c u l a r m e c h a n i c s m o d e l w h e n i n c l u s i o n c o m p l e x e s are studied [4]. T h e s a m e r e m a r k s are v a l i d also for the case o f a large m o l e c u l e w h e r e d i s t a n t f r a g m e n t s o f o n e a n d the s a m e m o l e c u l e h a v e the f r e e d o m to i n t e r a c t in a w a y as if t h e y belonged to different molecules. 3 Molecular m e c h a n i c s m o d e l s h a v e b e e n d e v e l o p e d m o s t l y to m a t c h specific p u r p o s e s ; 4 thus, it s h o u l d n o t b e s u r p r i s i n g if t h e y d i s p l a y a d i v e r s e b a l a n c e o f intera c t i o n s in s o m e cases. T h e m o l e c u l a r m e c h a n i c s conformational analysis of tylosin presents such an example. T h e p r e s e n t c o m p u t a t i o n a l s t u d y is a i m e d at finding low-energy conformations of 1 and 2 u s i n g the o n l y a v a i l a b l e (to us) s t a n d a r d s e a r c h t e c h n i q u e [15] i m p l e m e n t e d in the M a c r o M o d e l and BatchMin V4.0 molecular modelling p r o g r a m s [16]. T h e f o l l o w i n g q u e s t i o n s c a n be formulated that deserve elucidation in this c o m p u t a t i o n a l study: 1. D o 1 a n d 2 h a v e the s a m e o r d e r i n g o f c o n f o r m a t i o n s d i s p l a y e d b y the g e o m e t r i e s o f the 16-membered ring? 2. A r e t h e r e o n l y o n e or t w o l o w - e n e r g y c o n f o r m a t i o n s o f 2, as r e p o r t e d b y o t h e r a u t h o r s [1], a n d w h a t are the f o r c e s w h i c h r e s t r i c t the p o t e n t i a l for flexibility o f the 1 6 - m e m b e r e d ring to o n l y o n e or t w o s t r o n g l y f a v o u r e d forms? 3. Is the crystal c o n f o r m a t i o n o f 2 a m o n g the m o s t s t a b l e f o r m s o f the 1 6 - m e m b e r e d r i n g s o f 1 a n d 2 in s o l u t i o n ?
3Crystal calculations using molecular mechanics are an appropriate test for the ability of a force field to handle intermolecular and distant intramolecular van der Waals interactions. Even tiny modifications of the van der Waals parameters can affect significantly the balance of nonbonded interactions in large molecules or crystals. A comparison between three generic parametrizations, MM2 [6], MM2' [7] and MM3 [8], carried out with the KESSHOU program [9] (see below), can serve to illustrate the problem [10]. Calculations of the n-hexane crystal using the MM3 van der Waals potentials for atom-atom interactions closely reproduce observed unit cell constants. This is not unexpected, because crystal structural data had been extensively used for the parametrization of the MM3 force field [8]. MM2 produces about 15.0 ~3 smaller unit cell volume than the one found by experiment [11]. The MM2' parametrization [7], which was specifically developed with the intention of improving the performance of MM2 in reproducing rotational barriers, gives even larger deviations from experimental crystal data; i.e., the balance of nonbonded interactions has been sacrificed in MM2' for the purpose of improving performance in reproducing rotational barriers. It is well known that MM2 enhances the stabilization of the conformations when van der Waals interactions form the prevailing contribution to the energy balance of intramolecular interactions [3,12]; thus, using MM2 for studying the conformations of protonolide 2 was not the best choice the authors could have made [1]. Besides, modelling of hydrogen bonding was also not included [1]. KESSHOU [9] was recently further developed to handle also electrostatic interactions [9c]. MM3 crystal calculations with full molecular geometry optimization, using CHELPG electrostatic potential derived atomic charges [13], gave perfect reproduction of the packing of the 1.4dicyanobenzene molecules in the crystal ttriclinic crystal lattice, space group PL Z = 1, crystal radius 20.0 A, 373 unit cells considered; full geometry optimization of the reference molecule). The packing of the C60 molecules in the low-temperature phase of crystalline C~ is also exactly reproduced (cubic, Pa3, Z = 4, crystal radius 50.0 A, 332 unit cells). These results certify the improved quality of the MM3 van der Waals potential functions and parameters for modelling nonbonded atom-atom interactions. 4The present status of molecular mechanics force field development was recently reviewed by the author [14].
123
P.M. lvanov/Journal of Molecular Structure 440 (1998) 121-130
2. Computational details The Monte Carlo conformational searcher [15,16] in BatchMin was used (MCMM command) with the AMBER* force field as a model, and the parameters for CHC13 were adopted for modelling the solvent with Still's GB/SA method [17]. Extended cutoff distances for the nonbonded (van der Waals and electrostatic) interactions were set to 20.0 ,~. For protynolide, 2, the maximum number of MC steps was set to 6000 using the global search mode as a parameter for the MCMM command. The minimum and maximum number of degrees of freedom changed in an MC step were set to 2 and 11, respectively. The MC starting structure selection command (MCSS) was executed with use-directed mode (argl = 2) and an energy window of 25.0kJ tool <. Optimized structures with relative energies less than 25.0 kJ mol -j were kept during the search (the DEMX command), while structures with relative energy (after 200 minimization iterations) higher than 50.0 kJ mol -~ were excluded from further energy minimization. All 28 non-hydrogen atoms were used for comparison when eliminating duplicate minima. Twelve endocyclic torsion angles were varied. The four torsional angles of the C = C - C = C - C O fragment were not included in the list of variable torsion angles. The default values for a minimum and maximum dihedral angle variation were adopted. Conjugate gradients, followed by Newton-Raphson minimization, were applied in the final calculations of the 198 unique structures stored after the MCMM search was completed. We followed the same computational scheme for 1. Twenty-one torsional angles were assumed to vary in the case of tylosin. The number of heavy atoms used for comparison when eliminating duplicate minima is 52. These comprise the atoms of the 16-membered ring and the most important atoms from the side chains. The list of 21 variable torsional angles includes the same dihedrals as in the case of the protynolide (2), plus all dihedrals defining the rotational isomerism about side chains bonds. Endocyclic dihedrals of the six-membered rings were not included in the list of variable torsional angles; thus ring closure conditions were imposed only for the 16-membered ring. The starting conformation had the observed conformation of the
macrocycle of 2 in the crystal. The maximum number of MC steps was set to 10 000 using the global search mode as a parameter for the MCMM command. The search ended with 494 unique structures stored. Only three of them satisfied the criteria for good convergence. Thus, all stored structures were further subjected to consecutive conjugate gradients and Newton-Raphson optimizations. Restrictions were not imposed to keep the starting trans configurations around the C ( 1 0 ) = C ( l l ) and C(12)=C(13) double bonds during the MCMM search. Optimized low-energy structures were obtained for the 12-trans and 12-cis isomers of 2 and are designated as 2t and 2c, respectively. All low-energy structures obtained for 1 have cis configuration around the C(12)=C(13) bond. The global minimum conformation of 12-trans-1 is with 19.0 kJ mol -~ higher energy than the global minimum conformation of 12-cis-1. Only the conformations of the 12-cis isomer of 1 will be discussed. These conformations are designated lc.
3. Results and discussion Computed conformational data for lc are summarized in Table 1. Structural and energy data of 2t and 2c are presented in Table 2. These are relative energies of the low-energy conformations obtained with the AMBER* force field, as well as a letter-code representation of important dihedral angles. The numbering of the atoms of tylosin is displayed in Fig. 1. For the sake of compactness of the presentation we used a letter-code designations of dihedral angles, as given in Fig. 2. A letter G designates a dihedral angle in the range from 20 ° to 100 °, G from 100° to 150 °, A from 150° to - 150 °, etc. .20°
-1~0
20°
--
~50°
Fig. 2. Designation of dihedral angles.
124
P.M. lvanov/Journal of Molecular Structure 440 (1998) 121-130
Table 1 Computed conformational data of the 12-cis isomer (le) of tylosin (1)" 1-2-3-4-5-6-7-8-9-
I0-11-12-13-14-15-0-1
G' G' A G ' G' A G ' G A A ~ ' G' A G' G' A G' G A A ~ ' G' A G' G' A G' G A A ~ ' G' A G' G' A G' G A A ~ ' G' A G' G' A G' G A A ~ ' G' g G' G' g G' G g g ~G'AG'G'GG'AAAG'E~A ~ ' G' A G' G' A G' G A A ~ ' G' A G' G' A G' G A A
I
a
n
A A A A A g
E G G' EGG' EGG' EGG' EGG' E~G'
G' G' G' G' G' G' G A E ~ G ' G' A E ~ G ' G' rr
A A A A A A A A A
Xi
X2
X3
X4
X5
G
G G G G G G G G G'
G' ~' ~' ~ G' ~' G' ~' ~'
A A A G A G g A A
G G' ~ A G G ~' ~ G
~' G ~' G' m
AE (kJ mol )) 0.0 4.9 7.7 7.8 8.9 10.0 10.7 I1.9 12.0
Conformation
le.1 le.2 lc.3 lc.4 lc.5 lc.6 le.7 lc.8 to.9
.
X~=3-4-O-1 X2=3-4-O-H X3=3 - 2 - O - M e X 4 = 1 5 - 1 4 - 2 3 - O Xs = 1 4 - 2 3 - O - l " . Other side chain dihedral angles have the same values in all conformations: 2 - 3 - O - H (G'), 4 - 5 - O - 1' (G'), 5 - 6 - 1 9 - 2 0 (G'), 1 4 - 1 5 - 1 6 - 1 7 (A), 6 - 1 9 - C ( 2 0 ) = O (E), 5 - O - 1' - 2 ' (A), l ' - 2 ' - O - H (G), 2 ' - 3 ' - N - M e ( 7 ' ) (G), 4"-3"Y/-O-H (G'), 2 3 - O - 1 ' " - O (G'L 4 " - 3 " - O - M e (G'), 3 " - 4 " - O - H (G').
Table 3 contains dihedral angles and calculated vicinal coupling constants for some pairs of protons in the macrocyclic ring of low-energy conformations of lc, 2t and 2c. Relative energies of the most probable conformations estimated with different molecular mechanics force fields and semiempirical molecular orbital methods are given in Table 4. Stereographic views [18] of some conformations of lc are displayed in Fig. 3. Analogous presentations of conformations of 2t and 2e are given in Fig. 4. 3.1. Protynolide (2)
The computed lowest-energy structures of 2t and 2c are practically of the same energy. Except for the configuration around the C(12)~C(13) double bond, these two structures differ mostly in the dihedral angle C(7)-C(8)-C(9)-C(10), -59.0 ° and 104.0 ° in 2t.1 and 2c.1, respectively. Six-membered hydrogen bonded rings are formed in 2t.l and 2c.1 and these include the C(3) hydroxyl group and the C(1)=O oxygen. The more favourable torsional energy contribution (AEtors = 7,1 kJ mol -I) for the trans isomer, 2t.1, is balanced by advantageous van der Waals and electrostatic interactions in the cis form, 2c.1. The hydrogen bond distances O - H . . . O = C in 2t.1 and 2e.1 are about 2.1 A. The computed lowest-energy conformation of 2t, conformation 2t.1, indeed corresponds to the conformation of 2 found in the crystal [2]. The other six local minima have relative energies ranging from 4.6
to 8.8 kJ tool -l. The total number of conformations with relative energies less than 12.5 kJ mol -I is 10. Conformation 2t.2 has the same type of intramolecular hydrogen bonding as 2t.1. The preference of conformation 2t.1 is due to more advantageous angle bending and van der Waals interactions. Next in energy is a conformation with the C(3) hydroxyl group hydrogen bonded to the oxygen atom of the 16-ring (2t.3: AE = 5.1 kJ mol-l). Conformation 2t.4 is stabilized by the formation of a six-membered hydrogen bonded ring with the participation of the C(3) and C(5) hydroxyl groups of 2 AE = 5.9kJmol-l). The lowest-energy local minimum without hydrogen bond, 2t.6, is 8.5 kJ mo1-1 above the global minimum conformation. Hydrogen bonding through the ring (2t.10: O-H...O--C(9); 8-membered hydrogen bonded ring) is accompanied by a 12.5 kJ mol -l increase of the energy of the corresponding conformation. The solvation energy contribution does not affect markedly the computed ordering of conformations. Twenty-four conformations with relative energies less than 12.5 kJ moV I were found for 2c. Conformations 2c.3 and 2c.4 have exactly the same geometry of the macrocyclic ring and the same type of hydrogen bonding as conformation 2c.1. The lowest-energy local minimum without hydrogen bond, 2c.5, is 6.3 kJ mol -~ above the global minimum conformation 2c.1. Conformation 2c.2 (Fig. 4) deserves a brief comment: this conformation displays a completely different geometry for the 16-membered ring, as well as a different pattern for intramolecular hydrogen
125
P.M. lvanov/Journal of Molecular Structure 440 (1998) 121-130
Table 2 Computed conformational data of the 12-trans (2t) and 12-cis (2c) isomers of protynolide (2) 1 - 2 - 3 - 4 - 5 - 6 - 7 - 8 - 9 - 1 0 - 1 | - 12-13-14-15-0-1 12-trans-2 A A A G' G' A A A A G' G' A G' A A G' G' A G G' A G' G' A ~ A A ~ ' G' A G G' A G' G' G ~ G' A ~ G' ~ ~ ' G' G A A G' G A A A G G ~ A A G' G' A 12-cis-2 G A A G' G' G A G' G A A G G A A G' G' G ~ A A G' G' ~ G G' A G' G' G ~ ' A A G' G' ~ G G' A G' G' G ~ G' A G' G' ~ ~ G' A G' G' ~ G'G'A G'G'~ ~ G' A G' G' A G' A A G' G' G A G' A G' G' ~ G' G' A G' G' G ~ G' A G' G' ~ G G' A A A G' G A A a G G ~ G' A G' G' G A G' A G' G' G A G' A G' G' A G' G' A G' G' A A G' G A A G A G' A G' G' G G G' A A G G
G' G' G' G' G' G' G' G' G' G'
AE a (kJ mol-I)
Conformation
G' G' G' G' G' G' G' G ~' G'
A A G ~ ~ A G ~ A A
A A A A A A A A A A
A A G ~ A A A A A A
A A A A A A A A A A
G ~ G ~ ~ ~ ~ ~ ~ ~
G' G' G' G' G' G' G' G' G' G'
G G G' G ~ G ~ ~ G G
A A A A A A A A A A
0.0 4.6 5.1 5.9 6.8 8.5 8.8 12.3 12.3
2t.1 2t.2 2t.3 2t.4 2t.5 2t.6 2t.7 2t.8 2t.9
12.5
2t.lO
G' G A G G' G G' G A G G' G A G A G A G G'G G' G ~' G G' G' G' G G' G G' G G' G G' ~ A G G' G ~ ' G' A G A G A G
A E A A E A E E E A A A E A A G' a a E A E E E E
A A A A A A A A A A A A A A A A A A A A A A A A
A A A A A A A A A A G' ~ A A A A A A A A A A A A
E E E E E E E E E E E E E E E E E E E E E E E E
G G ~ G G G ~ G G G ~' ~ ~ G ~ G ~ G G G ~ G G G
G' G A G G' G G' G A G G' G' A G A G A G G'G'A G' G' G' G' G' G' G' G' G' G' G' G G' G G' ~ ' A G G' G' G' G A G A G A G
A A A A A A A A A
O. 1 2.4 3.0 4.2 6.4 6.7
2c.1 2c.2 2c.3 2c.4 2c.5 2c.6 2c.7 2c.8 2c.9 2c.10 2c.ll 2c.12 2c.13 2c.14 2c.15 2c.16 2c.17 2c.18 2c.19 2c.20 2c.21 2c.22 2c.23 2c.24
m
A A A A A A A A A A A A A A
7.1 7.5 8.4 8.7 8.8 9.4 9.6 9.8 10.0 10.1 10.4 11.3 11.4 11.4 11.7 11.8 12.1 12.3
aThe energies are given relative to the energy of the global minimum conformation of the 2t isomer (conformation 2t.l). b o n d i n g . A n e i g h t - m e m b e r e d h y d r o g e n b o n d e d ring is f o r m e d in this c a s e w i t h the p a r t i c i p a t i o n o f C ( 1 ) = O and the h y d r o x y l g r o u p at C(5). T h e c o n f o r m a t i o n has the m o s t d i s a d v a n t a g e o u s a n g l e b e n d i n g and van der W a a l s e n e r g y c o n t r i b u t i o n s , and the m o s t f a v o u r a b l e torsional and e l e c t r o s t a t i c terms. The O - H . . . O = C d i s t a n c e in 2c.2 is o n l y 1.82 A. A n a n a l o g o u s l o w - e n e r g y c o n f o r m a t i o n w a s not f o u n d for the 2t i s o m e r d u r i n g the M C M M search. Such a conformation of 2t may, however, be found if o t h e r s e a r c h p r o c e d u r e s or o t h e r m o l e c u l a r m e c h a n i c s m o d e l s are used.
3.2. T y l o s i n (1)
T h e c o n f o r m a t i o n s d i s c u s s e d in this s e c t i o n r e f e r to the i s o m e r w i t h c i s c o n f i g u r a t i o n a r o u n d the C ( 1 2 ) = C ( 1 3 ) d o u b l e b o n d , i.e. i s o m e r l c . O n l y nine c o n f o r m a t i o n s with relative e n e r g i e s less than 12.0 kJ m o l -j w e r e f o u n d for l c with the A M B E R * f o r c e field (Table 1). The d i h e d r a l 2 - 3 - 4 - 5 has the o b s e r v e d e x t e n d e d a n t i c o n f o r m a t i o n o f 2 [2] for this f r a g m e n t , and this c o m p l i e s with the four n e i g h b o u r ing all-equatorial s u b s t i t u t i o n s o f a n i n e - m e m b e r e d h y d r o g e n b o n d e d ring f o r m e d b y the C(3) h y d r o x y l
126
P.M. Ivanov/Journal of Molecular Structure 440 (1998) 121-130
lc.1
lc.mm3
Fig. 3. Stereographic views [ 18] of the lowest-energy conformations of lc computed with the AMBER* (lc.1) and MM3 (le.mm3) force fields.
group and the aldehyde oxygen (C(3)-OH . . . O = C H C H z C ( 6 ) , d i s t a n c e 1.80 A in the g l o b a l m i n i m u m c o n f o r m a t i o n ) . All n i n e c o n f o r m a t i o n s have the same pattern of intramolecular hydrogen b o n d i n g . E i g h t o f t h e m h a v e the s a m e c o n f o r m a t i o n o f the 1 6 - m e m b e r e d ring a n d d i f f e r f r o m e a c h o t h e r b y r o t a t i o n s a b o u t side c h a i n b o n d s . A local m i n i m u m w i t h r e l a t i v e e n e r g y A E = 10.7 kJ tool -l
(conformation lc,7) has different dihedrals about the C ( 6 ) - C ( 7 ) a n d C ( 1 4 ) - C ( 1 5 ) - O bonds. Conf o r m a t i o n l c . 7 h a s the m o s t f a v o u r a b l e n o n b o n d e d (van der Waals and electrostatic) interactions, but disadvantageous torsional and solvation energy c o n t r i b u t i o n s . E x c e p t for the d i f f e r e n c e in the c o n f i g u r a t i o n a b o u t the C ( 1 2 ) = C ( 1 3 ) d o u b l e b o n d , n o o n e a m o n g the A M B E R * c o m p u t e d c o n f o r m a t i o n s
Table 3 Dihedral angles (degrees) and calculated vicinal coupling constants (Hz) for protons in the macrocyclic ring of low-energy conformations of le, 2t and 2c Vicinal pair
H 2a./H3
H 2b/H3 H3/H4
Ha/H5 Hs/H6 H6]H 7a
H 6/I-I7b H7a/H 8 HTb/H8 ~Ref. [ll.
lc.1
2t.l
2c.l
Jaxp
2c.2
~b
Jcatc
~b
Jc.tc
~
Jcalc
~b
J~l~
179 65 50 -173 87 -166 79 92 -153
11.5 3.4 2.8 10.8 0.7 11.7 1.2 1.0 10.1
169 -70 58 178 67 -171 73 65 180
10.7 1.1 1.7 10.6 1.0 12.1 1.6 2.4 12.4
166 -78 71 -168 72 -140 105 70 -174
10.5 1.0 0.7 10.8 0.7 7.7 1.9 1.9 12.3
-163 82 --175 -58 -166 -177 --68 -65 179
11.3 1.6 10.7 1.7 9.5 12.3 2. I 2.4 12.4
10.3 0.8 < 1 9.1 1.3 10.1 3.4 5.0 11.0
P.M. lvanov/Journal of Molecular Structure 440 (1998) 121-130
127
Table 4 Relative energies of some conformations of le, 2t and 2e estimated with molecular mechanics force fields and semiempirical molecular orbital methods Conformation
AE~ (kJ mol -I) AMBER*
MM3
AM 1
PM3
0.0 4.9 15.2
14.6 10.0 0.0
2t.1
0.0
0.0
2c.1 2c.2 2c.3 2c.4
0.1 2.4 3.0 4.2
8.4 15.1 7.1 c
8.4 6.7 0.0 5.0 1.3 11.7 0.0 2.5
9.2 6.7 0.0 0.0 4.2 17.6 2.5 6.3
lc.l 1c.2 le.mm3 h
~Total energies of global minimum conformations (kJ mol-I): lc ( - 186.3 (AMBER*), 375.3 (MM3), - 1216.7918 x 105 (AM1), -1138.2315 x 1 0 3 (PM3)); 2t ( - 34.3 (AMBER*), 165.7 (MM3), - 461.4403 x 103 (PM3)); 2e ( - 489.5350 × 103 (AM1)). ~I'his conformation is not given in Table I. ~Converged to conformation 2c.3. of l c with relative energies less than 20.0 kJ mo1-1 has the same geometry o f the macrocycle as the observed ring conformation of 2t in the crystal [2].
3.3. Comparison of the computed data with geometries proposed from measured coupling constants and NOE enhancements [1] Table 3 presents dihedral angles and computed vicinal coupling constants [19] for some pairs of protons from the macrocyclic ring of the global m i n i m u m conformations l e A , 2t.1 and 2 e . l , as well as conformation 2e.2. In the rightmost column of the table are the measured couplings. Conformations lc.1 and 2t.1 have calculated vicinal coupling constants in accordance with the observed ones [ 1]. The macrocyclic ring conformation of the computed lowest-energy form of I t is very similar to structure lc.1 (a letter-code designation G ' G ' A G ' G ' A G ' G ' A A A A G G ' G ' A ) . Thus, the measured coupling constants cannot be used unequivocally to make an analogy for the preferred geometries of 1 and 2. The computed distances for pairs of protons in 2t.1 most closely match the NOE data [1]: H(3)/H(6) (2.2~,), and H(10)/H(21) (2.4 A). The global minimum conformation of the cis isomer, 2e.1, has the distances between protons H(3)/H(6) (2.1 A), H(3)/ H ( l l ) (3.3 A), while conformation 2c.2 is the only one with a short contact o f 2.2 ,~ for the pair H(8)/ H(11). The A M B E R * global minimum conformation
lc.1 has only one pair of protons in close proximity, H(3)/H(6) (2.1 A).
3.4. Force field dependence of the results. Semiempirical molecular orbital calculations of some low-energy conformations of lc, 2t and 2c All A M B E R * local minima of l c with relative energies less than 16.0 kJ mo1-1 were subsequently minimized also with the MM3* [16] force field parameters. The results confirmed the expectation for possible significant dependence of the relative energies of the conformations on the molecular mechanics force field used. A local minimum with A M B E R * relative energy 1 5 . 2 k J m o l -~ now has the lowest energy with MM3*. This conformation (designated l c . m m 3 in Table 4) and the conformations l e a and lc.2 were further calculated with the original MM3 version of Allinger's force field [8]. The same calculations were carried out also for 2t.1 and the four lowest energy conformations of 2c obtained with AMBER*. These results are compared in Table 4 with the relative energies of the same structures estimated with the semiempirical molecular orbital parametrizations AM1 [20] and PM3 [21], using M O P A C 93.0 [22]. MM3 and the semiempirical molecular orbital methods are consistent in predicting the conformation l c . m m 3 to have the lowest energy (Table 4). The three methods reverse the order of preference for
128
P.M. lvanov/Journal of Molecular Structure 440 (1998) 121-130
2t.1
2c.1
2c.2 Fig. 4. Stereographic views 118] of low-energy conformations of 2t and 2e.
P.M. lvanov/Journal of Molecular Structure 440 (1998) 121-130
conformations lc.1 and lc.mm3, computed with the AMBER* force field. The preference of lc.1 with AMBER* is determined by more favourable torsional and electrostatic contributions. A completely different balance of interactions is displayed by the MM3 force field. The angle bending and torsional terms have the decisive role in lowering the energy of conformation lc.mm3 when the MM3 force field is used. The computed vicinal coupling constants for pairs of protons in the macrocyclic ring of lc.mm3 are also in accordance with the corresponding measured couplings [1], and this conformation has the same short distances for pairs of protons as structure 2c.1. In accordance with the AMBER* results, MM3 and PM3 predicted the conformation 2t.1 to have lower energy than the four conformations of 2c considered. Only AM 1 gives a significantly different order of stability for the set of structures of 2t and 2c given in Table 4. In summary, the analysis of the conformational data obtained from the search for low-energy conformations of lc, 2t and 2c with the AMBER* force field, and the comparative studies of some local minima also with MM3, AM 1 and PM3, suggest that macrocyclic ring conformations having the hydrogen bonding patterns of lc.1 and lc.mm3 also have to be considered among the probable conformations of tylosin in solution.
4. Conclusions A Monte Carlo conformational search technique was used with the AMBER* force field as a model to search for low-energy conformations of tylosin, 1, and protynolide, 2. Selected optimized structures obtained from the search were also further studied with the MM3 force field [8] and semiempirical molecular orbital methods [20-22]. Structural and energy data were obtained for the 12-cis isomer of 1 (lc), and the 12-trans (2t) and 12-cis (2c) isomers of 2. The diene system at C(10), C(ll), C(12) and C(13) of 2 has the t r a n s - t r a n s configuration [2]. The NOE experiments show that the two C=C double bonds of 1 have the s - t r a n s arrangement [1 ]. The computed data for lc, 2t and 2c provide grounds for conclusions to be drawn about probable conformations of 1 and 2 in solution. Irrespective of the configuration around the C(12)=C(13) double bond, the 12-atoms fragment
129
of the 16-membered ring, not containing the diene system, has enough flexibility to adopt different low-energy hydrogen bonded conformations. The answers to the questions formulated at the beginning can be briefly summarized as follows. The geometry of the 16-membered ring of the AMBER* computed low-energy conformations of le, as well as the lowest-energy conformation of lt, is not the same as the geometry of the macrocycle in the crystal of 2 [2], nor does it resemble the geometry of the low-energy conformations of 2t and 2c that were computed with the same force field. Intramolecular hydrogen bonding is among the important factors which determine the preferred geometry of the macrocycle. The pattern of intramolecular hydrogen bonding is also different in the most probable conformations of 1 and 2 found with the AMBER* force field. The lowest energy conformation of 2t, conformation 2t.1, indeed corresponds to the conformation found for 2 in the crystal [2]. The computed relative energies of the conformations of the isomers of 1 and 2 strongly depend on the molecular mechanics force field used. Thus, in addition to the most probable conformation lc.1 of lc found from the Monte Carlo search with the AMBER* force field, we propose a conformation with the 16-ring conformation and hydrogen bonding pattern of lc.mm3 (Fig. 3), also to be considered as a possible candidate for a low-energy conformation of 1. Other conformations may also be found if other search procedures are applied [5], or if more advanced molecular mechanics models are used (see [23] for some recent developments). Testing different values for the arguments of the MCMM (Monte Carlo Multiple Minimum), TORS (variable TORSion selection) and RCA4 (Ring Closure Atoms(4)) commands [15,16], and/or adding also the TORC (TORsional Constraint) command [15,16] of the search procedure used in our study, could also influence the results. Such a systematic approach, requiring significant computational resources, is, however, beyond our reach at present.
Acknowledgements Financial support from the Direcci6 General de Universitats from Generalitat de Catalunya, through
130
P.M. Ivanov/Journal of Molecular Structure 440 (1998) 121-130
a grant for a visiting professorship at the Universitat Aut6noma de Barcelona, Spain (1995), is gratefully acknowledged. Thanks are due also to Prof. Carlos Jaime for providing computational resources, and to Prof. E. Osawa for facilities for preparing the manuscript.
References
[10] [11] [12] [13]
[1] S. Simova, G. Ivanova, Magn. Reson. Chem. 34 (1996) 255. [2] S. Omura, H. Matsubara, A. Nakagawa, A. Fukusaki, T. Matsumoto, J. Antibiot. 33 (1980) 915. [3] P.M. Ivanov, J. Mol. Struct., in press. [4] (a) P.M. Ivanov, C. Jaime, Anal. Quim. Int. Ed. 92 (1996) 13. (b) P.M. Ivanov, C. Jaime, J. Mol. Struct. 377 (1996) 137. [5] (a) H. Goto, E. Osawa, J. Mol. Struct. (Theochem) 285 (1993) 157. (b) H. Goto, E. Osawa, J. Chem. Soc., Perkin Trans. 2, (1993) 187. (c) F. Villamagna, M.A. Whitehead, J. Chem. Soc., Faraday Trans. 96 (1994) 47. (d) J. Koca, S. Perez, A. Imberty, J. Comput. Chem. 16 (1995) 296. (e) L. DosenMicovic, Tetrahedron 51 (1995) 6789. (f) S.B. Engelsen, J. Koca, I. Braccini, C.H. du Penhoat, S. Perez, Carbohydr. Res. 276 (1995) 1. (g) A. Smellie, S.D. Kahn, S.L. Teig, 25 (1995) 285. (h) J.C. Meza, R.S. Judson, T.R. Faulkner, A.M. Treasurywala, J. Comput. Chem. 17 (1996) 1142. (i) I. Kolossvary, W.C. Guida, J. Am. Chem. Soc. 118 (1996) 5011. (j) C.S. Wang, J. Comput. Chem. 18 (1997) 277. [6] (a) N.L. Allinger, J. Am. Chem. Soc. 99 (1977) 8127. (b) U. Burkert, N.L. Allinger, Molecular Mechanics, American Chemical Society, Washington, DC, 1982. [7] C. Jaime, E. Osawa, Tetrahedron 39 (1983) 2769. [8] (a) N.L. Allinger, Y.H. Yuh, J.H. Lii, J. Am. Chem. Soc. 111 (1989) 8551; 8566; 8576. (b) MM3(92): Technical Utility Corporation, 235 Glen Village Court, Powell, OH 43065, USA. [9] (a) P.M. Ivanov, J.M. Rudzinski, E. Osawa, results presented at the 13th Symposium on Chemical Information and Computer Sciences, 28-29 November 1990, Toyohashi,
[14]
[15]
[16]
[ 17] [18]
[19] [20] [21] [22]
Japan. (b) E. Osawa, P.M. Ivanov, J.M. Rudzinski, KESSHOU: Molecular mechanics program for molecular and crystal structure optimization, presented in part in: E. Osawa, Potential energy calculations of crystals, Ch. 1. I of Reactivity in Molecular Crystals, Y. Ohashi (Ed.), Kodansha, Tokyo and VCH, Weinheim, 1993, pp. 1-9. (c) P.M. Ivanov, J. Mol. Struct., in press. P.M. Ivanov, unpublished results, 1991. N. Norman, H. Mathisen, Acta Chem. Scand. 15 (1961) 1755. P.M. Ivanov, T.G. Momchilova, I.G. Pojarlieff, J. Comput. Chem. 12 (1991) 427. C.M. Breneman, K.B. Wiberg, J. Comput. Chem. 11 (1990) 361. (a) P.M. Ivanov, Jpn Chem. Program Exch. Newslett. 7(3) (1995) 3. (b) P.M. Ivanov, Jpn Chem. Program Exch. Newslett. 7(4) (1996) 13. (a) G. Chang, W.C. Guida, W.C. Still, J. Am. Chem. Soc. 111 (1989) 4379. (b) M. Saunders, K.N. Houk, Y.D. Wu, W.C. Still, M. Lipton, G. Chang, W.C. Guida, J. Am. Chem. Soc. 112 (1990) 1419. F. Mohamadi, N.G.J. Richards, W.C. Guida, R. Liskamp, C. Caufield, G. Chang, T. Hendrickson, W.C. Still, J. Comput. Chem. 11 (1990) 440. W.C. Still, A. Tempczyk, R.C. Hawley, T. Hendrickson, J. Am. Chem. Soc. 112 (1990) 6127. M.N. Burnett, C.K. Johnson, ORTEP-III: Oak Ridge thermal ellipsoid plot program for crystal structure illustrations, Oak Ridge National Laboratory Report ORNL-6895, 1996. C.A.G. Haasnoot, F.A.A.M. De Leeuw, C. Altona, Tetrahedron 36 (1980) 2783. M.J.S. Dewar, E.G. Zoebisch, E.F. Healy, J.J.P. Stewart, J. Am. Chem. Soc. 107 (1985) 3902. J.J.P. Stewart, J. Comput. Chem. 10 (1989) 221. MOPAC 93.00: J.J.P. Stewart, Fujitsu Ltd, Tokyo, Japan, 1993.
[23] (a) AMBER: W.D. Cornell, P. Cieplak, Ch.I. Bayly, I.R. Gould, K.M. Merz, Jr, D.M. Ferguson, D.C. Spellmeyer, Th. Fox, J.W. Caldwell, P.A. Kollman, J. Am. Chem. Soc. 117 (1995) 5179. (b) MM4: N.L. Allinger, K. Chen, J.H. Lii, J. Comput. Chem. 17 (1996) 642. (c) Merck Molecular Force Field: T.A. Halgren, J. Comput. Chem. 17 (1996) 616.