THEOCH 5283
Journal of Molecular Structure (Theochem) 419 (1997) 121–131
Quantum chemical modelling of thyroid hormone analogues Wiesław Nowak a,*, Andrzej Wojtczak b a
Institute of Physics, N. Copernicus University, ul. Grudzia¸dzka 5, 87-100 Torun´, Poland Faculty of Chemistry, N. Copernicus University, ul. Gagarina 7, 87-100 Torun´, Poland
b
Received 1 November 1996; accepted 11 February 1997
Abstract The thyroid hormone thyroxine (T4) and its analogues are important factors in regulating a range of physiological processes, among them tissue growth and differentiation. The interaction of the hormone with the transport protein transthyretin (TTR) has been the subject of X-ray studies (A. Wojtczak et al., Acta Cryst., D52 (1996) 758) and early molecular mechanics simulations (J. Blaney et al., J. Am. Chem. Soc., 104 (1982) 6424). In order to obtain a better insight into the mechanism of hormone recognition, we have performed quantum chemical studies of T4, 3,39-diiodothyronine, 39,59-dinitrothyronine, EMD21388, T4Ac and 3,39,59-triiodothyronine using standard ZDO (AM1, PM3) methods. For T4, the DFT method was used additionally. The structures of isolated molecules were optimized, and charge distributions and electron densities were calculated. The optimized structures are in reasonable agreement with the available X-ray data. Results of our calculations, especially related to electrostatics, will be used in simulations of the dynamics of the hormone–TTR complex. q 1997 Elsevier Science B.V. Keywords: Thyroid hormones; Thyroxine; AM1 method; PM3 method; Quantum chemical calculations; Electrostatic potential; Geometry optimization
1. Introduction The thyroid hormone thyroxine (T4) and the product of its monodeiodination, triiodo-l-thyronine (T3), regulate tissue growth and the differentiation, synthesis and metabolism of lipids and proteins as well as tissue oxygen consumption. The hormone is delivered to target receptors throughout the general circulation as a ligand of three transport proteins, among them transthyretin (TTR). While bound in TTR, the thyroxine core in a skewed conformation is positioned between aliphatic side chains of Ala and Leu. Its phenolic hydroxyl group interacts with the Ser side chains from two * Corresponding author. E-mail:
[email protected]
TTR subunits, and the alanyl moiety of the hormone forms a pair of salt bridges to Glu and Lys near the binding site entrance [1,2]. The protein crystal structures of TTR complexes have revealed alternative modes of thyroxine analogue binding [2–4]. Halogen-like substituents are positioned in the hydrophobic pockets between anti-parallel beta strands of TTR. T4 iodine substituents form alternative interactions to nucleophilic carbonyl oxygen atoms of the surrounding strands. The resulting alternative ˚ , and in some cases positions differ by almost 1 A the presence of additional water molecules mediating the binding interactions is suggested [3,5]. The position of the hormone analogue is affected by tyrosyl ring substitution and ether bridge flexibility. In this way, the mode of the hormone substitution affects
0166-1280/97/$17.00 q 1997 Elsevier Science B.V. All rights reserved. PII S 0 16 6- 1 28 0 (9 7 )0 0 25 1 -0
W. Nowak, A. Wojtczak/Journal of Molecular Structure (Theochem) 419 (1997) 121–131
Fig. 1. Compounds studied.
122
W. Nowak, A. Wojtczak/Journal of Molecular Structure (Theochem) 419 (1997) 121–131
the ability to form polar interactions to Glu and Lys, differentiating the binding affinity. Genetic disorders causing abnormalities in the hormone binding to transport proteins result in hypoor hyperthyroxinemia and severe disfunction of the organisms of affected patients. A molecular level understanding of ligand recognition and binding specificity is a prerequisite condition for an effective treatment. Our research on the binding of thyroxine and its analogues to TTR also indicated the possibility of multiple ligand binding modes [2]. In order to understand the molecular basis of ligand recognition in transthyretin, we have initiated the molecular modelling of ligands and TTR. Firstly, we examined the rigidity of the thyroxine ether bridge using the molecular dynamics technique with the CHARMM force field [6]. Secondly, the minimized structures and molecular electrostatic potentials of T4, 3,5,39,59-tetraiodothyroacetic acid (T4Ac), 3,39,59triiodothyronine (rT3), 3,39-diiodothyronine (T2), 39,59-dinitro-N-acetylthyronine (DNNAT), and 39,59dibromo-3-methyl-6,49-dihydroxyflavone (EMD21388) have been calculated using standard quantum-chemical methods (mainly AM1 as implemented in the InsightII package [7]). Results for both X-ray and minimized structures were compared, and the effects
123
of the conformations on the electrical properties of the molecules were studied. These data are compared with those of the AMBER force field obtained by Kollman and co-workers [8,9]. Our results will be used in simulations of alternative binding modes of T4 in TTR tetramers.
2. Methods The AM1 method was developed for structural and energetic studies of large organic molecules [10]. Despite its known shortcomings [11], it often provides useful and reliable information on the geometries of flexible systems (see, for example, Ref. [12]). In the present paper, the AM1 method has been applied for the first time for the modelling of thyroxine hormone analogues presented in Fig. 1. All initial structures obtained from the scratch (‘‘Builder’’ module of InsightII [7]) EMD, T4Ac, rT3 or crystal data were subject to complete energy minimization using the EF method of the mopac6 package [13] as implemented in the InsightII v.3.5 software [7]. In all cases the PRECISE option was used, and the frequencies in the minimized structures were calculated to exclude transition states.
Fig. 2. Atom numbers and definitions of torsional angles.
124
W. Nowak, A. Wojtczak/Journal of Molecular Structure (Theochem) 419 (1997) 121–131
Fig. 3. AM1 calculated bond lengths in EMD.
T4 was also subject to classical CHARMM force field minimization using the X-PLOR package [14]. Additionally, the DFT DMOL code [7] was used for calculations of the electrostatic potential of the T4 molecule in the AM1 optimized conformation. The molecules were also subject to geometry optimization using the PM3 method [15]. Unless otherwise indicated, all partial charges were calculated using Mulliken population analysis. The electrostatic potential (ESP) charges were obtained using the Merz–Kollman method of fitting the electrostatic potential on the Williams surfaces of the studied molecules [16]. Fig. 5. Charge distribution in T4.
3. Results 3.1. Structure of the molecules The AM1 minimized structures are presented in the form of mopac archive files which may be obtained as described in Section 6. The XMOL software allows for interactive studies of details of the geometry. These structures are also
shown in Figs. 3–10, together with the calculated AM1 (Mulliken) charge distributions. The values of the torsional degrees of freedom that are most important for the general shape of the hormones are presented in Table 1. Standard definitions of angles are shown in Fig. 2, and the terminology used for discussions of the conformations is the same as that used in Ref. [3]. Selected geometry parameters obtained from the theoretical calculations are compared with the X-ray and literature data presented in Table 1. The AM1 geometry of the EMD analogue is shown in Fig. 3 (bond lengths) and in Fig. 4 (bond angles). 3.2. Charge distributions
Fig. 4. AM1 calculated bond angles in EMD.
There are numerous ways of obtaining the charge distribution in organic molecules using quantum chemical methods. In Figs. 5–10 we present selected charges calculated within the AM1 method for the
125
W. Nowak, A. Wojtczak/Journal of Molecular Structure (Theochem) 419 (1997) 121–131 Table 1 Comparison of the geometry of thyroxine analogues a Parameter
f (deg) f9 (deg) x 1 (deg) x 2 (deg) w (deg) C–I phe ˚) (A C–I tyr ˚) (A a (deg)
T4
rT3 Protein, ˚c 2.0 A
AMB
T4Ac
AM1
e
AM1
e
T2
DNNAT
AM1
PM3
77.5 25.2 48.6 74.8 −162.3 2.02
83.5 13.0 43.9 76.7 −174.3 1.97
86.5 14.1 29.8 82.6 158.8 2.09
108.3 −29.2 −50.0 162.2 −159.3 2.10
102.9 −13.7 −144.7 97.9 93.8 2.10
90.0 0.0 – – – 2.08
83.3 25.7 47.0 73.7 −162.6 2.02
105.0 −19.9 −61.9 50.9 – 2.02
22.6 59.5 −57.2 107.4 172.6 2.02
2.02
1.97
2.09
2.09
2.08
2.08
2.02
2.02
2.02
117.1
121.3
Crystal
d
AM1
116.2
XPL
b
120.4
123.6
125.0
115.9
116.2
116.5
Protein, ˚f 2.0 A 93.8 19.2 −20.6 115.8 171.0 2.20 g
AM1 39.1 43.0 52.3 96.9 165.4 –
94.8 65.3 −49.8 96.0 124.1 –
–
–
116.2
118.2
2.07 127.0
Protein, ˚c 2.2 A
a
AM1 denotes the AM1 minimized structure; XPL denotes the X-PLOR minimized structure [14], CHARMM force field; and AMB denotes the AMBER force field minimized structure. Definition of torsional angles (see Fig. 2): f, C3–C4–O4–C19; f9, C4–O4–C19–C29; x 1, N8– C8–C7–C1; x 2, C8–C7–C1–C2; w, N8–C8–C9–O10. a, C4–O4–C19 Bond angle; C–I phe, phenolic ring carbon–iodine bond length; C–I tyr, tyrosyl ring carbon–iodine bond length. b Ref. [1]. c Ref. [2]. d Ref. [8]]. e Crystal structures in TTR or molecular crystals unknown. f Ref. [4]. g ˚ resolution. Deformation in the crystallographic refinement of 2.0 A
Fig. 6. Charge distribution in T4Ac.
Fig. 7. Charge distribution in rT3.
126
W. Nowak, A. Wojtczak/Journal of Molecular Structure (Theochem) 419 (1997) 121–131
Fig. 10. AM1 charge distribution in EMD.
by Blaney et al. [8] in their AMBER modelling for T4 are also shown.
4. Discussion Fig. 8. Charge distribution in T2.
4.1. Structure of the molecules minimized structures. All AM1 charges for each compound may be extracted from the mopac *.arc files. The charges relating to the most important atoms and functional groups of the analogues are presented in Table 2. In this table the DFT charges obtained with the BLYP functional (DMOL code [7]) and those used
Fig. 9. Charge distribution in DNNAT.
The theoretical prediction of the geometries of large, flexible molecules is a non-trivial task. A comparison of calculated and experimental X-ray structures gives useful hints as to whether the semiempirical AM1 method can be used in further modelling of the thyroxine hormones. Thyroxine is a derivative of the amino acid tyrosine. Therefore the geometry of T4 and its analogues is discussed in terms of torsional angles generally used to describe the conformation of the polypeptide chain and amino acid side chain conformations. Below, a short discussion of the optimized structures is presented for each compound separately. 4.1.1. T4 (thyroxine) Both crystal structures, X-PLOR, AMBER, PM3 and AM1 minimized structures have a skewed conformation of the ether bridge close to the ideal (908,08). This conformation might be expected to be predominant due to the steric efects of tyrosyl iodine atoms. The tyrosyl x 1 angle is found to be + gauche in AM1, PM3 and X-PLOR, − gauche in the T4 crystal structure, and −1458 in the TTR–T4 complex, reflecting one of the three most probable conformations of the side chain. The tyrosyl x 2 value is about 908 in the AM1 and X-PLOR minimization, as well as in the protein complex, while in the T4 crystal structure it
127
W. Nowak, A. Wojtczak/Journal of Molecular Structure (Theochem) 419 (1997) 121–131 Table 2 Calculated partial charges on selected atoms Atom (group)
T4 AM1
T4 PM3 a
T4 DFT b
T4 DFT c
T4 AMBER d
T4Ac
rT3
T2
DNNAT
EMD
O4 (ether) O49 (phenol) H (phenol) I39 (phenol) I59 (phenol) I3 (tyrosyl) I5 (tyrosyl) O10 (carboxy) O11 (carboxy) N8 (amino) H (amino) H (amino) H (amino)
−0.152 −0.230 0.234 0.162 0.185 0.179 0.200 −0.454 −0.526 −0.060 0.203 0.283 0.230
−0.109 −0.195 0.203 0.058 0.082 0.084 0.116 −0.460 −0.577 0.675 −0.011 0.013 0.115
−0.104 −0.173 0.121 0.060 0.091 0.066 0.103 −0.310 −0.331 −0.006 0.153 0.173 0.138
−0.496 −0.636 0.503 0.037 0.058 0.056 0.081 −0.474 −0.494 −1.119 0.491 0.574 0.503
−0.24 −0.40 0.40 −0.07 −0.07 −0.07 −0.07 −0.543 −0.543 −0.364 0.332 0.332 0.332
−0.149 −0.231 0.234 0.160 0.184 0.182 0.188 −0.357 −0.311 – – – –
−0.158 −0.231 0.234 0.161 0.184 0.192 – −0.453 −0.528 −0.061 0.231 0.284 0.232
−0.144 −0.233 0.227 0.183 – 0.185 – −0.443 −0.535 −0.062 0.282 0.280 0.224
−0.145 −0.199 0.277 −0.125 – – – – – – – – –
– −0.223 0.239 0.095(Br) 0.075(Br) – – – – – – – –
a
PM3 optimized geometry. AM1 geometry used, Hirshfeld partition method. c AM1 geometry used, Mulliken partition method. d AMBER parameters used, data from Ref. [8]. b
is 1628, the latter probably being due to packing forces in the crystal lattice. The angle corresponding to the main chain w angle is found to be −1608 in both AM1 and T4 crystal structures and 1608 in X-PLOR, the conformation from the allowed but not the most favoured region of the Ramachandran plot. In the protein complex, w = 948 corresponds to the b conformation.
The difference reflects the contribution of two salt bridges formed by the T4 alanyl moiety to Glu-54 and Lys-15 in the protein environment of TTR. ˚) The C–I distances in the AM1 structure (2.02 A are slightly shorter than those in the crystal structures ˚ ). An ether bridge or predicted by X-PLOR (2.09 A angle of about 1208 is similar for all methods. The
Fig. 11. A comparison of AM1 and X-PLOR optimized structures of T4 with X-ray molecular crystal data and X-ray data for the T4–TTR complex.
128
W. Nowak, A. Wojtczak/Journal of Molecular Structure (Theochem) 419 (1997) 121–131
T4 structures obtained from the T4 crystals, crystals of the complexes of T4 with TTR, X-PLOR and the AM1 methods are compared in Fig. 11. 4.1.2. T4Ac (tetraiodothyroacetic acid) The ether bridge conformation found using the AM1 method (skewed) is consistent throughout the whole series. The tyrosyl moiety conformation −608/ 508 is the most frequent of those found in proteins. The x 2 value of 508 is slightly smaller than the expected 908 value, but it is still in an acceptable range for a crowded protein environment with tight packing interactions. 4.1.3. rT3 (3,39,59-triiodothyronine) The AM1 results gave a skewed conformation of the ether bridge (268/838) consistent with the crystal structure of T4 [1] and the ligand conformation found in TTR complexes of 3,39-T2 [4] and T4 [2]. Such a conformation is caused by the restricted flexibility of the ether bridge due to steric effects of the bulky iodine substituents. The tyrosyl moiety conformation is + gauche (478/748) and is not the most favourable for the side chains in proteins. However, the same conformation is calculated with AM1 for T4 and DNNAT. This energy minimum is allowed for an isolated amino acid with no steric hindrance from the adjacent residue in the polypeptide chain. The AM1 value of the analogue of the main chain w angle of −162.68 is not favoured in proteins. However, it is in an acceptable range. 4.1.4. T2 (3,39-diiodothyronine) ˚ and the ether The phenolic C–I distance of 2.20 A bridge angle of 1278 result from deformations in the ˚ resolution crystallographic refinement with 2.0 A PROLSQ. The ether bridge conformation in the TTR complex structure is skewed (908,08), while for the isolated molecule minimization with AM1 gives a structure between anti-skewed and twisted. This reflects the increased conformational flexibility for the analogue with the monosubstituted tyrosyl ring. The conformations of the tyrosyl moiety in the protein complex and in the AM1 minimized structure are similar (−gauche). An w angle of 1708 is identical in both the crystal structure and the AM1 results.
4.1.5. DNNAT (3,39-dinitro-N-acetylthyronine) The ether bridge conformation (958,658) is found to be close to the perpendicular (908,908) in the TTR complex. The AM1 calculations result in a twisted structure (458,458). The difference is probably due to protein binding effects compared to isolated molecule minimization. The different conformations reflect the increased flexibility of the bridge in the hormone analogue not substituted in the tyrosyl ring. In this analogue, the tyrosyl side chain conformations 528/ 978 (AM1) and −508/968 (protein complex) correspond to two preferred conformations in the proteins. The AM1 structure is much less probable, and may suggest that a local minimum in the rotational space was found. The w angles of 1658 in AM1 and 1248 in the protein complex are similar to those observed in T4. The difference seems to be a systematic effect of the isolated molecule minimization. The ether bridge angle is almost identical in the theoretical calculation and in the crystal structure. 4.1.6. EMD21388 The most important bond lengths and angles are presented in Figs. 3 and 4 respectively. Selected geometrical AM1 data are compared with the EMD crystal structure [17] in Table 3. For the numbering of atoms used in this table, Fig. 12 should be consulted. As one can see, the AM1 minimized geometry is in very good agreement with the protein crystal structure; however the AM1 C–Br bond lengths are 0.026 Table 3 A comparison of selected AM1 internal coordinates of EMD21388 with the X-ray data a Coordinate
AM1
X-ray
˚) C4–Br5 (A ˚) C9–Br10 (A ˚) C1–C27 (A ˚) C19–O20 (A ˚) C6–O7 (A ˚) C15–O13 (A ˚) C21–O22 (A ˚) C27–O22 (A ˚) C27–C28 (A ˚) C28–C29 (A
1.873 1.875 1.467 1.239 1.367 1.377 1.385 1.388 1.358 1.480 43.93 117.35 109.54
1.887 1.901 1.497 1.234 1.350 1.360 1.378 1.378 1.331 1.500 45.31 118.92 107.90
O22–C27–C1–C2 (deg) C21–O22–C27 (deg) O22–C27–C1 (deg) a
For numbering of the atoms, see Fig. 12.
W. Nowak, A. Wojtczak/Journal of Molecular Structure (Theochem) 419 (1997) 121–131
129
Fig. 12. The numbering of atoms in EMD.
˚ shorter than those observed in the and 0.014 A crystal. The torsional angle of 45.88 is close to the 43.98 angle observed in the crystal, which suggests that the crystal packing does not involve any special strain imposed on this torsion. From a comparison of the structures, we conclude that the AM1 model gives a reasonable description of the geometry of the EMD analogue. 4.2. Charge distribution Specific electrostatic interactions may be responsible for the negative cooperativity observed in the TTR tetramers. Therefore it is very important to study the effects of different structures and conformations of the ‘‘tail’’ tyrosyl group or its equivalents on the charge distribution in the ‘‘head’’ phenol ring moiety. The most important partial charges are collected in Table 2. In all iodinated analogues the AM1 charges on I39 and I59 are similar (+0.16 and +0.18) and do not depend on details of the tyrosyl group stucture. A slight asymmetry between I39 and I59 results from the interaction of I59 with the hydrogen atom of the hydroxyl group. The charges of the O49 (hydroxyl group) and the corresponding H atom are also very similar in all the derivatives, including EMD. A slightly higher charges of the H atom in DNNAT indicate a greater acidity of the –OH group in this analogue. The PM3 method gives lower positive charges on the iodine atoms (+0.05, +0.08) but very similar charges to the AM1 method on the hydroxyl group. The AM1 charges on the bromine atoms in EMD are smaller (0.095 and 0.075) than those on the iodine atoms, so we expect that the halogen interactions of EMD with protein nucleophilic groups should be smaller than the interactions of iodinated compounds. Since in the first modelling of T4 interactions
Fig. 13. A comparison of AM1 electron density distribution in T4 (black) and DNNAT (grey).
with proteins, AMBER parameters with negative values of the charges located on all iodine atoms have been used (see Table 2), perhaps the electrostatic interactions were not correctly represented in that model. The DFT partial charges obtained for T4 using Mulliken population analysis are quite different from those obtained based on the Hirshfeld method, the later being more consistent with those obtained by semiempirical calculations. It is worthwhile to note that in all zwitterionic molecules (T4, T2 and rT3) we observe a clear asymmetry of the AM1/PM3 charge distribution on the oxygen atoms in the carboxy group. This observation indicates that O10 and O11 are not equivalent with respect to their hydrogen-bonding abilities since the ‘‘cis’’ orientation of the COO − and NH +3 groups leads to the preferential interaction of one of the oxygen atoms with the amino group. The total charges calculated for the NH +3 group are identical in all charged compounds, but the PM3 and AM1 methods give different distributions of the charges on the N and H atoms. The sum of the partial charges of the atoms belonging to the nitro group of DNNAT is negative (−0.125) while all calculated partial charges for iodine atoms are positive. Since the size of the iodine substituent
130
W. Nowak, A. Wojtczak/Journal of Molecular Structure (Theochem) 419 (1997) 121–131
Fig. 14. Electrostatic potential charges based on the AM1 charge distribution in DNNAT.
and the nitro group is similar (in terms of the van der Waals radius), the detailed differences in the electron density distribution must be responsible for the different binding properties of T4 (T2) and
DNNAT. Indeed, the electron density contoured at the same level reveals relatively lower values on iodine atoms compared with the nitro group substituent (see Fig. 13 below). This confirms the role of protein nucleophilic group interactions with ligand halogens in ligand recognition and effective binding in TTR. The calculated distribution of the electron density is consistent with the observation based on the protein crystal structures that iodine atoms of T4 interact with nucleophilic carbonyl groups of the polypeptide main chains, while similarly positioned nitro substituents of DNNAT hydrogen bond to the serine side chain hydroxyl groups. The electrostatic potential derived charges, based on the AM1 electron density distribution, gave reasonable results only for DNNAT, i.e. the only analogue without iodine or bromine substituents. The selected ESP charges are presented in Fig. 14 and all values may be obtained from the specially edited mopac archive file dnnatesp.arc from the WWW version of the paper. In the other cases, the ESP InsightII option led to very high (1.5–3.5) values for the point charges located on the phenol ring carbon atoms. In our opinion, these results of the fitting procedure indicate a deficiency in the standard point charge model of electrostatics in iodinated hormones. Perhaps higher moments should be included or dummy potential sites should be introduced to describe correctly lone pairs of iodine and bromine atoms [16]. The high asymmetry of the electrostatic potential in the vicinity of the phenol ring is demonstrated in Fig. 15. There is also the possibility that we have encountered a case in which the standard least-squares fitting procedure is inadequate (see, for example, Ref. [18]). Thus the problem of ESP charges in iodinated derivatives requires further study.
5. Conclusions
Fig. 15. Electrostatic potential contours (in hartrees) based on the DFT charge distribution in T4. 0.001, grey; 0.010, black; dots, van der Waals radii.
The AM1 method is suitable for structural studies of analogues of thyroxine hormones. Minimized structures compare well with the available X-ray data. Calculated electron density distributions indicate basically different types of interactions with
W. Nowak, A. Wojtczak/Journal of Molecular Structure (Theochem) 419 (1997) 121–131
protein receptors in DNNAT and iodinated analogues. The AM1 charges, despite their ‘‘reasonable’’ values, should perhaps be used only with caution in classical MD modelling of these hormones. Large gradients of molecular electrostatic potentials near the iodinated phenol group lead to ‘‘exotic’’ values of ESP AM1 charges. The DFT study of tyroxine analogues, which is now in progress in our laboratory, should give a better-grounded model of electrostatics in these systems.
6. Supplementary material The AM1 minimized structures are presented in the form of mopac archive files that may be viewed using, for example, the XMOL code from the location of the electronic version of this paper. http://www.phys.uni.torun.pl/~wiesiek/ECCC3/ Paper 5l.html
Acknowledgements Support from the Polish State Committee for Scientific Research, grant no. 6 P04A 032 11, project BiMol (FNP) and UMK grants no. 328-F and 389-F, is acknowledged. The authors also thank Mr. K. Wejer for his assistance in HTML editing and Mr. G. Bakalarski for DFT calculations.
131
References [1] V. Cody, Acta Cryst., Sect. B 37 (1981) 1685. [2] A. Wojtczak, V. Cody, J. Luft, W. Pangborn, Acta Cryst., Sect. D 52 (1996) 758. [3] P. de la Paz, J.M. Burridge, S.J. Oatley, C.C.F. Blake, in: C.R. Beddell (Ed.), The Design of Drugs to Macromolecular Targets, J. Wiley & Sons Ltd., London, 1992, pp. 119-172. [4] A. Wojtczak, J. Luft, V. Cody, J. Biol. Chem. 267 (1992) 353. [5] C.C.F. Blake, S.J. Oatley, Nature 268 (1977) 115. [6] W. Nowak, (1995) unpublished results. [7] InsightII v.2.3.7, Biosym Technologies, Inc., USA, 1995. [8] J. Blaney, P.K. Weiner, A. Dearing, P.A. Kollman, E.C. Jorgensen, S.J. Oatley, J.M. Burridge, C.F.F. Blake, J. Am. Chem. Soc. 104 (1982) 6424. [9] T.A. Andrea, S.W. Dietrich, W.J. Murray, P.A. Kollman, E.C. Jorgensen, S. Rothenberg, J. Med. Chem., 22 (1979) 221. J.M. Blaney, E.C. Jorgensen, M.L. Connolly, T.E. Ferrin, R. Langridge, S.J. Oatley, J.M. Burridge, C.F. Blake, J. Med. Chem., 25 (1982) 785. [10] M.J.S. Dewar, E.G. Zoebisch, E.F. Healy, J.J.P. Stewart, J. Am. Chem. Soc. 107 (1985) 3902. [11] J.P. Stewart, in: K.B. Lipkowitz, D.B. Boyd (Eds.), Reviews in Computational Chemistry, Vol. 3, VCH Publishers, Inc., 1991, New York, p. 45. [12] W. Nowak, M. Wierzbowska, J. Mol. Struct. (Theochem) 386 (1996) 223. [13] J.J.P. Stewart, J. Comput.-Aided Mol. Des. 4 (1990) 1. [14] A. Brunger, X-PLOR Manual, Yale University Press., 1992. [15] J.J.P. Stewart, J. Comput. Chem. 10 (1989) 209. [16] D.E. Williams, in: K.B. Lipkowitz, D.B. Boyd (Eds.), Reviews in Computational Chemistry, Vol. 3, VCH Publishers, Inc., 1991, New York, pp. 219-271. [17] V. Cody, J. Luft, M. McCourt, K. Irmscher, Struct. Chem. 2 (1991) 601. [18] M.M. Francl, Ch. Carey, L.M. Chirlian, D.M. Gange, J. Comput. Chem. 17 (1996) 367.