Insights into Ligand Selectivity in Estrogen Receptor Isoforms

J. Phys. Chem. B 2008, 112, 2719-2726

2719

Insights into Ligand Selectivity in Estrogen Receptor Isoforms: Molecular Dynamics Simulations and Binding Free Energy Calculations Juan Zeng,† Weihua Li,† Yaxue Zhao,† Guixia Liu,*,† Yun Tang,*,† and Hualiang Jiang†,‡ Laboratory of Molecular Modeling and Design, School of Pharmacy, East China UniVersity of Science and Technology, Shanghai, 200237, People’s Republic of China, and Shanghai Institute of Materia Medica, Chinese Academy of Sciences, Shanghai, 201203, People’s Republic of China ReceiVed: October 15, 2007; In Final Form: December 7, 2007

Estrogen receptor R (ERR) and β (ERβ) are ligand activated transcription factors that have different physiological functions and differential expression in certain tissues. The ligand binding domain of ERβ shares 58% sequence identity with that of ERR. However, in the binding pocket there are only two relative residue substitutions. This high similarity at the active site is a great challenge for designing selective estrogen receptor modulators. ERβ is shown to be related to several diseases. To understand the molecular basis of ERβ selectivity, molecular dynamics simulations were carried out for both ERR and ERβ complexes. Our simulations revealed the conformational changes at the active site of the ERs and the difference of affinity with ligand. The electrostatic repulsion between the Sδ atom of ERR Met421 and the acetonitrile group nitrogen atom of the ligand led to unfavorable binding. The repulsion resulted in the conformational change of the side chain of ERR Met421, which changed the conformation of both Leu346 and Phe425. These residues changes expanded the volume of binding cavity, which led to unstable binding of the ligand. In addition, ERβ Met336 formed more hydrophobic contacts with the ligand relative to corresponding residue ERR Leu384. Furthermore, the binding free energy analysis was shown to be correlated with the previous results determined by experiment. At last, free energy decomposition evidently indicated the contributions of key residues. The present results could help explain the mechanism of ERβ selectivity and may be considered in the design of subtype-selective ligands.

1. Introduction Estrogens are hormones that have many important physiological and pharmacological activities that include critical roles in the development and function of the female reproductive system, in the maintenance of bone mineral density,1,2 and in the control of blood lipid levels.3,4 Estrogens exert most of their biological effects via interaction with the estrogen receptor (ER), a member of the nuclear hormone receptor superfamily. To date, it is known that there are two ERs, ERR and ERβ. The first estrogen receptor (now called ERR) was cloned in 1986,5 and at that time it appeared that only one ER existed. However, the discovery of a second estrogen receptor subtype (ERβ) in 1996 led to an intense interest in elucidating ERβ function6,7 and determining which aspects of estrogen biology it mediated. The expression pattern of the two subtypes differs. Although widely found in many tissues, ERβ is predominately expressed in ovarian granulosa cells, lung, bladder, and prostate, whereas ERR is the predominant ER expressed in the uterus, kidney, and ovarian theca cells.8-10 A recent work, using propylpyrazole triol, the ERR selective agonist as efficacious as estradiol, suggested that many of the classical effects of estrogens were mediated via ERR.11 What about ERβ? There is a pressing need for ERβ selective compounds for new therapies, as well as understanding of the * Corresponding authors. Tel: +86-21-64251052. Fax: +86-2164253651. E-mail: (G.L.) [email protected]; (Y.T.) ytang234@ ecust.edu.cn. † East China University of Science and Technology. ‡ Chinese Academy of Sciences.

physiological function of ERβ. Fortunately, a highly selective ERβ agonist ERB-041 was discovered recently, and with the identification of it it was discovered that ERβ might exert its function in inflammatory diseases such as intestinal and joint inflammation.12,13 Moreover, selective ERβ agonist was recently found to have novel activity in rodent benign prostatic hyperplasia models14 and to be able to inhibit breast cancer cell proliferation.15 ERR and ERβ have distinct functions and differential expression in certain tissues. These differences have stimulated the search for subtype-selective ligands. ER isoforms share a common topology with a striking degree of conservation in the active site. The ligand binding domain of ERR and ERβ structures has 58% sequence identity.16 Actually, there are only two amino acid differences among the residues that are in close proximity to the ligand in the binding pocket with Met336 in ERβ being replaced by Leu384 in ERR and Ile373 in ERβ being replaced by Met421 in ERR. Although this slight difference in the binding cavity poses a great challenge in developing ER subtype-selective ligands, a great number of compounds of selective ERβ agonists with impressive activity have been found, such as benzoxazole/benzofunan,17 benzopyran,14,18 biphenyls,19 indazoles,20 Chromenoquinoline,21 tetrahydrochrysenes,22,23 and diarylpropionitriles.24 Furthermore, several crystal structures of human ERR and ERβ isoforms with various classes of ligands have been determined.14,17,23,25,26 These findings are very useful for identifying the structural variations and their implication on the design of selective ligands. On the basis of structural and biochemical data, the mechanism of selectivity of ERβ has been postulated.17,23,27-29

10.1021/jp710029r CCC: $40.75 © 2008 American Chemical Society Published on Web 02/12/2008

2720 J. Phys. Chem. B, Vol. 112, No. 9, 2008

Figure 1. Chemical structure of WAY-244.

In the present study, to clarify the mechanism of ligand selectivity, we carried out molecular dynamics simulations of ERR and ERβ complexed with WAY-244 (see Figure 1, a selective ERβ ligand with IC50 values of 14.0 nM for ERβ and 1062.5 nM for ERR17) and subsequent binding free energy analysis as well as the calculation of free energy decomposition. Our simulations indicated that the ERβ selectivity of the ligand attributed to not only its interactions with the only two variable residues at the active site, Met336 in ERβ and Met421 in ERR, but also the influence of environmental residues at the active pocket. Our results further shed light on its selectivity and aid in the development of selective and potent ligands of the two isoforms. 2. Materials and Methods 2.1. Molecular Dynamics Simulations. Before carrying out molecular dynamics (MD) simulations, quantum chemical calculations were executed for WAY-244 to derive the atom charges utilized in MD simulations. After geometry optimization, the electrostatic potential of WAY-244 was preliminarily calculated at the B3LYP/6-31G(d,p) level using the Gaussian03 program.30 The atomic charges were then fitted to the electrostatic potentials using the restrained electrostatic potential fit procedure,31 as implemented in the AMBER 9 antechamber module.32 The simulation models used in the present study were constructed based on the X-ray crystal structures of WAY-244 complexed with ERβ and ERR (PDB code: 1X78 and 1X7E respectively).17 The missing atoms were added by using the molecular modeling software Sybyl 7.0.33 One crystallization water molecule remained, which was involved in a hydrogenbonding network between ERβ residues Glu305 and Arg346 (ERR residues Glu353 and Arg394) and the ligand, and other water molecules were removed. All the minimizations and MD simulations were performed with Sander module of AMBER 9 package.34 Each initial structure for WAY-244 complexed with ERβ and ERR was modeled using LeaP module. The AMBER ff99 force field35 was used as the parameters for van der Waals and bonded energy terms basically, and the general AMBER force field36 was used as the parameters for WAY-244. Hydrogen atoms were then added to the system. All ionizable residues were set at their default protonation states at a neutral pH. To ensure overall electroneutrality of the system, an appropriate number of counterions were added. The entire system was subsequently surrounded with a box of TIP3P water molecules37 with a margin of 10.0 Å along each dimension. All covalent bonds to hydrogen atoms were constrained using the SHAKE algorithm.38 The Particle Mesh Ewald method39 was applied to calculate longrange electrostatic interactions. The cutoff distance for the longrange electrostatic and the van der Waals energy terms was set at 10.0 Å. Periodic boundary conditions were applied to avoid edge effects in all calculations. The solvated models were first minimized via two steps. At first, harmonic constraints were applied to the complex with a

Zeng et al. strength of 500 kcal/(mol‚Å2). And then, all atoms were allowed to move freely. In each step, energy minimization was executed by the steepest descent method for the first 500 steps and the conjugated gradient method for the subsequent 4500 steps. Afterward, the system was heated from 0 to 300 K in 20 ps and equilibrated in 16 ps with the protein and ligand coordinates being constrained until the total energy of the model reached a stable state. Subsequently, the whole solvated model was subjected to a free simulation for 40 ps. Finally, a production run for 5 ns was performed using NPT ensemble at a constant temperature of 300 K and a constant pressure of 1 atm. 2.2. Binding Free Energy Analysis. MM-PB(GB)SA approach,40-43 which is implemented in the AMBER program, was applied to compute the binding free energy (∆Gbind) between the ligand and each protein. At first, a total of 100 snapshots were taken from the last 1 ns trajectory with an interval of 10 ps. All counterions and water molecules were stripped. For each snapshot, the protein-ligand binding free energy was calculated using the following equation:

∆Gbind ) Gcomplex - [Gprotein + Gligand]

(1)

where Gcomplex, Gprotein, and Gligand are the free energies of the complex, the protein, and the ligand, respectively. Each of them can be estimated as follows:

∆G ) ∆GMM + ∆Gsol - T∆S

(2)

where ∆GMM is the molecular mechanics free energy, ∆Gsol is the solvation free energy, and - T∆S represents the entropy term. The molecular mechanics free energy is calculated by the equation:

∆GMM ) ∆Gele + ∆Gvdw

(3)

where ∆Gele and ∆Gvdw represent the electrostatic and van der Waals interactions, respectively. The solvation free energy is composed of two components

∆Gsol ) ∆Gele,sol + ∆Gnonpol,sol

(4)

∆Gele,sol is the polar contribution to solvation, and ∆Gnonpol,sol is the nonpolar solvation term. The former could be obtained by solving the Poisson-Boltzmann equation44 for MM-PBSA (molecular mechanics Poisson-Boltzmann surface area) method or generalized Born models45 for MM-GBSA (molecular mechanics-generalized Born surface area), whereas the latter term was determined using

∆Gnonpol,sol ) γSASA + b

(5)

where γ, representing the surface tension, and b, being a constant, were set to 0.0072 kcal/(mol‚Å2) and 0, respectively, as in the work of Still and co-workers.46 SASA is the solventaccessible surface area (Å2) determined using the linear combination of pairwise overlaps model.47 In this study, we calculated the binding free energies of both complexes using the MM-PBSA method. With regard to entropic effects, structural reorganization and solvent entropy are important contributions to enzyme-ligand binding. Normalmode analysis is useful to estimate entropic changes of the solute molecule, but the calculation is considered problematic and timeconsuming.48,49 Moreover, this approach does not take the solvent entropy into consideration. For the above reasons, the entropy contribution to the binding free energy was neglected in this study. To obtain a detailed view of the protein-ligand binding, the interaction energies were decomposed to binding free energies

Insights into Ligand Selectivity in Estrogen Receptor Isoforms

J. Phys. Chem. B, Vol. 112, No. 9, 2008 2721

Figure 2. Time dependence of the RMSD of protein backbone atoms (a) and all atoms (b). The RMSD of ERR and ERβ complex is shown in black and red, respectively.

on a single residue. MM-GBSA was used for this purpose. Free energy decompositions were for molecular mechanics and solvation energies but not for entropies. The details of the method, which are the same as those of MM-PBSA, are described above. 3. Results and Discussion 3.1. Overall Structural Changes. To assess the dynamic stability of two estrogen receptor-ligand complexes during the molecular dynamics simulation, root-mean-square deviation (RMSD) values of protein backbone atoms on ERR and ERβ relative to the initial X-ray crystal structure along the entire MD trajectory were monitored, as shown in Figure 2a. In both cases, the complex structure appeared to have reached a stable state after 4 ns equilibration, where the RMSD value converged to a value around 1.7 Å. For each system, the time series of the RMSD of all atoms was also calculated and shown in Figure 2b. The protein atoms did not significantly deviate from the starting structure after 4 ns of the simulation. Therefore, our subsequent energy analyses in each case were based on the MD trajectory truncated between 4 and 5 ns. To investigate the dynamic features in ligand binding, rootmean-square fluctuation (RMSF) was also monitored along the MD trajectories of both ERR and ERβ complexes. The RMSF reflects the mobility of a certain residue around its average position. For both systems, the RMSF values of each residue derived from the 4-5 ns window are given in Figure 3. The results indicated that large fluctuations of residues occurred in loops connecting secondary structure elements. Such residues belonged to the surface regions of the protein, which were exposed to the solvent. Hence, it was not surprising that the region had a large RMSF value. 3.2. Structural Changes at the Active Site. To analyze conformational changes at the active site, the average structures of both estrogen receptor-ligand complexes from the equilibration were compared with the crystal structures. The average structures were obtained by averaging 100 snapshots, which were taken from the last 1 ns on the MD trajectory with an interval of 10 ps. The average structures of both systems showed

Figure 3. Average fluctuation of each residue observed in MD simulation for ERR complex (a) and ERβ complex (b).

high similarity in protein backbones. Both average structures were fitted to the crystal structures based on the coordinates of the main chain atoms N, CR, and C. The ERR complex was found to exhibit a slightly larger RMSD value of 1.64 Å, whereas the RMSD value for ERβ complex was 1.61 Å. Figure 4a,b shows the superposition of the crystal structure with the average structure. Furthermore, to conveniently compare the structural differences at the active site in both ERR and ERβ complexes, the average structures from the equilibration of the two systems were superimposed, which is shown in Figure 4c. 3.2.1. Conformational Changes of Two Variable Residues. Figure 4 obviously shows that the structure of the ERβ complex was maintained very well. In contrast, there was larger structural deviation from the crystal structure for the ERR model. The average structure of ERR complex superimposed on the crystal structure is shown in Figure 4a. There was a particularly large difference in orientation of the side chain of Met421. Moreover, a rotation occurred at the methylthio group of the side chain. The rotation resulted from the electrostatic repulsion between the acetonitrile group N atom of the ligand and the Sδ atom of Met421. It was the electrostatic repulsion that resulted in the side chain of Met421 turning its conformation in the active site to allow the acetonitrile group of the ligand to align optimally. However, the replacement of ERR Met421 by ERβ Ile373 could eliminate the repulsive interaction, because isoleucine was nonpolar and shorter and occupied less volume than methionine.27 To relieve the unfavorable repulsive interaction, the Sδ atom of ERR Met421 moved away from the acetonitrile group. Compared with ERR Met421, there were no large conforma-


Zeng et al.

Figure 4. Average structures from the last 1 ns of the MD trajectories of the complexes superimposed on the X-ray structures (a) ERR complex and (b) ERβ complex. And both average structures are superimposed (c). Heavy atoms of the ligands and selected neighboring key residues are shown. For panels a and b, the residues of the crystal structures are shown in magenta; the residues of MD structures are shown in atom type. For panel c, ERR residues are shown in purple, while ERβ residues are shown in atom type.

Figure 5. The distance between the center of mass of ERR Leu384 (red)/ERβ Met336 (black) and B-ring of the ligand.

tional changes seen in the corresponding residue ERβ Ile373 whose side chain slightly swung during the whole MD simulations. As for the other variable amino acid, ERR Leu384/ERβ Met336, to further know its dynamic change several key distances were monitored during the MD simulations. The plot of the distance between the center of mass of ERR Leu384/ ERβ Met336 and the B-ring of the ligand is shown in Figure 5. As seen in Figure 5, the plots of both isoforms fluctuated significantly before 1700 ps. After 1700 ps, however, the distance between the center of mass of ERβ Met336/ERR Leu384 and the B-ring of the ligand was apparently much closer with their average value of about 6.3 and 7.4 Å, respectively. It was obvious that ERβ Met336 formed more van der Waals contacts and stronger interactions with the ligand compared to ERR Leu384. ERR Leu384 and ERβ Met336 are both hydrophobic amino acids. Yet, compared with ERR Leu384, ERβ Met336 contacted with the ligand with a stronger hydrophobic interaction; additionally, Methionine-aromatic interactions were thought to stabilize the protein.50-52 So it was possible that the interaction of the B-ring with ERβ Met336 relative to ERR Leu384 could lead to improved ERβ selectivity, which was in accord with previous studies.17,24,29 3.2.2. Conformational Changes of Other Residues. Conformational changes were not only observed at the two variable amino acids, but also other amino acids in the active site. As viewed in Figure 4a, the methyl group of Met421 remarkably moved toward the acetonitrile group of the ligand. However, the orientation of methyl group interfered with the interaction

Figure 6. The distance between the N atom of the acetonitrile group of the ligand and the center of mass of the benzene ring of ERR Phe425 (red) and ERβ Phe377 (black).

between Phe425 and the ligand. Phe425 had to rearrange in response to the steric clash, and then its side chain significantly moved away from the active site. As to the ERβ complex, the conformation of Phe377 of the average structure presented only a little slight deviation from that of the crystal structure. To compare the difference between the ERR and ERβ complex, we monitored the dynamic changes during the molecular dynamics simulation. The distance between the acetonitrile group N atom of the ligand and the center of mass of the benzene ring of ERR Phe425/ERβ Phe377 with time throughout the simulation is shown in Figure 6. As seen in Figure 6, in general the distance between the N atom and the center of mass of the benzene ring in the ERR complex was much larger than that in the ERβ complex. Moreover, the distance fluctuation of the ERβ complex was much smaller than that of the ERR complex. On the other hand, in the case of the ERR model, the orientation of the methyl group of Met421 could also influence the interactions between the ligand and Leu346. To reduce steric conflicts with the side chain of Met421, the entire structure of Leu346 modestly moved away from the ligand. In comparison, as for ERβ model, the side chain of Leu298 moved toward the ligand. Consequently, ERβ Leu298 exhibited more hydrophobic interactions with the ligand compared with ERR Leu346, and ERβ Leu298 had a relatively high affinity with the ligand. In addition, one could see from Figure 4a that the side chain of ERR Phe404 moderately shifted in the same direction as Phe425 shifted. This conformational change might be resulted from the



TABLE 1: Hydrogen Bonds in Each Modela hydrogen bond donor ERR

direct hydrogen bonds water-mediated hydrogen bonds

ERβ

direct hydrogen bonds water-mediated hydrogen bonds

ligand ligand Arg394 Wat8543 Arg394 ligand ligand Arg346 Wat243 Arg346 Wat243

O28 O12 NH2 O NH2 O28 O12 NH2 O NH2 O

acceptor H6 H11 HH22 H2 HH22 H6 H11 HH22 H1 HH22 H2

Glu353 His524 ligand Glu353 Wat8543 Glu305 His475 ligand Glu305 Wat243 Glu305

OE2 ND1 O28 OE1 O OE2 ND1 O28 OE1 O OE2

distance (Å)b

occupancy (%)c

2.576 2.912 3.002 2.732 3.079 2.598 2.929 3.061 2.711 3.011 3.172

100.0 99.1 71.2 169.8 89.0 100.0 99.1 45.9 93.0 79.80 52.6

a The listed donor and acceptor pairs satisfy the criteria for the hydrogen bond over 40.0% of the time during the last 1 ns of simulation. bThe average distance between hydrogen acceptor atom and hydrogen donor atom in the investigated time period. cOccupancy is in unit of percentage of the investigated time period.

switch of the ligand in the binding cavity. The shift of the side chain of Phe404 occurred to refrain from stereospecific clash with the ligand. In contrast, in ERβ model the corresponding residue Phe356 almost maintained its conformation. The structural difference between ERR Phe404 and ERβ Phe356 could be also observed in Figure 4c. As a result of these residues changes, the local binding pocket formed by residues Leu346, Leu402, Phe404, Met421, Ile424, Phe425, and Leu428 was expanded wider, which might lead to the unsteadily binding with the ligand. To make better contacts, the acetonitrile group of the ligand exerted its effort to fit this hydrophobic binding pocket, which induced the whole structure of the ligand fluctuated in the active site. All the above interactions and differences led to unfavorable binding with the ligand for the ERR model, and further contributed to higher selectivity of the ligand toward ERβ over ERR. In addition, according to previous studies17 conservative amino acids created hydrogen bonds with the ligand, which played a crucial role in stabilizing the binding conformation. Thus, we examined the snapshots during the last 1 ns of the simulation. The formation of a hydrogen bond was defined in terms of distance and orientation. The combination of donor D, hydrogen H, and acceptor A atoms with a DsH...A configuration was regarded as a hydrogen bond when the distance between donor D and acceptor A was shorter than 3.5 Å and the angle AsHsD was larger than 120.0°. Table 1 lists direct or one water molecule-mediated hydrogen bonds between the estrogen receptor and the ligand in each model with occupancy more than 40.0%. ERR His524 and ERβ His475 solidly formed a hydrogen bond with the ligand, respectively. Besides, in both ERR and ERβ systems, one water molecule mediated a hydrogen-bonding network among ERR residues Glu353 and Arg394 (ERβ residues Glu305 and Arg346) and the ligand. The average structure of this water molecule could be observed in Figure 4. For both models, this water molecule was fully present over the entire MD simulation process. 3.3. Conformational Changes of Ligand in ERR and ERβ Models. Apart from analysis of structural changes at the active site, conformational changes of the ligand were also examined to compare the difference of the interactions between the ligand and both isomers of ERs. According to average structure and snapshots from the MD trajectories, the ligand in the ERR complex underwent large conformational changes at the active site relative to that in the ERβ complex during the entire simulation, which could be seen in Figure 4. In addition, to assess the deviation of the ligand relative to their initial structure, the RMSD values during the whole MD simulation were calculated. Figure 7 shows RMSD changes of the ligand in both

Figure 7. RMSD values of all ligand atoms in ERR complex (shown in red) and ERβ complex (shown in black).

systems with respect to simulation time. It could be observed from Figure 4a that the conformation of the ligand in the ERR complex compared with the crystal structure significantly changed. It was primarily due to the shift of the plane of the ligand. In contrast, the extent of such shift in the ERβ complex was relatively small (see Figure 4b). In the case of the ERR complex, the entire structure moved to some extent relative to the crystal structure. This conformational mobility of the ligand in the ERR complex was suggested to cause unstable binding at the binding pocket. The conformational change was probably due to the rotation of the side chain of Phe425, which principally altered the shape of the local binding pocket formed by side chains of residues Leu346, Leu402, Phe404, Met421, Ile424, Phe425, and Leu428. Moreover, the ligand moved outward from the side chain of ERR Leu384, which resulted in reduced van der Waals contacts, while in the ERβ complex the ligand could interact with the side chain of Met336 via more hydrophobic contacts. At the same time, for the ERR complex the orientation of A-ring of the ligand led to unfavorable steric collisions with Phe404. Thus, the side chain of Phe404 moved away from the ligand, and the conformational change might decrease the interactions between the ligand and Phe404. His524 moved along with the ligand to form a stable hydrogen bond with the ligand. Additionally, regarding the ligand in the ERR complex, a comparison of the structure of the ligand in the crystal structure and that in simulation structure showed that another large conformational change occurred at the acetonitrile group, including changes in position, orientation, and conformation. This was probably due to the electrostatic repulsion between the Sδ atom of Met421 and the acetonitrile


Zeng et al.

Figure 8. Structure from 3660 ps (a) and 4790 ps (b) of the MD trajectories of the ERR model superimposed on the X-ray structures. Heavy atoms of the ligands and selected neighboring key residues are shown. For crystal structures, the residues are shown in magenta; for MD structures, the residues are shown in atom type.

N atom of the ligand. In the presence of Met421, the acetonitrile was restricted to some extent and unable to align itself optimally. So the conformation of the ligand changed in turn to alleviate the repulsion between the ligand and Met421. On the contrary, this electrostatic repulsion did not occur between the ligand and the nonpolar residue ERβ Ile373. It could be concluded from the analysis above that the conformational change of the ligand could lead to a significant decrease in ERR binding affinity relative to ERβ. The plot of RMSD changes of the ligand in both systems (Figure 7) showed that there existed larger RMSD values for the ligand in ERR complex throughout the MD simulation, which suggested that the binding of the ligand with ERR was less stable during the 5 ns simulation than that with ERβ due to large conformational changes. As for the ERR complex, the RMSD value of the ligand rose sharply from 3 to 5 ns, with the increment of about 0.4 to 0.6 Å. To further explore the protein-ligand interactions, two specific snapshots with large RMSD, namely, the structures at 3660 and 4790 ps among the MD trajectory, were extracted. And comparisons between these two structures and the corresponding crystal structures were carried out. Both of the MD simulation structures were superimposed to the crystal structure with the coordinates of the main chain atoms N, CR, and C (see Figure 8). By comparing the structures of each ligand in the MD structure with that in the crystal structure (Figure 8), it was found that the plane moieties in the ERR complex (Figure 8a) seemed almost the same as those of the average structure (Figure 4a), while acetonitrile group atoms in both structures were all found to exhibit large conformational changes. Compared to the crystal structure, the acetonitrile group N atom was situated at the position where it was distant from ERR Met421 Sδ atom, and the distance between acetonitrile group N atom of the ligand and ERR Met421 Sδ atom was 6.14 for 3660 ps structure and 6.39 Å for 4790 ps structure, respectively. This conformational change could still be ascribed to the repulsive interaction between the Met421 side chain and the ligand acetonitrile group. Therefore, from the snapshots of MD trajectory we could conclude that large RMSD values in the ligand were probably due to the conformational change of the acetonitrile group. With this question, the dihedral angle formed by C2, C3, C13, and C14 in the ligand (see Figure 1) was monitored along the entire MD trajectory of the ERR model. Figure 9 shows the dihedral angle C2-C3-C13-C14 of the ligand versus simulation time in the MD trajectory. Through close investigation, it was observed that the curve tendency of dihedral angle was consistent with the change in RMSD values, which indicated

Figure 9. Dihedral angle C2-C3-C13-C14 of WAY-244 versus simulation time in the MD trajectory. The atomic numbering of WAY244 is shown in Figure 1.

that large conformational change of the ligand in the ERR complex mainly occurred at the acetonitrile group. 3.4. Binding Free Energy Analysis. To compare the binding affinity of the ligand with both isoforms, binding free energy (∆Gbind) was calculated using MM-PBSA method. Table 2 shows the results, including the overall binding free energy and all of the energy terms given by the MM-PBSA method for each compound. According to Table 2, a higher value of ∆Gbind for the ERR complex than that for the ERβ complex corresponded to a decrease in the binding affinity. Furthermore, the difference in the calculated binding free energy between the ERR and ERβ complexes was 2.74 kcal/mol. Particularly, the biochemical data17 indicated that the IC50 value of WAY-244 in the ERR complex was about 75.8 times higher than that in the ERβ complex. For direct comparison with calculated affinities, ∆G was estimated approximately from IC50 values via ∆G ≈ -RTlnIC50. Accordingly, the difference in free energy between the ERR and ERβ complexes was 2.57 kcal/mol. Consequently, our energetic results that suggested the ERβ model was more stable than the ERR model were greatly compatible with that determined by experiment. Because entropic contribution was excluded, the calculated binding free energy was lower than that derived from the experiment. Nevertheless, it was still of significance to compare their relative magnitude. A detailed analysis suggested that, major contributions favorable to binding were van der Waals and electrostatic energies, whereas polar solvation energies opposed the binding. On the other hand, nonpolar solvation energies, which were



TABLE 2: Binding Free Energy Calculations for ERr and ERβ Complexesa ERR ERβ

∆Gele

∆Gvdw

∆Gele,sol

∆Gnonpol,sol

∆GMM

∆Gsol

∆Gbind

∆Gbind,exp

-48.96 ( 3.83 -40.85 ( 4.94

-32.84 ( 3.19 -37.77 ( 2.84

64.65 ( 4.78 58.68 ( 8.96

-4.90 ( 0.09 -4.85 ( 0.11

-81.80 -78.62

59.75 53.83

-22.05 ( 4.63 -24.79 ( 8.62

-7.93 -10.74

a ∆Gele: electrostatic interaction. ∆Gvdw: van der Waals interaction. ∆Gele,sol: polar contribution to solvation. ∆Gnonpol,sol: nonpolar contribution to solvation. ∆GMM: ∆Gele + ∆Gvdw. ∆Gsol: ∆Gele,sol + ∆Gnonpol,sol. ∆Gbind: binding free energy in the absence of entropic contribution. ∆Gbind,exp: 17 experimental data given as IC50 for both ligands. For direct comparison to calculated affinities, conversion to ∆G estimated by ∆G ≈ -RT ln IC50. All energies are in kcal/mol.

TABLE 3: Energy Contributions of ERr and ERβ Residues to the Binding of WAY-244a

a

ER

residue

∆Gele

∆Gvdw

∆Gele,sol

∆Gnonpol,sol

∆Gbind

R β R β R β R β R β R β R β R β R β

Leu346 Leu298 Glu353 Glu305 Leu384 Met336 Arg394 Arg346 Phe404 Phe356 Met421 Ile373 Phe425 Phe377 His524 His475 ligand ligand

-1.31 ( 0.27 -1.35 ( 0.30 -12.64 ( 1.37 -12.44 ( 1.41 0.06 ( 0.06 0.17 ( 0.15 -3.05 ( 1.02 -1.49 ( 1.60 -0.36 ( 0.32 -0.76 ( 0.18 0.01 ( 0.35 0.28 ( 0.14 -0.15 ( 0.23 -0.25 ( 0.27 -4.61 ( 0.80 -4.24 ( 0.79 -24.48 ( 1.92 -20.43 ( 2.47

-2.08 ( 0.35 -2.40 ( 0.39 1.23 ( 0.97 1.12 ( 0.88 -0.51 ( 0.18 -1.44 ( 0.25 -0.11 ( 0.26 -0.17 ( 0.25 -1.55 ( 0.37 -1.74 ( 0.31 -1.15 ( 0.35 -1.01 ( 0.27 -0.37 ( 0.17 -0.45 ( 0.15 -0.26 ( 0.63 -0.39 ( 0.65 -16.42 ( 1.60 -18.88 ( 1.42

1.15 ( 0.22 1.01 ( 0.21 10.00 ( 0.75 9.44 ( 0.98 -0.07 ( 0.06 0.01 ( 0.12 1.69 ( 0.69 0.57 ( 1.08 0.82 ( 0.33 1.00 ( 0.15 0.11 ( 0.31 -0.24 ( 0.11 0.19 ( 0.18 0.28 ( 0.21 2.75 ( 0.26 2.59 ( 0.28 34.17 ( 1.02 31.45 ( 1.54

-0.20 ( 0.05 -0.23 ( 0.04 -0.02 ( 0.01 -0.22 ( 0.01 -0.08 ( 0.03 -0.16 ( 0.03 -0.02 ( 0.01 -0.01 ( 0.01 -0.15 ( 0.04 -0.16 ( 0.03 -0.11 ( 0.03 -0.08 ( 0.04 -0.01 ( 0.01 -0.02 ( 0.01 -0.06 ( 0.02 -0.06 ( 0.02 -3.46 ( 0.04 -3.47 ( 0.04

-2.44 ( 0.50 -2.97 ( 0.41 -1.43 ( 0.43 -1.91 ( 0.58 -0.61 ( 0.18 -1.43 ( 0.25 -1.49 ( 0.39 -1.10 ( 0.49 -1.24 ( 0.35 -1.66 ( 0.31 -1.14 ( 0.30 -1.04 ( 0.28 -0.33 ( 0.22 -0.44 ( 0.16 -2.18 ( 0.45 -2.10 ( 0.47 -10.19 ( 1.37 -11.33 ( 1.50

All energies are in kcal/mol.

responsible for the burial of SASA upon binding, contributed slightly favorably. In both models, hydrogen bond and hydrophobic interactions were primarily considered to contribute to the binding energy.17 Hence, the results of our simulations were consistent with the results of previous experimental analyses. One could see from Table 2 that the electrostatic energy in the ERR complex was 8.11 kcal/mol less than that of the ERβ complex. However, the favorable electrostatic interactions were counteracted by the solvation effect. After offset, the difference in the electrostatic interactions decreased to 2.14 kcal/mol. The van der Waals interaction energies between WAY-244 and the ERR/ERβ were -32.84 and -37.77 kcal/mol, respectively. Obviously, the ERR model caused a large increase of the van der Waals energy. The calculated difference in van der Waals interaction could partly explain why the conformational change of the ligand of ERR model at binding pocket led to its reduced binding affinity. 3.5. Free Energy Decomposition. To further examine the protein-ligand interactions and the contribution of each residue, the binding free energy in the ERR and ERβ complexes was decomposed on per residue with MM-GBSA approach. Table 3 lists the energy contributions of certain structurally important residues for the protein-ligand interactions in both complexes. As can be seen, the residues with most favorable contributions to the binding free energy in both systems were residues Leu298/ 346 and His475/524. It was observed that their contribution to the binding free energy was more than -2 kcal/mol. Special attention had been paid to those residues with relatively large difference of the contributions to binding free energies, such as ERR Leu346/ERβ Leu298, ERR Leu384/ERβ Met336, and ERR Phe404/ERβ Phe356. For both ERR Leu346/ERβ Leu298 and ERR Leu384/ERβ Met336, the contribution of van der Waals energy was mainly responsible for the difference. And for ERR Phe404/ERβ Phe356, both electrostatic energy and van der Waals energy participated in yielding the energetic difference. From these results, we could conclude that hydrophobic

interactions played an important role in the binding affinity with ERβ over ERR. The contribution of ligand segments to binding was larger for the ERβ model (-11.33 kcal/mol) than that for the ERR model (-10.19 kcal/mol), which was in accord with the analysis above that the conformational change of ligand of the ERR model at binding pocket led to its reduced binding affinity. 4. Conclusions Molecular dynamics simulations and free energy calculations were performed on both ERR-WAY-244 and ERβ-WAY-244 complexes. MD simulation results showed that it was the interaction of the acetonitrile group of the ligand with Met421 in ERR corresponding to Ile373 in ERβ that caused the major conformational difference between these two isoforms. The electrostatic repulsive interaction between the acetonitrile group nitrogen atom of the ligand and the Sδ atom of ERR Met421 resulted in large conformational change of the side chain of Met421 itself with methyl group pointing to the ligand. To alleviate steric conflicts, the side chain of Phe425 significantly moved outward from the active pocket. In addition, the orientation of the side chain of Met421 also caused residue Leu346 to moderately move outward the binding cavity. Furthermore, due to the conformational changes of these residues, the active pocket was expanded, which may lead to the binding of the ligand in ERR to be less stable than that in ERβ. Division of the binding free energy into individual components showed that van der Waals and electrostatic terms were favorable to binding, whereas polar solvation opposed binding for both systems. The result by energy decomposition revealed that hydrophobic interactions primarily contributed to the ERβ selectivity over ERR. All these findings can be utilized in the design of new potent selective ERβ ligands. Acknowledgment. This work was supported by National Natural Science Foundation of China (Grants 20572023 and

2726 J. Phys. Chem. B, Vol. 112, No. 9, 2008 30600785), Shanghai Rising-Star Program (Grant 07QA14016), Shanghai Pujiang Program (Grant 05PJ14034), Shanghai Key Basic Research Project (Grant 05JC14092), and Key 863 HighTech Program (Grant 2006AA020404). References and Notes (1) Westerlind, K. C.; Wronski, T. J.; Ritman, E. L.; Luo, Z. P.; An, K. N.; Bell, N. H.; Turner, R. T. Proc. Natl. Acad. Sci. U.S.A. 1997, 94, 4199-4204. (2) Rapuri, P. B.; Gallagher, J. C.; Haynatzki, G. J. Clin. Endocrinol. Metab. 2004, 89, 4954-4962. (3) Mendelsohn, M. E.; Karas, R. H. N. Engl. J. Med. 1999, 340, 18011811. (4) Mendelsohn, M. E. Am. J. Cardiol. 2002, 89, 12E-17E; discussion 17E-18E. (5) Green, S.; Walter, P.; Kumar, V.; Krust, A.; Bornert, J. M.; Argos, P.; Chambon, P. Nature 1986, 320, 134-139. (6) Kuiper, G. G. J. M.; Enmark, E.; Pelto-Huikko, M.; Nilsson, S.; Gustafsson, J. A. Proc. Natl. Acad. Sci. U.S.A. 1996, 93, 5925-5930. (7) Mosselman, S.; Polman, J.; Dijkema, R. FEBS Lett. 1996, 392, 49-53. (8) Couse, J. F.; Lindzey, J.; Grandien, K.; Gustafsson, J. A.; Korach, K. S. Endocrinology 1997, 138, 4613-4621. (9) Fitzpatrick, S. L.; Funkhouser, J. M.; Sindoni, D. M.; Stevis, P. E.; Deecher, D. C.; Bapat, A. R.; Merchenthaler, I.; Frail, D. E. Endocrinology 1999, 140, 2581-2591. (10) Kuiper, G. G. J. M.; Carlsson, B.; Grandian, K.; Enmark, E.; Haggblad, J.; Nilsson, S.; Gustafsson, J. A. Endocrinology 1997, 138, 863870. (11) Harris, H. A.; Katzenellenbogen, J. A.; Katzenellenbogen, B. S. Endocrinology 2002, 143, 4172. (12) Harris, H. A.; Albert, L. M.; Leathurby, Y.; Malamas, M. S.; Mewshaw, R. E.; Miller, C. P.; Kharode, Y. P.; Marzolf, J.; Komm, B. S.; Winneker, R. C.; Frail, D. E.; Henderson, R. A.; Zhu, Y.; Keith, Jr. J. C. Endocrinology 2003, 144, 4241-4249. (13) Malamas, M. S.; Manas, E. S.; McDevitt, R. E.; Gunawan, I.; Zhang, X. B.; Collini, M. D.; Miller, C. P.; Dinh, T.; Henderson, R. A.; Keith, J. C., Jr.; Harris, H. A. J. Med. Chem. 2004, 47, 5021-5040. (14) Norman, B. H.; Dodge, J. A.; Richardson, T. I.; Borromeo, P. S.; Lugar, C. W.; Jones, S. A.; Chen, K. Y.; Wang, Y.; Durst, G. L.; Barr, R. J.; Rafizadeh, C. M.; Osborne, H. E.; Amos, R. M.; Guo, S.; Boodhoo, A.; Krishnan, V. J. Med. Chem. 2006, 49, 6155-6157. (15) Cvoro, A.; Paruthiyil, S.; Jones, J. O.; Tzagarakis-Foster, C.; Clegg, N. J.; Tatomer, D.; Medina, R. T.; Tagliaferri, M.; Schaufele, F.; Scanlan, T. S.; Diamond, M. I.; Cohen, I.; Leitman, D. C. Endocrinology 2007, 148, 538-547. (16) Pike, A. C. W.; Brzozowski, A. M.; Hubbard, R. E.; Bonn, T.; Thorsell, A.-G.; Engstrom, O.; Ljunggren, J.; Ohman, L.; Gustafsson, J. A.; Carlquist, M. EMBO J. 1999, 18, 4608-4618. (17) Manas, E. S.; Unwalla, R. J.; Xu, Z. B.; Malamas, M. S.; Miller, C. P.; Harris, H. A.; Hsiao, C. L.; Akopian, T.; Hum, W. T.; Malakian, K.; Wolfrom, S.; Bapat, A.; Bhat, R. A.; Stahl, M. L.; Somers, W. S.; Alvarez, J. C. J. Am. Chem. Soc. 2004, 126, 15106-15119. (18) Richardson, T. I.; Norman, B. H.; Luqar, C. W.; Jones, S. A.; Wang, Y.; Durbin, J. D.; Krishnan, V.; Dodge, J. A. Bioorg. Med. Chem. Lett. 2007, 17, 3570-3574. (19) Edsall, J. R. J.; Harris, H. A.; Manas, E. S.; Mewshaw, R. E. Bioorg. Med. Chem. 2003, 11, 3457-3474. (20) Angelis, M. De.; Stossi, F.; Carlson, K. A.; Katzenellenbogen, B. S.; Katzenellenbogen, J. A. J. Med. Chem. 2005, 48, 1132-1144. (21) Vu, A. T.; Campbell, A. N.; Harris, H. A.; Unwalla, R. J.; Mamas, E. S.; Mewshaw, R. E. Bioorg. Med. Chem. Lett. 2007, 17, 4053-4056. (22) Sun, J.; Meyers, M. J.; Fink, B. E.; Rajendran, R.; Katzenellenbogen, J. A.; Katzenellenbogen, B. S. Endocrinology 1999, 140, 800-804. (23) Shiau, A. K.; Barstad, D.; Radek, J. T.; Meyers, M. J.; Nettles, K. W.; Katzenellenbogen, B. S.; Katzenellenbogen, J. A.; Agard, D. A.; Greene, G. L. Nat. Struct. Biol. 2002, 9, 359-364. (24) Sun, J.; Baudry, J.; Kztzenellenbogen, J. A.; Katzenellenbogen, B. S. Mol. Endocrinol. 2003, 17, 247-258. (25) Shiau, A. K.; Barstad, D.; Loria, P. M.; Cheng, L.; Kushner, P. J.; Agard, D. A.; Greene, G. L. Cell 1998, 95, 927-937.

Zeng et al. (26) Tanenbaum, D. M.; Wang, Y.; Williams, S. P.; Sigler, P. B. Proc. Natl. Acad. Sci. U.S.A. 1998, 95, 5998-6003. (27) Hsieh, R. W.; Rajan, S. S.; Sharma, S. K.; Guo, Y.; Desombre, E. R.; Mrksich, M.; Greene, G. L. J. Biol. Chem. 2006, 281, 17909-17919. (28) Hillisch, A.; Peters, O.; Kosemund, D.; Muller, G.; Walter, A.; Schneider, B.; Reddersen, G.; Elger, W.; Fritzemeier, K. H. Mol. Endocrinol. 2004, 18, 1599-1609. (29) Manas, E. S.; Xu, Z. B.; Unwalla, R. J.; Somers, W. S. Structure 2004, 12, 2197-2207. (30) Frisch, M. J.; Trucks, G. W.; Schlegel, H. B.; Scuseria, G. E.; Robb, M. A.; Cheeseman, J. R.; Montgomery, J. A., Jr.; Vreven, T.; Kudin, K. N.; Burant, J. C.; Millam, J. M.; Iyengar, S. S.; Tomasi, J.; Barone, V.; Mennucci, B.; Cossi, M.; Scalmani, G.; Rega, N.; Petersson, G. A.; Nakatsuji, H.; Hada, M.; Ehara, M.; Toyota, K.; Fukuda, R.; Hasegawa, J.; Ishida, M.; Nakajima, T.; Honda, Y.; Kitao, O.; Nakai, H.; Klene, M.; Li, X.; Knox, J. E.; Hratchian, H. P.; Cross, J. B.; Bakken, V.; Adamo, C.; Jaramillo, J.; Gomperts, R.; Stratmann, R. E.; Yazyev, O.; Austin, A. J.; Cammi, R.; Pomelli, C.; Ochterski, J. W.; Ayala, P. Y.; Morokuma, K.; Voth, G. A.; Salvador, P.; Dannenberg, J. J.; Zakrzewski, V. G.; Dapprich, S.; Daniels, A. D.; Strain, M. C.; Farkas, O.; Malick, D. K.; Rabuck, A. D.; Raghavachari, K.; Foresman, J. B.; Ortiz, J. V.; Cui, Q.; Baboul, A. G.; Clifford, S.; Cioslowski, J.; Stefanov, B. B.; Liu, G.; Liashenko, A.; Piskorz, P.; Komaromi, I.; Martin, R. L.; Fox, D. J.; Keith, T.; Al-Laham, M. A.; Peng, C. Y.; Nanayakkara, A.; Challacombe, M.; Gill, P. M. W.; Johnson, B.; Chen, W.; Wong, M. W.; Gonzalez, C.; Pople, J. A. Gaussian 03, revision C.02; Gaussian, Inc.: Wallingford, CT, 2004. (31) Cieplak, P.; Cornell, W. D.; Bayly, C.; Kollman, P. A. J. Comput. Chem. 1995, 16, 1357-1377. (32) Wang, J.; Wang, W.; Kollman, P. A.; Case, D. A. J. Mol. Graph. Model. 2006, 25, 247-260. (33) Sybyl, Version 7.0; Tripos Inc.: St. Louis, MO, 2004. (34) Case, D. A.; Darden, T. A.; Cheatham, T. E., III; Simmerling, C. L.; Wang, J.; Duke, R. E.; Luo, R.; Merz, K. M.; Pearlman, D. A.; Crowley, M.; Walker, R. C.; Zhang, W.; Wang, B.; Hayik, S.; Roitberg, A.; Seabra, G.; Wong, K. F.; Paesani, F.; Wu, X.; Brozell, S.; Tsui, V.; Gohlke, H.; Yang, L.; Tan, C.; Mongan, J.; Hornak, V.; Cui, G.; Beroza, P.; Mathews, D. H.; Schafmeister, C.; Ross, W. S.; Kollman, P. A. AMBER 9; University of California: San Francisco, 2006. (35) Wang, J.; Cieplak, P.; Kollman, P. A. J. Comput. Chem. 2000, 21, 1049-1074. (36) Wang, J.; Wolf, R. M.; Caldwell, J. W.; Kollman, P. A.; Case, D. A. J. Comput. Chem. 2004, 25, 1157-1174. (37) Jorgensen, W. L.; Chandrasekhar, J.; Madures, J.; Impey, R. W.; Klein, M. L. J. Chem. Phys. 1983, 79, 926-935. (38) Ryckaert, J. P.; Ciccotti, G.; Berendsen, J. C. J. Chem. Phys. 1977, 23, 327-341. (39) Darden, T.; York, D.; Pedersen, L. J. Chem. Phys. 1993, 98, 1008910092. (40) Massova, I.; Kollman, P. A. Perspect. Drug DiscoVery Des. 2000, 18, 113-135. (41) Wang, J.; Morin, P.; Wang, W.; Kollman, P. A. J. Am. Chem. Soc. 2001, 123, 5221-5230. (42) Fogolari, F.; Brigo, A.; Molinari, H. Biophys. J. 2003, 85, 159166. (43) Kollman, P. A.; Massova, I.; Reyes, C.; Kuhn, B.; Huo, S.; Lee, M.; Lee, T.; Duan, Y.; Wang, W.; Donini, O.; Cieplak, P.; Srinivasan, J.; Case, D. A.; Cheatham, . T. E. Acc. Chem. Res. 2000, 33, 889-897. (44) Gilson, M. K.; Sharp, K. A.; Honig, B. H. J. Comput. Chem. 1987, 9, 327-335. (45) Tsui, V.; Case, D. A. Biopolymers 2001, 56, 275-291. (46) Still, W. C.; Tempczyk, A.; Hawley, R. C.; Hendrickson, T. J. Am. Chem. Soc. 1990, 112, 6127-6129. (47) Weiser, J.; Shenkin, P. S.; Still, W. C. J. Comput. Chem. 1999, 20, 217-230. (48) Cheatham, T. E., III; Srinivasan, J.; Case, D. A.; Kollman, P. A. J. Biomol. Struct. Dyn. 1998, 16, 265-280. (49) Kuhn, B.; Kollman, P. A. J. Med. Chem. 2000a, 43, 3786-3791. (50) Pal, D.; Chakrabarti, P. J. Biomol. Struct. Dyn. 2001, 19, 115128. (51) Pranata, J. Bioorg. Chem. 1997, 25, 213-219. (52) Viguera, A. R.; Serrano, L. Biochemistry 1995, 34, 8771-8779.

Insights into Ligand Selectivity in Estrogen Receptor Isoforms

Recommend Documents