Molecular docking and dynamics simulations in the ligand binding domain of steroid hormone receptors

H. van der Goot (Editor) Trends in Drug Research III © 2002 Elsevier Science B.V. All rights reserved

57

Molecular Docking and Dynamics Simulations in the Ligand Binding Domain of Steroid Hormone Receptors Milou Kouwijzer-f- and Jordi Mestres§ i Molecular Design & Informatics, N.V. Organon, 5340 BH Oss, The Netherlands § Computational Medicinal Chemistry, Organon Laboratories Ltd., Newhouse, Lanarkshire MLl 5SH, Scotland

Introduction During the process of protein-ligand recognition, the structures of both the protein and the ligand undergo conformational changes as a result of a mutual induced fit aiming at optimising their interaction. Although ligand flexibility is now efficiently accounted for in docking calculations [1-8], proper treatment of protein flexibility still remains a critical issue for docking methods due to the complexity of sampling and scoring and the amount of computer power required [9-17]. The inability to account for protein flexibility can sometimes have a significant effect on the results of docking calculations. On one hand, it is known from crystal structures of proteins cocrystallized with different ligands that the backbone and/or the side chains of the residues in the binding pocket adapt to the ligand bound. This situation is illustrated in Figure la for two crystal structures of the estrogen receptor a (ERa). On the other hand, there is an intrinsic ambiguity in the orientation of heteroatoms for some residues (particularly threonine, histidine, glutamine and asparagine) due to the fact that differences in the electron density of carbon, nitrogen, and oxygen atoms are often difficult to appreciate. This situation is illustrated in Figure lb for the two independent molecules in the crystal structure of the progesterone receptor. It is thus clear that the result of a docking calculation depends on such details in the protein structure, and that scoring functions cannot correct for the inflexibility of the protein. One way of taking into account the flexibility of a protein is using molecular dynamics (MD) simulations. Because of the computer time requirements, it is not applicable for large numbers of ligands but it can be used for binding mode assessment activities. The aim of the present work is to evaluate the performance of MD simulations when applied to the ligand binding domain of steroid hormone receptors and compare the results with those obtained from flexible ligand docking methods that make use of a single rigid crystal structure.

58

Figure 1. Side chain adjustments and ambiguities in the binding pockets of steroid hormone receptors, a (left): superposition of ERa cocrystallized with estradiol [23] (light gray) or DBS [24] (dark gray), b (right): superposition of two PR structures [26] (molecule A in light gray and B in dark gray). The choice of a single rigid structure for docking purposes will clearly affect the results.

Methodology For the MD simulations, only the ligand and the residues that have at least one atom within 8 A of the ligand were allowed to move (i.e. -800 out of -4000 atoms). From a 100 ps run at 400 K (after a heating phase of 10 ps) coordinate sets were saved every picosecond. After the MD simulation, these coordinate sets were energy minimized (still taking into account the constraints) and the average interaction energy from these 100 frames was calculated. The orientation having the lowest average interaction energy was taken as the most favorable orientation. All calculations were performed with QUANTA/CHARMm [18] and charges were taken from templates. A 15 A cutoff radius was used, with a shift and a switch function between 11 and 14 A for the Coulomb and van der Waals energy, respectively. For the minimizations and MD runs, a distance-dependent dielectric constant 8 of 4r was used, whereas for the calculation of the interaction energy, 8 was set to r [19]. Docking calculations were performed with DOCK 4.0 [4]. Sphere centers were generated with SPHGEN [20]. For computational reasons, initial sphere sets were reduced to a number of around 30, ensuring appropriate coverage of the binding pocket. Energy (AMBER [21]) and chemical scoring grids of 0.3 A were generated using a united atom model and an 8 of 4r with a 10 A cut-off radius. Mol2 files of the ligands were generated with SYBYL [22] using Gasteiger charges. Virtual screening calculations were performed using uniform sampling with a maximum of 50 orientations and 5 seeds, whereas higher sampling values of 500 orientations and 25 seeds were set for binding mode

59 calculations. The minimal anchor size was set to four atoms and the maximum number of steric clashes between protein and ligand was set to three. Energy/chemical score minimization of maximally 100 iterations to a convergence of 0.1 kcal/mol was performed to each docked ligand. A database of 1000 diverse compounds selected from the Available Chemicals Directory (ACD) was used for virtual screening purposes. For validation purposes, MD simulations in the ligand-binding domains of steroid hormone receptors with their co-crystallized ligands were performed, namely, ERa__lere (estrogen receptor a, structure from pdb code lere) with estradiol [23], ERa_3erd (structure from pdb code 3erd) with diethylstilbestrol (DES) [24], ERp with genistein [25], and the progesterone receptor (PR) with progesterone [26]. In addition, MD simulations of ERa structures with the ligand that was co-crystallized in another structure were also performed, namely, ERa^lere with DES, and ERa_3erd with estradiol. The former calculation was repeated for 1000 ps at 300 K to validate the high temperature used in the simulations. The aim of these calculations was mainly to assess the effect and the extent of the accommodation of the binding pocket to a ligand other than the native co-crystallized ligand. Finally, as an application example, a prediction of the orientation of genistein in the ERa binding pocket based on MD simulations is presented (the crystal structure of the ERa-genistein complex has not yet been determined experimentally). Due to the high computer time requirements for the MD simulations, sampling of ligand orientations in the binding pocket was limited to the most probable orientations. The four starting orientations for the ligands containing a steroidal core are represented below, orientation 1 being the one observed in the crystal structures containing estradiol [23] and progesterone [26]. Analogous orientations were used for the non-steroidal ligands, DES and genistein.

1

2

3

4

Final structural details: in the ERo/p receptors, the proton of the (uncharged) histidine residue close to the steroidal D ring was placed on Ne2 and, in the PR structure (molecule A), the labels of atoms O E I and NE2 in Gln-725 were exchanged, in accordance with the hydrogen bonding scheme published [26] and the second molecule in the asymmetric unit.

60 Molecular docking In order to assess to which extent the use of a single rigid crystal structure affects the final results of a docking calculation, the structures of ERa cocrystallized with estradiol [23] and DES [24] were used in a virtual screening exercise. As shown above in Figure 1, the binding pocket of the two structures differs mainly in the conformation of four residues, namely, Met343, Thr347, Met421, and His524. While DES pushes away Met343 and His524 with respect to the position they adopt with estradiol, estradiol induces a significant conformational change to Met421 with respect to the position it adopts with DES. More interestingly, while in lere the hydroxy 1 of Thr347 points towards the ligand cavity, in 3erd it points away from the ligand cavity. All these conformational changes will define essentially a different accessible space and interaction pattem in the binding pocket and thus result potentially in a different ranking for the molecules in a database. The correlation between the energy and chemical score ranks obtained by docking 1000 diverse molecules from the ACD in the two rigid ERa structures (lere and 3erd) is presented in Figure 2. As can be visually observed, a better correlation between molecule ranks is found when using the chemical score (r^=0.82) than when using the energy score (r^=0.67). In this respect, the chemical score seems to be less dependent than the energy score on the conformational changes of the residues defining the binding pocket. This is a particularly interesting statement that would require though further investigation to other protein binding pockets to assess its generality. In conclusion, it has been shown that the use of different single rigid structures does have a significant effect in the final outcome of docking calculations and that efficient treatment of protein flexibility during the course of a docking calculation would be highly desired. BOO

'."^.s-S&l^'

••••••.• • • •••

700 600

••"••.Js^^A-

111 500 Ul •- 400

k'^-^"'0

100

200

300

'••• 300

400 3ER0

SOO

200

600

700

800

»30

100 0<

pr\. }

100

•': ••••-% • 200

300

•. 400

500

600

700

800

90

3ERD

Figure 2. Correlation between the energy score ranks (left) and the chemical score ranks (right) obtained by docking 1000 diverse molecules from the ACD in the ligand binding domains defined by the ERa crystal structures with PDB entries lere and 3erd.

61 Molecular dynamics simulations Simulations of the ERa ligand-binding domain as defined by the PDB entry 3erd crystal structure with estradiol, DES, and genistein were recently reported in the literature [27]. In the work presented here, MD simulations of the ERa, ERP, and PR ligand-binding domains were performed and the results are discussed in the remainder of the paper. The first validation test is to check whether the experimentally observed binding mode is found as the most favorable orientation from MD simulations of a protein crystal structure with its cocrystallized ligand. The results of this exercise are presented in Figure 3. It can be seen that in all cases orientation 1 (see above for orientation numbering), which corresponds to the orientation observed in the crystal structure, is indeed the most favorable. This shows that MD simulations are capable of discriminating the correct from alternative orientations of the ligand in the binding pocket. Only for DES in ERa__3erd, differences between alternative solutions are small, in agreement with the highly symmetrical nature of this ligand. For the other cases, orientation 1 is clearly preferred over alternative orientations. The next step is to check whether the experimentally observed binding mode is reproduced from MD simulations when a ligand is docked in the binding pocket of a non-native protein crystal structure. Figure 4 shows the results for the binding of estradiol and DES in their counterpart crystal structures. Again, the experimentally observed orientation is found to be the most favorable, although the amino acid side chains in the binding pocket were initially not optimally adapted to the ligands. This proves that the conformations of the residues in the binding pocket can indeed accommodate for other ligands than the one with which the crystal structure was determined. When the MD simulation for the ERa_lere complex with DES was run for 1000 ps at 300 K the results did not significantly differ from those obtained from that for 100 ps at 400 K.

62 Estradiol in ERa (1ere)

Progesterone in PR (1a28) 1

1

I

o

I LU V

.1

•45

-45

cd

1

•55

i

-55

-65

65

-75

7C

i

1

1

-

t 1

•

1

1 2 3 4 Genistein in ERp {1ql
1 2 3 4 Diethylstilbestrol in ERa (3erd)

O

E cd

LU V

1 2 3 4 orientation number

-75

1 2 3 4 orientation number

Figure 3. Average interaction energies obtained from MD simulations of different ligands in four orientations in the ligand-binding domain of steroid hormone receptors. A comparison of the results obtained from MD and DOCK calculations for the docking of estradiol and DES in their native and non-native protein structures (lere and 3erd) is presented in Table 1. RMSD values are calculated between the most favorable orientation found with MD or DOCK and the experimentally determined ligand orientation. As can be observed, results for the two native protein-ligand structures (estradiol in lere and DES in 3erd) are reasonably accurate from both methods (RMSD<0.60 A), MD giving slightly lower RMSD values than DOCK (0.05 A for estradiol and 0.28 A for DES). However, MD results for the two non-native protein-ligand structures (DES in lere and estradiol in 3erd) are significantly better than DOCK results (0.74 A for estradiol and 0.32 A for DES). The poorer results obtained by DOCK in the non-native cases are indicative of docking calculations using rigid protein structures.

63 Estradiol in ERa (3erd) 1

Diethylstilbestrol in ERa (lere)

1

-45

-45 h L

o

E Iri o c

UJ"

-55 -65

V 7 C

-

T

-if 1

1

-55 -65

1

-75

1 2 3 4 orientation number

1 2 3 4 orientation number

Figure 4. Average interaction energies from MD simulations of estradiol in ERa_3erd and DES in ERa_lere in the same four orientations as in Figure 3.

Table I. RMSD (A) for the non-hydrogen atoms between the optimal ligand orientation obtained from MD or DOCK and the experimentally determined ligand orientation. Plain and italic numbers correspond to native and non-native protein-ligand structures, respectively. ERa_.lere

ERa_3erd

MD

DOCK

MD

DOCK

Estradiol

0.52

0.57

0.44

0.76

DES

0.34

1.08

0.32

0.60

Simulations on the preferred orientation of genistein in the ligand-binding domain of ERa were also performed. In the crystal structure of genistein bound to ERP [25], helix 12 adopts the antagonistic conformation. In the present MD simulations, helix 12 in ERa was retained in the agonistic conformation (as in lere and 3erd). The basis for this assumption is twofold: first, genistein is a partial agonist in ERp but, given its rather small size, it is very well possible that helix 12 in ERp with genistein bound can adopt both agonistic and antagonistic conformations [25] and, second, although the affinity of genistein for ERa is lower than for ERp, it appears to be a full agonist (even a slight superagonist) in ERa [28].

64 The results of the MD simulations of genistein in four orientations in ERa are shown in Figure 5. As can be observed, the results are consistent with orientation 4 being the preferred orientation, irrespective of the ERa crystal structure used originally. However, interestingly enough, the predicted binding mode of genistein in ERa (orientation 4) differs from the experimentally observed binding mode in ERP (orientation 1), as illustrated in Figure 6. Since the crystal structure of the ERa-genistein complex is not available, the scope of this prediction cannot be established yet. However, since the results from both ERa_lere and ERa_3erd agree and since the same methodology applied to genistein in ERp gives indeed the correct result, we are confident that the present prediction may have a high degree of reliability.

Genistein in ERa (1ere) 1

1

1

Genistein in ERa (3erd) 1

1

-45

-45

Id -55

-55

lij" -65

-65

1

1

1

1

1

"o

g V

-75

1

1

1

1

1 2 3 4 oritentation number

-75

I

I

1 2 3 4 orientation number

Figure 5. Average interaction energies from MD simulations of genistein in ERa.

65

Figure 6. Left: experimentally determined orientation of genistein in ERp [25], including hydrogen bonding partners (Glu-305, Arg-346, a water molecule, and His-475). Right: predicted binding orientation in ERa (including hydrogen bonding partners Glu-353, Arg-394 and His-524; water was not included)

Conclusions One of the remaining problems with docking ligands in receptors is the flexibility of the protein. Standard docking programs do not take protein flexibility into account in an efficient manner yet. The dynamical calculations described here do take protein flexibility properly into account, as shown by the correct reproduction of the experimentally determined ligand binding modes in the crystal structures of steroid hormone receptors. This is also the case when a ligand is simulated in the binding pocket of a crystal structure that was cocrystallized with another ligand. Of course, these calculations take much more computer time than standard docking programs (one orientation takes about 12 hours on a SGI octane with a RIOOOO processor). The methodology is therefore not to be used for large numbers of ligands. However, it can be highly valuable in structure-based drug design problems, when one is interested in the binding mode of a small series of structurally diverse ligands for which the binding mode is not known. Finally, the relevance of the prediction that the binding mode of genistein in ERa is different than in ERP will have to await experimental confirmation.

References [1] M. Rarey, B. Kramer, T. Lengauer, G. Klebe, 7. MoL Biol. 261 (1996) 470. [2] W. Welch, J. Ruppert, A.N. Jain, Chem. Biol 3 (1996) 449. [3] G. Jones, P. Willett, R.C. Glen, A.R. Leach, R. Taylor, J, MoL Biol, 267 (1997) 727. [4] T. Ewing, I.D. Kuntz. /. Comput. Chem. 18 (1997) 1175.

66 [5] G.M. Morris, D.S. Goodsell, R.S. Halliday, R. Huey, W.E. Hart, R.K. Belew, A.J. Olson, J. Comput. Chem. 19 (1998) 1639. [6] S. Makino, I.D. Kuntz. J. Comput. Chem. 19 (1998) 1834. [7] C.A. Baxter, C.W. Murray, D.E. Clark, D.R. Westhead, M.D. Eldridge, Proteins 33 (1998) 367. [8] J. Wang, P.A. KoUman, I.D. Kuntz, Proteins 36 (1999) 1. [9] A.R. Leach, J. Mot. Biol. 235 (1994) 345. [10] R.M.A. Knegtel, I.D. Kuntz, CM. Oshiro, J. Mol. Biol. 266 (1997) 424. [11] J. Apostolakis, A. Pliickthun, A. Caflisch, J. Comput. Chem. 19 (1998) 21. [12] B. Sandak, H.J. Wolfson, R. Nussinov, Proteins 32 (1998) 159. [13] V. Schnecke, C.A. Swanson, E.D. Getzoff, J.A. Trainer, L.A. Kuhn, Proteins 33 (199S) 14. [14] M. Zacharias, H. Sklenar, J. Comput. Chem. 20 (1999) 287. [15] J.Y. Trosset, H.A. Scheraga, J. Comput. Chem. 20 (1999) 412. [16] M. Mangoni, D. Roccatano, A. Di Nola, Proteins 35 (1999) 153. [17] A.C. Anderson, R.H. O'Neil, T.S. Surti, R.M. Stroud, Chem. Biol. 8 (2001) 445. [18] CHARMm version 25.2, QUANTA 98, Molecular Simulations Inc., San Diego, USA. [19] P.D.J. Grootenhuis, P.J.M. van Galen, Acta Cryst. D 51 (1995) 560. [20] I.D. Kuntz, J.M. Blaney, S.J. Oatley, R. Langridge, T.E. Ferrin, J. Mol. Biol. 161 (1982) 269. [21] S.J. Weiner, P.A. KoUman, D.T. Nguyen, D.A. Case, J. Comput. Chem. 7 (1986) 230. [22] SYBYL version 6.7, Tripos Inc., St. Louis, USA. [23] A.M. Brzozowski, A.C.W. Pike, Z. Dauter, R.E. Hubbard, T. Bonn, O. Engstrom, L. Ohman, G.L. Greene, J.-A. Gustafsson, M. Carlquist, Nature 389(1997)753. [24] A.K. Shiau, D. Barstad, P.M. Loria, L. Cheng, P.J. Kushner, D.A. Agard, G.L. Greene, Cell 95 (1998) 927. [25] A.C. Pike, A.M. Brzozowski, R.E. Hubbard, T. Bonn, A.G. Thorsell, O. Engstrom, J. Ljunggren, J.-A. Gustafsson, M. Carlquist, EMBO J. 18 (1999)4608. [26] S.P. Williams, P.B. Sigler, Nature 393 (1998) 392. [27] B.C. Oostenbrink, J.W. Pitera, M.M.H. van Lipzig, J.H.N. Meerman, W.F. van Gunsteren, / Med. Chem. 43 (2000) 4594. [28] Barkhem T, Carlsson B, Nilsson Y, Enmark E, Gustafsson J & Nilsson S. Mol. Pharmacol. 54(1998) 105-112.