Computational Insights into the Mechanism of Ligand Unbinding and

Jul 7, 2009 - Discovery of Potent Ligands for Estrogen Receptor β by Structure-Based Virtual Screening. Jie Shen , Chengfang Tan , Yanyan Zhang , Xi ...
1 downloads 0 Views 3MB Size
Download PDF
10436

J. Phys. Chem. B 2009, 113, 10436–10444

Computational Insights into the Mechanism of Ligand Unbinding and Selectivity of Estrogen Receptors Jie Shen,† Weihua Li,† Guixia Liu,† Yun Tang,*,† and Hualiang Jiang†,‡ School of Pharmacy, East China UniVersity of Science and Technology, 130 Meilong Road, Shanghai 200237, China; and Drug DiscoVery and Design Center, Shanghai Institute of Materia Medica, Chinese Academy of Sciences, 555 Zuchongzhi Road, Shanghai 201203, China ReceiVed: April 24, 2009; ReVised Manuscript ReceiVed: June 9, 2009

Estrogen receptors (ER) belong to the nuclear receptor superfamily, and two subtypes, ERR and ERβ, have been identified to date. The differentiated functions and receptor expressions of ERR and ERβ made it attracted to discover subtype-specified ligands with high selectivity. However, these two subtypes are highly homologous and only two residues differ in the ligand binding pocket. Therefore, the mechanism of ligand selectivity has become an important issue in searching selective ligands of ER subtypes. In this study, steered molecular dynamics simulations were carried out to investigate the unbinding pathways of two selective ERβ ligands from the binding pocket of both ERR and ERβ, which demonstrated that the pathway between the H11 helix and the H7∼H8 loop was the most probable for ligand escaping. Then potentials of mean force for ligands unbinding along this pathway were calculated in order to gain insights into the molecular basis for energetics of ligand unbinding and find clues of ligand selectivity. The results indicated that His524/475 in ERR/ERβ acted as a “gatekeeper” during the ligand unbinding. Especially, the H7∼H8 loop of ERβ acted as a polar “transmitter” that controlled the ligand unbinding from the binding site and contributed to the ligand selectivity. Finally, the mechanism of ligand selectivity of ER subtypes was discussed from a kinetic perspective and suggestions for improving the ligand selectivity of ERβ were also presented. These findings could be helpful for rational design of highly selective ERβ ligands. Introduction The estrogen receptors (ERs) are ligand-activated transcription factors belonging to a superfamily named nuclear receptors (NRs), which play crucial roles in the human body, including the growth and differentiation of the reproductive system, the central nervous system, and the skeletal system.1-3 Like other members of the NR superfamily, the ER is composed of three domains, including an N-terminal domain, a DNA-binding domain (DBD), and a C-terminal ligand binding domain (LBD). Ligands such as endogenous 17β-estradiol (E2) exert their regulatory effects via binding to the LBDs of ERs, which induces conformational changes in the protein structure and produces the physiological effects of the receptor.4 To date, two highly homologous subtypes of ER were identified, namely ERR and ERβ. They share 98% amino acid sequence identity in their DBD and 56% sequence identity in their LBD. Studies of knockout mice have shown that the two subtypes have distinct functions and are differentially expressed in certain tissues.1 ERR is predominantly expressed in the breast, uterus, and vagina, whereas ERβ is widely expressed in several tissues, including the central nervous system, the cardiovascular system, the gastrointestinal system, and the immune system.5,6 The differentiated functions and receptor expressions not only suggest that tissues could be differentially targeted with subtypeselective ligands but also stimulate the search for subtypeselective ligands that can elicit tissue- or cell-specific ER activity. Previous studies suggested that ERβ selective ligands * To whom correspondence should be addressed. E-mail: ytang234@ ecust.edu.cn. † East China University of Science and Technology. ‡ Chinese Academy of Sciences.

may offer some of the benefits of hormone replacement therapy and could be potential drugs for some diseases such as chronic intestinal and joint inflammation without increasing the risk of breast or uterine cancer.7-11 Therefore, discovering ERβ selective ligands has attracted much attention in recent years. At present, however, the mechanism underlying the ligand selectivity of ERs remains unclear, which seriously hampers the discovery and design of selective ER ligands. Crystal structures of the LBD of ERR and ERβ reveal that only two residues differ in the binding pockets: ERR Leu384 and Met421 are substituted by Met336 and Ile373 in ERβ, respectively. Therefore, many previous studies tried to interpret the possible mechanism of ligand selectivity based on these two pairs of residues. It has been suggested that ERR Met421 was the major contribution to ERR selectivity,12-14 whereas ERβ selectivity could be ascribed to Met336.15-18 However, in addition to these two residues, others may also contribute to the ligand selectivity. For example, the studies on ER subtype chimeras generated by DNA shuffling showed that some residues in the third helix (H3) of the ERβ LBD determined the ERβ selectivity of diarylpropionitrile (DPN),15 and some residues from the eighth helix (H8) to the middle part of the eleventh helix (H11) were required for ERR selectivity of propylpyrazole triol (PPT).19 A recent study of our group indicated that hydrophobic residues such as Leu346/298 and Phe404/356 of ERR/ERβ also contributed to the ligand selectivity.20 Previous studies also showed that the substrate specificity may be influenced not only by the interaction in the binding site but also by the binding or unbinding processes.21-23 Therefore, such dynamic processes might also contribute to the ligand selectivity. Nevertheless, currently little information is available on ligands binding to or unbinding from NRs, since the ligand-binding

10.1021/jp903785h CCC: $40.75  2009 American Chemical Society Published on Web 07/07/2009

Ligand Unbinding and Selectivity of Estrogen Receptors

Figure 1. Structures of ER ligands: (a) Way-244 and (b) Genistein.

pocket of NR is buried in the hydrophobic core of the proteins with no obvious entry or exit route. Lots of previous works mainly focused on the ligand unbinding processes of different NRs using different MD protocols.24-31 For example, a study of retinoic acid receptors (RARs) using enhanced sampling molecular dynamics (ESMD) suggested that ligand unbinding was connected with detachment of H12, hence confirming the “mousetrap” mechanism.24 Kosztin and colleagues investigated three possible pathways of retinoic acids unbinding from RARs using steered molecular dynamics simulations (SMD) and suggested one pathway for ligand binding and a “back door” pathway for ligand unbinding.25 In addition, the unbinding pathways of RAR have also been explored by random expulsion molecular dynamics (REMD), with the result that displacing H12 was not necessary for ligand unbinding.27 Polikarpov et al. performed different molecular dynamics simulations including locally enhanced sampling (LES),26 SMD,28 and multipoint SMD32 to explore the ligand binding and unbinding of the thyroid hormone receptor’s (TR) LBD. Three possible pathways were observed from the LES study,26 and then details of the three pathways were investigated by SMD, with the conclusion that ligand dissociation was not favored through H12.28 Furthermore, the ligand entry into the LBD of TR was also studied using multipoint SMD (MP-SMD).32 Recently, Sonoda and Polikarpov et al. suggested that different ER ligands could escape from the human ERR (hERR) LBD along four possible unbinding pathways in the LES MD simulations.30 Two of the four pathways involved repositioning of H12 in accordance with the “mousetrap” mechanism,33 and the two others involved the separation of H8 and H11, which were also found in other NRs.27,28 However, LES failed to identify which pathway was the most feasible for ligand unbinding due to lack of quantitative information. In this work, SMD simulations were performed to investigate the dynamic processes of two ERβ selective ligands (WAY24416 and Genistein5) unbinding from both ERR and ERβ, which aimed to explore the most possible unbinding pathway as well

J. Phys. Chem. B, Vol. 113, No. 30, 2009 10437 as its influences on ligand selectivity. WAY-244 is a derivative of benzofuran with 75.8-fold selectivity to ERβ over ERR, and Genistein is a plant-derived nonsteroidal compound with the ERβ selectivity of 8-40-fold (Figure 1). In SMD, an external force was applied to the bound ligand in order to facilitate the ligand’s unbinding from the receptors along four different predetermined pathways identified by LES MD simulations separately (Figure 2).30 To further understand the energetics of ligand unbinding and its contribution to ligand selectivity, potentials of mean force (PMFs) were calculated from a number of repeated SMD trajectories based on the Jarzynski’s equality.34 The PMF results of ligand unbinding correlated well with the ligand selectivity. Furthermore, the differences between ligand unbinding from ERR and ERβ were also found to have contributions to the ligand selectivity. Methods and Models System Preparation. The X-ray crystal structures of the ERR and ERβ bound with Way-244 extracted from the Protein Data Bank (PDB entries: 1X7E and 1X78)16 were used as initial structures for models S1 and S2. And the structures of ERR and ERβ bound with Genistein (PDB entries: 1X7R and 1X7J)17 were used as initial structures for models S3 and S4 (Table 1). Missing side chains of residues were modeled with Sybyl 7.0 (Tripos, Inc., St. Louis, MO), and missing residues were modeled with Modeller 9.235 using 1QKU and 1QKM as templates. The protonation states of histidines were determined using Propka 2.0.36 The ligands were optimized at the Hartree-Fock level with the 6-31G(d,p) basis set using Gaussian 03 (Gaussian Inc., Pittsburgh, PA), and then restrained electrostatic potential (RESP) charges were calculated using the B3LYP/cc-pVTZ quantum mechanical method. The Amber ff03 force field37 and general Amber force field (GAFF)38 were applied to proteins and ligands, respectively. Topology and parameter files were generated using the LEaP program, and ligand force field parameters of ligands were developed by the Antechamber module39 in the program AmberTools1.2. Molecular Dynamics Setup and Equilibrium. All MD and SMD simulations were carried out using the Amber10 package40 with the same protocols. Each simulation system was surrounded by a truncated octahedron periodic box of TIP3P water molecules41 with a margin of 10.0 Å along each dimension, and sodium ions were added to maintain charge neutrality. All covalent bonds to hydrogen atoms were constrained using the SHAKE algorithm.42 Electrostatic interactions were calculated

Figure 2. Four different pathways for ligand unbinding. (The initial position of the ligand is colored in gray.)

10438

J. Phys. Chem. B, Vol. 113, No. 30, 2009

Shen et al.

TABLE 1: Models for Both ERr and ERβ Complexes with Different ERβ Selective Ligands CR of residue in protein f pulled atom in ligand model

ligand

protein

PDB entry

resolution (Å)

path 1

path 2

S1 S2 S3 S4

WAY-244 WAY-244 Genistein Genistein

hERR hERβ hERR hERβ

1X7E 1X78 1X7R 1X7J

2.80 2.30 2.00 2.30

I424 f C5 I376 f C5 M427 f C10 M379 f C10

R412 f C1 F377 f C1 F425 f C2 F377 f C2

using the particle-mesh Ewald (PME) algorithm43 with a cutoff of 10 Å applied to Lennard-Jones interactions. Periodic boundary conditions were applied to avoid edge effects in all calculations. Before the MD production, 500 steps steepest-descent minimization and 1500 steps conjugated gradient minimization were applied to the solvent and the entire model, respectively. Then the system was heated up from 0 to 300 K gradually over 30 ps using the NVT ensemble with the solutes restrained by a weak harmonic potential. Afterward, equilibrations totaling 140 ps were carried out in the NPT ensemble via three steps: the solutes were restrained while the waters and counterions were equilibrated in the first 20 ps; then the side chains of the protein were relaxed in the next 20 ps; and all of the restraints were removed in the last 100 ps. Finally, a 5-ns MD simulation was conducted at 1 atm and 300 K under the NPT ensemble with a time step of 2 fs. Temperature was controlled using Langevin dynamics. The final equilibrium states of the MD productions were saved as starting points for the followed SMD studies. Coordinates were saved every 1 ps during the entire process. Steered Molecular Dynamics Simulation. SMD simulations were performed under the NPT condition with Langevin dynamics, and external restraints were applied to pull the ligand out of its initial binding site. For each system, the initial pulling directions were defined according to pathways obtained from the previous study (Figure 2).30 For each pathway, the pulling direction was determined with a vector connecting a specified CR atom in the protein with the pulled atom in the ligand (Table 1). In the steering processes, the specified CR atom and three additional CR atoms of Leu453, Ile482, and Arg503 in ERR or Lys401, Val433, and Leu458 in ERβ were restrained in order to avoid protein translation and rotation. A constant pulling velocity of 0.02 Å/ps was used, which was considerably slow compared to the cases of other SMD studies. According to the stiff spring approximation,44,45 a sufficiently stiff spring constant of 10 kcal · mol-1 · A-2 was used. The pulling force F at time t in the SMD was calculated as follows:

F(t) ) 2k(Vt - x(t))

(1)

where k is the spring constant of the constraint, V is the pulling velocity, and x(t) is the position of the ligand at time t. The work W was calculated by integrating the force over the pulled trajectory:

W(x(t)) )

∫0x(t) F(t) dx(t)

(2)

The trajectories were saved for every 1 ps, and the calculated force and work were saved for every step. Each trajectory along the individual pathway was repeated 10 times with different pseudorandom number seeds. For the specified pathway which was used to calculate the PMF, the trajectories were subsequently repeated for more times in order to reduce the statistical error due to insufficient sampling in the PMF calculation.

path 3 L387 L339 L387 L339

f f f f

C1 C1 C3 C3

path 4 E380 f C6 S338 f C6 E380 f C1 E332 f C1

Potentials of Mean Force Construction. The Jarzynski’s equality established a connection between the nonequilibrium work and equilibrium free energies, which has been widely used to estimate the PMF or nonequilibrium free energy difference from the work.46-51 The formula of Jarzynski’s equality is

exp(-β∆F) ) 〈exp(-βW0fλ)〉ave

(3)

where W0fλ is the work performed from state 0 to state λ, β ) 1/(κBT), κB is Boltzmann constant, and T is the temperature. In order to use this equation directly, a number of trajectories with sufficiently small W0fλ values are needed to confirm the accuracy. However, such trajectories are rarely sampled in practice, and therefore, some approximate formulas and methods were introduced to deal with such problems from the perspective of systematic error and statistical error.44,45,52 In this work, the cumulative integral (CI) extrapolation method developed by Ytreberg and Zuckerman was used,52 which uses an integral block-averaged work calculated from the results of the Jarzynski’s equality (eq 3). It has been proved to be accurate in reducing the required data by 5-40-fold and works well in practical cases.49,52 Here, a 25-ps (0.5-Å) sampling window was used in order to obtain converged PMF data. The scripts of the CI extrapolation method were downloaded from the Web site of Zuckerman’s lab: http://www.ccbb.pitt.edu/Faculty/zuckerman/ software.html. Results Prior to the SMD simulations, all four systems were subjected to the MD equilibrium. All the systems have reached equilibrium after 5 ns MD simulations (Figure S1, Supporting Information). The final snapshots at 5 ns MD simulations were saved as starting points for subsequent SMD simulations. Details of Ligand Unbinding Pathways. The ligands were pulled out of the binding site along four pathways, as shown in Figure 2. Path 1 is located between the H11∼H12 loop and the N-terminal part of the H3 helix, which corresponds to the mousetrap mechanism. Path 2 involves the displacement of the H12 helix, where the ligand escaped above the H11∼H12 loop. Path 3 is located between the H7∼H8 loop and the H11 helix. In path 4, ligands escaped between the H7 helix and the C-terminal part of the H2∼H3 loop, which caused the unfolding of H7. Here, the restraint on the ligand was spherical, and only the distance between the atoms in the protein and in the ligand was changed. Therefore, ligand would find a way to escape according to the initial predefined direction. Ten trajectories were repeated for each pathway. Hence, a total of 40 SMD simulations were performed on each system. The trajectories of ligand escaping from the receptor successfully along the initial direction within 1 ns were accepted. The average work and the number of accepted trajectories of the ligands unbinding along each pathway were shown in Table 2.

Ligand Unbinding and Selectivity of Estrogen Receptors

J. Phys. Chem. B, Vol. 113, No. 30, 2009 10439

TABLE 2: Average Work of Ligand Unbinding for Each Pathway path 1

path 2

modela

n

〈W〉

n

S1 S2 S3 S4

8 0 4 0

80.84

0 0 8 7

102.83

path 3

path 4

〈W〉

n

〈W〉

n

〈W〉

111.50 78.09

10 10 10 8

65.25 85.18 66.93 72.01

0 0 6 3

77.54 85.88

a A total of 40 SMD trajectories were carried out in each model. n is the number of trajectories with ligand unbinding along the pathway successfully; 〈W〉 is the average work for the n trajectories. b

From Table 2, it is easy to see that path 1 was only observed for the ERR system (model S1 and S3); only Genistein could unbinding along paths 2 and 4 (models S3 and S4); and path 3 was observed for all the four models. Furthermore, ligands escaped from either ERR or ERβ along path 3 with the least average work and the highest probability among four models. Ligands Unbinding along Path 1. Path 1 was only observed for ligand successfully unbinding from ERR. This pathway is lined by the residues in the H11∼H12 loop and the H3 helix. During the ligand unbinding, the N-terminal part of H3 and the C-terminal part of H11 first moved apart and some of the contacts between them were broken. Then, the R-helix from Glu339 to Leu349 in the N-terminus of H3 bent outward, which caused a deformation between Ler349 and Arg352 and thus formed a cavity. After the ligand successfully escaped, the N-terminal part of H3 moved back. At the same time, the C-terminus of H11 together with the H11∼H12 loop moved back to their initial positions. In the case of Way-244 unbinding from ERR (model S1), 8 of 10 trajectories were repeated successfully. During the ligand unbinding, Way-244 pushed the N-terminal region of H3 and thus moved it outward from 200 to 600 ps. The first force peak appearing from 250 to 400 ps was mainly due to the disruption of the initial interactions between Way-244 and the residues in the binding site, including Glu353 and Arg394, which hydrogenbonded with the phenolic hydroxyl of Way-244 and a highly ordered water molecule. The hydrogen bond between the water and Arg394 broke after 400 ps. The second force peak appearing in 600 ps was mainly due to the breakage of the hydrogen bond between the phenolic hydroxyl and Glu353; meanwhile, the side chain of Glu353 moved back and the hydrogen bond with the water was also broken (Figure 3a, black line). Genistein bound with ERR more tightly than Way-244; thus, more work was required to pull it out. Only 4 of 10 trajectories were repeated successfully. The first force peak appearing at ∼200 ps was mainly caused by the breakage of the hydrogen bond between the Genistein 2-OH group and the Nδ atom of His524. The hydrogen-bonding network involved with the phenolic hydroxyl of Genistein was not broken until 400 ps. A water molecule mediated hydrogen-bonding network was formed by the phenolic hydroxyl and the residues His524 and Gly521. As Genistein moved out, increased force was required to break them. This corresponded to the second force peak from 500 to 700 ps (Figure 3a, blue line). Ligands Unbinding along Path 2. Path 2 is neighboring with path 1, which was above the H11∼H12 loop. It was only observed for Genistein escaping from both ERR and ERβ. For the ERR-Genistein complex (model S3), Genistein rotated clockwise ∼90° (from top view) on its prolate axis from 200 to 300 ps before unbinding. Meanwhile, the hydrogen bond between the Genistein 2-OH group and the Nδ atom of His524 was broken. The hydrogen-bonding network involved in the

Genistein phenolic hydroxyl group began to break after 400 ps. The force peak at 500 ps was due to the breakage of the hydrogen bond between the phenolic hydroxyl group and Glu353 (Figure 3b, blue line). In contrast to the case for path 1, the water molecule in the binding site did not leave the pocket along with the ligand. For the ERβ-Genistein complex (model S4), Genistein also contrarotated ∼90° on its prolate axis from 200 to 400 ps before ligand unbinding. Meanwhile, the major interactions in the binding pocket were corrupted, corresponding to the force valley at ∼300 ps. As Genistein moved out of the binding pocket, the polar functional groups on the dual ring part formed a strong interaction with the backbone atoms of the H11∼H12 loop after 500 ps; at that time, a new hydrogen-bonding network formed, involving the residues Val485 and Asn483. These interactions became a major hindrance, which prevented the ligand from escaping. For this reason, the force increased and remained on a high level until the ligand left the binding site (Figure 3b, green line). Ligands Unbinding along Path 3. Path 3 was observed for all four models. The length of this pathway was the shortest among the four pathways, and the ligands could escape out of the receptors within 800 ps. In addition, ligands could successfully escape along this pathway in almost all the repeated trajectories and the average work was much lower than those of other pathways. In the case of Way-244 unbinding from the ERR (model S1), there was little hindrance along this pathway. It was different from other models that the C-terminal part of H11 moved outward and the pathway was opened before the ligand unbinding. This phenomenon was also observed during the MD simulation and leading to a high fluctuation in this part, which could be mostly ascribed to the lower binding affinity of Way244 with ERR.20 With the ligand being pulled out, the force increased slowly to break the original interactions between Way244 and the residues in the binding site of ERR. The hydrogenbinding network formed by the phenolic hydroxyl group and the residues between Glu353 and Arg394 via a water molecule began to corrupt at ∼200 ps, which led to a decrease in the pulling force. After Way-244 left the binding pocket, the ligand pushed outward the C-terminal part of the H7∼H8 loop and the side chain of His524, expanding the cavity for the ligand out (Figure 3c, black line). Compared to ERR, Way-244 bound with ERβ more tightly and had 75.8-fold selectivity over ERR. Therefore, Way-244 was totally confined in the binding pocket of ERβ (model S2). The force increased dramatically in the first 200 ps to counterbalance the strong interactions between Way-244 and the residues in the binding pocket in the process of the ligand unbinding (Figure 3c, red line). The initial interactions began to lose at ∼200 ps and totally disappeared after 300 ps, and therefore, the force dropped quickly. In addition, the steric clash between the acetonitrile group of Way-244 and the residues Ile373 and Ile376 caused the movement of the H7∼H8 loop and breakage of the hydrogen bond between the Nε atom of His475 and the carbonyl group of Glu371. The dihedral angle between acetonitrile and the benzofuran plane flipped in some trajectories to reduce the clash. Then a hydrogen bond was formed between the benzofuran hydroxyl of Way-244 and Glu371. With the ligand moving out, another weak hydrogen bond formed between the acetonitrile group and Gly372 at ∼500 ps. The force peak at ∼600 ps was mainly due to the breakage of these two hydrogen bonds.

10440

J. Phys. Chem. B, Vol. 113, No. 30, 2009

Shen et al.

Figure 3. Average forces of different pathways for each model. (Black line for model S1, where Way-244 is unbinding from ERR; red line for model S2, where Way-244 is unbinding from ERβ; blue line for model S3, where Genistein is unbinding from ERR; and green line for model S4, where Genistein is unbinding from ERβ. The same colors for models S1-S4 were used in Figures 4, 5, and 7.)

Genistein bound with both ERR and ERβ with higher binding affinity. Hence, the protein conformations for both ERR and ERβ maintained firm and the ligands were buried in the binding pockets during the MD equilibrium. In the case of Genistein unbinding from ERR (model S3), the initial force peak at 0-200 ps was primarily due to the breakage of the hydrogen-bonding network in the binding site (Figure 3c, blue line). Meanwhile, the His524 side chain was pushed away by Genistein. The H7∼H8 loop and the N-terminal region of H8 also moved outward, forming a cavity for ligand escaping. After the ligand moved out of the binding pocket, a hydrogen bond between the carboxyl group of Glu419 and the 2-OH group of Genistein formed. In addition, a water molecule mediated hydrogen bond formed between the phenol hydroxyl group of Genistein and the carboxyl group of Leu346. Therefore, a relatively large force was required to break these interactions, which resulted in a large force peak from 400 to 500 ps. For the ERβ-Genistein complex (Model S4), the force profile was almost similar to that of the ERR-Genistein complex. The difference was that the force was maintained at a high level for a little longer time in the first 300 ps (Figure 3c, green line). This was mainly due to the stronger interactions of Genistein with ERβ than with ERR. The initial interactions between Genistein and the residues in the binding pocket began to decrease after 200 ps. With Genistein pulled out, a relatively large force was needed to compensate the clash with the residues in the H7∼H8 loop and the breakage of the hydrogen bond between the Nε atom of His475 and the carbonyl group of Glu371, which corresponds to the minor force peak from 250 to 300 ps. Subsequently, the side chain of His475 was pushed outward, and the C terminal part of H11 and the H7∼H8 loop departed, forming a cavity for ligand escaping. Because of its high polarity, the residues in the H7∼H8 loop hydrogen-bonded with the hydroxyl groups or the carbonyl group of Genistein, leading the force to drop at ∼400 ps. Then the force increased to counterbalance the hydrogen-bonding interactions from 550

to 700 ps. The ligand completely escaped out of the receptor after 800 ps. Ligand Unbinding along Path 4. Only the Genistein could escape from the receptors along this pathway; moreover, such processes were rarely sampled. For ERR (model S3), Genistein changed its orientation in the binding pocket in the first 350 ps with the force increased. After that, Genistein began to leave the binding pocket with the breakage of the hydrogen bonding network. Therefore, the force dropped quickly at ∼350 ps. As the ligand moving out, H7 was partially unwinding and moved outward. The force remained high mainly due to the strong polar contacts between the polar function groups on Genistein and the residues in H7 during ligand unbinding. The decreased force after 500 ps was mainly due to the reduction of steric clash and the contacts switching from 2-OH to 4-OH and 6-carbonyl (Figure 3d, blue line). In ERβ (model S4), the orientation of Genistein changed in the first 300 ps as well as ERR. The initial interactions between Genistein and ERβ were corrupted after 400 ps, corresponding to the dropping in force (Figure 3d, green line). As Genistein moved out, the N-terminal parts of H3 and H7 were unwinding and separated due to the steric clash. Therefore, a large force was required to pull Genistein out. In light of the results of force profiles (Figure 3), the average work, and the escaping probability (Table 2) of the different pathways, we assumed that path 3 is the most probable pathway for Way-244 and Genistein unbinding from both ERR and ERβ. Therefore, we only calculated the PMFs of the trajectories for ligands unbinding from path 3. Construction of the Potential of Mean Force along Path 3. To obtain reasonable PMF profiles of the nonequilibrium process of ligand unbinding from the ER, a number of external works have been sampled by repeating the SMD trajectories for each model many times. In total, 33 trajectories for model S1, 45 trajectories for model S2, 45 trajectories for model S3, and 45 trajectories for model S4 were performed, and the PMF

Ligand Unbinding and Selectivity of Estrogen Receptors

Figure 4. Histograms of work distributions for different ligands unbinding from both ERR and ERβ.

J. Phys. Chem. B, Vol. 113, No. 30, 2009 10441 size of Genistein was smaller than that of Way-244, and thereby, there was less steric clash, corresponding to a lower PMF. It should be noted that the PMF decreased a little after Genistein pulled out of the binding site. We found that besides the carboxyl of Glu371 hydrogen-bonded with the 2-OH group of Genistein, another hydrogen bond formed between the carbonyl group of Gly372 and the 4-OH group of the Genistein. These strong interactions contributed to the decrease of free energy of the system corresponding to the PMF trough from 12 Å to 15 Å. The PMF of Genistein unbinding from ERβ was converged to 56.9 kcal/mol, which was 2.3 kcal/mol higher than that for the case of ERR. The difference in the PMF of Way-244 unbinding from ERR and ERβ was 30.9 kcal/mol, which was significantly higher than that of Genistein. It was consistent with the different folds of selectivity of Way-244 (75.8-fold) and Genistein (8-40-fold). The results also confirmed that path 3 was the reasonable unbinding pathway for Way-244 and Genistein escaping from both ERR and ERβ. Discussion

Figure 5. PMFs of ligands unbinding along path 3.

profiles of ligands unbinding from both ERR and ERβ were calculated on the basis of Jarzynski’s equality using the CI extrapolation method. The work distribution for each model approximately satisfied the Gaussian distribution (Figure 4). The width of the work distribution was in a proper size as shown in Figure 4. Figure 5 shows the results of PMF profiles of the ligand unbinding reaction coordinate (RC) along path 3. For model S1, the PMF of Way-244 unbinding increased continuously in the first 5 Å of pulling (Figure 5, black line). At this point, the ligand began to leave the binding site with the breakup of the initial interactions and forming new interactions with the residues outside of the binding site. The PMF maintained a plateau from 5 Å to 7 Å, mainly due to the counterbalance of interactions switched. The PMF increased as the Way-244 pulled out and finally converged to 41.7 kcal/mol. In model S2, Way-244 bound tightly with ERβ and had 75.8fold selectivity over ERR. The PMF shows a steep increase in the first 6 Å of pulling, corresponding to the much more strong interactions of Way-244 with the residues in the binding site of ERβ (Figure 5, red line). Finally, the PMF of Way-244 unbinding from ERβ was 72.6 kcal/mol, which was much higher than that for the case of ERR. The PMF profile of Genistein unbinding from the binding site of ERR was similar with Way-244 in the first 10 Å of pulling (Figure 5, blue line). Because Genistein bound with ERR was more stable than Way-244, the unbinding pathway of the Genistein-ERR complex was buried unless the ligand was pulled out. As a result, the PMF increase after 10 Å was mainly due to the clashes between Genistein and the residues out of the binding pocket. The final PMF of Genistein unbinding from ERR was 54.6 kcal/mol. For model S4, the PMF plot was almost overlapped with that of model S2 in the first 3 Å (Figure 5, green line). For this reason, the behaviors of the ligands were almost the same in the initial steps of unbinding. However, the

Path 3 is the Most Probable Ligand Unbinding Pathway. In this study, two ERβ selective ligands, Way-244 and Genistein, were forced to leave the binding sites of both ERR and ERβ along four predefined pathways. The processes were repeated for 10 times, and the results were analyzed thoroughly and carefully. The results demonstrated that the average work for pulling the ligands unbinding along path 3 was much lower than other pathways among the four models. And almost all the repeated trajectories for this pathway were accepted (Table 2). The force profiles for pulling ligand out of the receptors along different pathways differed. For each individual model, the magnitude of the force peaks for path 3 was smaller than those of the other three pathways. In addition, this pathway was also observed in thyroid hormone receptors (TRs),27,28 which raised the possibility that path 3 was reasonable for ligands unbinding from ERs. Furthermore, the relative PMF values in different subtypes also suggested that path 3 was reasonable, since the PMFs of ERβ selective ligands unbinding along path 3 from ERβ were larger than those of ERR. Previous mutagenic experiments of ERR also suggested that mutations at Gly521, His524, Leu525, and Met528, which are located around path 3, had the greatest effects on E2-induced transcription.53 Key Role of ERr His524/ERβ His475 in Ligand Unbinding. For the ERβ selective ligands, the binding affinity with ERβ is higher than that with ERR. Therefore, a larger force is required to pull the ligands out of the binding pocket of ERβ. During the ligand unbinding, the residue His524/475 (ERR/ERβ) played important roles. The side chain of His524/475 (ERR/ ERβ) flipped as ligand unbound and opened the pathway for ligand escaping. The bottleneck of Way-244 and Genistein unbinding for both ERR and ERβ was the clashes between the ligand and the side chain of His524/475 (ERR/ERβ), which caused the force to increase and the PMF to increase significantly (Figures 3 and 5). In addition, the side chain of His524/475 formed hydrogenbonding interactions with Glu419/371 or Gly420/372, which made a major contribution to keeping the H11 and H7∼H8 loop closed and covering the ligand unbinding pathway. The mutagenic experiment of H524A showed 35-fold reducing responsiveness on E2-induced transcriptional activity, which has confirmed the importance of His524.53 Molecular dynamics simulation also suggested the role of His524 in maintaining the

10442

J. Phys. Chem. B, Vol. 113, No. 30, 2009

Shen et al.

Figure 6. Movements of the H7∼H8 loop during Genistein escaping from ERβ: green, RC ) 0 Å; yellow, RC ) 5.1 Å; magenta, RC ) 8 Å; orange, RC ) 11.9 Å; cyan, RC ) 15.5 Å.

protein structure in the biologically active agonist conformation.54 Furthermore, crystal structures of ERR reveal that the side chain of His524 in the apo-ERR Y537S flipped and hydrogen-bond with the backbone of H11.55 Consequently, His524/His475 performs like a “gatekeeper” covering the ligand unbinding pathway. Contributions of the H7∼H8 Loop to Ligand Selectivity on ERβ. In the cases of ligand unbinding from ERβ, the polar residues in the C-terminal part of the H7∼H8 loop including Glu371 and Gly372 also played an important role. Figure 6 shows mobility of the H7∼H8 loop during Genistein unbinding for ERβ. This region acted as a polar “transmitter” to control the ligand escaping. It moved outward in the first 8 Å of pulling as Genistein left the binding site. Then hydrogen bonds were formed between the polar functional groups on the ligands and the backbone atoms or the polar side chains of the H7∼H8 loop. Therefore, Genistein remained close to the C-terminal part of H7∼H8 loop from 12 Å to 16.2 Å. After that, Genistein separated from the H7∼H8 loop “transmitter” and escaped out of the ERβ with the breakage of these hydrogen bonds. Since the Genistein had more polar function groups in the dual-ring part than Way-244, the interactions between Genistein and the H7∼H8 loop were counteracted more than that of Way-244. Therefore, the PMF of model S4 showed a trough from 12 Å to 15 Å (Figure 5, green line). Compared with Genistein, Way244 has fewer polar functional groups and a larger branch; therefore, the interactions between Way-244 and the H7∼H8 loop were not sufficient to counteracted the steric clashes. Therefore, the PMF kept increasing after 10 Å of pulling (Figure 5, red line). Nevertheless, the slope was smaller than those of ERR (Figure 5, black line and blue line) due to such counteractions. Interestingly, ligands kept closer to the H7∼H8 loop in ERβ than they did in ERR during their unbinding. The distances between the ligand and the N-terminal part of H11 were monitored during the process of ligands unbinding as well as the distances between the ligand and the C-terminal part of the H7∼H8 loop, and the average distances were shown in Figure 7. This suggested that the ligand unbinding pathway in ERR was closer to the H11 helix, while the pathway in ERβ was closer to the H7∼H8 loop. This might be caused by two factors. First, Ser527 in H11 of ERR is substituted with Asn478 in ERβ. The steric clash between the ligand and ERβ Asn478 in H11 kept the ligand away from the H11. Yet the side chain of Ser527 in ERR is smaller, and thus the clash was also reduced. On the other hand, the Ile373 in the H8 N-terminus of ERβ is

Figure 7. Average distances between the mass center of Genistein and residues in receptors varied along the escaping pathway: (a) the distance between Genistein and the N-terminal part of H11 (residues 521-531 in ERR and residues 472-482 in ERβ); (b) the distance between Genistein and the C-terminal part of H7∼H8 loop (residues 419-420 in ERR and residues 371-372 in ERβ).

substituted with Met421 in ERR. Quantum chemical calculations proved that the electrostatic repulsive interaction of methionine was stronger than that of isoleucine,16 and therefore, the side chain of ERR Met421 prevented the ligands from approaching the H7∼H8 loop. The different behaviors of ligands unbinding from ERR and ERβ could contribute to the ligand selectivity. Since the ligand unbinding pathway in ERβ was close to the H7∼H8 loop, the transmitting function of the H7∼H8 loop was more significant in ERβ than in ERR. Therefore, the process of polar ligand unbinding from ERβ cannot be easily affected by other factors due to the controlling functions of the H7∼H8 loop. In other words, the correlation of ligand polarity and ligand unbinding from ERβ was more complicated than that from ERR. This was also consistent with the results of the 3D QSAR study, which concluded that more electrostatic field descriptors were needed in the ERβ CoMFA model than that of ERR.14 Implications for Designing Selective Ligands for ERβ. To improve the ligand selectivity for ERβ, we could either increase the ligand binding affinity with ERβ or reduce the binding affinity with ERR. From a kinetic perspective, this would correlate with the decrease in PMF for ligand unbinding from ERR or the increase in PMF for ligand unbinding from ERβ. The problem is how to make an optimum between the two aspects. The ligand unbinding pathway mainly consisted of hydrophobic residues. Therefore, it is possible that increasing the polar

Ligand Unbinding and Selectivity of Estrogen Receptors surface of the ligand would benefit dissociation. We could add polar functional groups such as halogen atoms or hydroxyls on the ligand or change the polarity of the skeleton. However, the polar H7∼H8 loop was close to the unbinding pathway in ERβ, which could form polar contacts with the dual-ring part of ligands. Therefore, ligand escaping from ERβ was controlled by the polar “transmitter” of the H7∼H8 loop. Once the polar functional groups on the dual-ring part of the ligand were decreased, then the polar contacts between the ligand and the “transmitter” were reduced. Meanwhile, the steric clashes were also reduced, which had more impact on ERR and caused the PMF of ligand escaping from ERR to decrease more than that of ligand escaping from ERβ. As a result, the PMF of ligand escaping from ERβ would increase more than that of ligand escaping from ERR. On the other hand, increasing the polarity of ligands by adding the polar functional groups on the phenol part of the ligand would affect more the process of ligand unbinding from ERR, while the ligand unbinding from ERβ was controlled by the “transmitter”. Therefore, the PMF of ligand unbinding from ERR decreased much more than that of ligand unbinding from ERβ. On the basis of these results, we suggested that the ERβ selectivity could be improved by cutting the polar function group on the dual-ring part or adding a polar function group on the phenol part of the ligand. Conclusions The present SMD simulations of two ERβ selective ligands, Way-244 and Genistein, unbinding along four different pathways from both ERR and ERβ, suggested that path 3 was the most probable unbinding pathway for ligand unbinding from both ERR and ERβ. The comparative results of PMF for ligands unbinding also confirmed the conclusion. During the ligand unbinding, His524/His475 in ERR/ERβ served as a “gatekeeper”. Especially the H7∼H8 loop acted as a polar “transmitter”, controlling the ligand unbinding from ERβ, which contributed to the ligand selectivity of ERβ. Finally, the implications for designing ERβ selective ligand were discussed. These findings could be used to interpret the mechanism of ER selectivity and could open the opportunities to find novel ERβ highly selective ligands. Acknowledgment. We thank Dr. Dave Case for his kindness to offer us the Amber10 package. This work was supported by the Program for New Century Excellent Talents in University (Grant No. NCET-08-0774), the 863 High-Tech Project (Grant No. 2006AA020404), the 111 Project (Grant No. B07023), the National Natural Science Foundation of China (Grant No. 30600785), and the Shanghai Rising-Star Program (Grant No. 07QA14016). Supporting Information Available: The rms deviations of the receptors and ligands for the four models in MD simulations are shown in Figure S1, which is available free of charge via the Internet at http://pubs.acs.org. References and Notes (1) Couse, J. F.; Korach, K. S. Endocr. ReV. 1999, 20, 358. (2) Pettersson, K.; Gustafsson, J. A. Annu. ReV. Physiol. 2001, 63, 165. (3) Zhao, C.; Dahlman-Wright, K.; Gustafsson, J. A. Nucl. Recept. Signaling 2008, 6, e003. (4) Kushner, P. J.; Agard, D. A.; Greene, G. L.; Scanlan, T. S.; Shiau, A. K.; Uht, R. M.; Webb, P. J. Steroid Biochem. Mol. Biol. 2000, 74, 311. (5) Kuiper, G. G.; Carlsson, B.; Grandien, K.; Enmark, E.; Haggblad, J.; Nilsson, S.; Gustafsson, J. A. Endocrinology 1997, 138, 863.

J. Phys. Chem. B, Vol. 113, No. 30, 2009 10443 (6) Korach, K. S.; Emmen, J. M.; Walker, V. R.; Hewitt, S. C.; Yates, M.; Hall, J. M.; Swope, D. L.; Harrell, J. C.; Couse, J. F. J. Steroid Biochem. Mol. Biol. 2003, 86, 387. (7) 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, J. C., Jr. Endocrinology 2003, 144, 4241. (8) Malamas, M. S.; Manas, E. S.; McDevitt, R. E.; Gunawan, I.; Xu, Z. 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. (9) Hillisch, A.; Peters, O.; Kosemund, D.; Muller, G.; Walter, A.; Schneider, B.; Reddersen, G.; Elger, W.; Fritzemeier, K. H. Mol. Endocrinol. 2004, 18, 1599. (10) 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. (11) Harris, H. A.; Bruner-Tran, K. L.; Zhang, X.; Osteen, K. G.; Lyttle, C. R. Hum. Reprod. 2005, 20, 936. (12) Bhat, R. A.; Stauffer, B.; Unwalla, R. J.; Xu, Z.; Harris, H. A.; Komm, B. S. J. Steroid Biochem. Mol. Biol. 2004, 88, 17. (13) Tuccinardi, T.; Bertini, S.; Martinelli, A.; Minutolo, F.; Ortore, G.; Placanica, G.; Prota, G.; Rapposelli, S.; Carlson, K. E.; Katzenellenbogen, J. A.; Macchia, M. J. Med. Chem. 2006, 49, 5001. (14) Salum, L. B.; Polikarpov, I.; Andricopulo, A. D. J. Chem. Inf. Model. 2008, 48, 2243. (15) Sun, J.; Baudry, J.; Katzenellenbogen, J. A.; Katzenellenbogen, B. S. Mol. Endocrinol. 2003, 17, 247. (16) Manas, E. S.; Unwalla, R. J.; Xu, Z. B.; Malamas, M. S.; Miller, C. P.; Harris, H. A.; Hsiao, C.; 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. (17) Manas, E. S.; Xu, Z. B.; Unwalla, R. J.; Somers, W. S. Structure 2004, 12, 2197. (18) 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. (19) Nettles, K. W.; Sun, J.; Radek, J. T.; Sheng, S.; Rodriguez, A. L.; Katzenellenbogen, J. A.; Katzenellenbogen, B. S.; Greene, G. L. Mol. Cell 2004, 13, 317. (20) Zeng, J.; Li, W.; Zhao, Y.; Liu, G.; Tang, Y.; Jiang, H. J. Phys. Chem. B 2008, 112, 2719. (21) Winn, P. J.; Ludemann, S. K.; Gauges, R.; Lounnas, V.; Wade, R. C. Proc. Natl. Acad. Sci. U.S.A. 2002, 99, 5361. (22) Ferrari, S.; Costi, P. M.; Wade, R. C. Chem. Biol. 2003, 10, 1183. (23) Li, W.; Liu, H.; Luo, X.; Zhu, W.; Tang, Y.; Halpert, J. R.; Jiang, H. Drug Metab. Dispos. 2007, 35, 689. (24) Blondel, A.; Renaud, J. P.; Fischer, S.; Moras, D.; Karplus, M. J. Mol. Biol. 1999, 291, 101. (25) Kosztin, D.; Izrailev, S.; Schulten, K. Biophys. J. 1999, 76, 188. (26) Martinez, L.; Sonoda, M. T.; Webb, P.; Baxter, J. D.; Skaf, M. S.; Polikarpov, I. Biophys. J. 2005, 89, 2011. (27) Carlsson, P.; Burendahl, S.; Nilsson, L. Biophys. J. 2006, 91, 3151. (28) Martinez, L.; Webb, P.; Polikarpov, I.; Skaf, M. S. J. Med. Chem. 2006, 49, 23. (29) Genest, D.; Garnier, N.; Arrault, A.; Marot, C.; Morin-Allory, L.; Genest, M. Eur. Biophys. J. 2008, 37, 369. (30) Sonoda, M. T.; Martinez, L.; Webb, P.; Skaf, M. S.; Polikarpov, I. Mol. Endocrinol. 2008, 22, 1565. (31) Perakyla, M. Eur. Biophys. J. 2009, 38, 185. (32) Martinez, L.; Polikarpov, I.; Skaf, M. S. J. Phys. Chem. B 2008, 112, 10741. (33) Moras, D.; Gronemeyer, H. Curr. Opin. Cell Biol. 1998, 10, 384. (34) Jarzynski, C. Phys. ReV. Lett. 1997, 78, 2690. (35) Sali, A.; Blundell, T. L. J. Mol. Biol. 1993, 234, 779. (36) Bas, D. C.; Rogers, D. M.; Jensen, J. H. Proteins 2008, 73, 765. (37) Duan, Y.; Wu, C.; Chowdhury, S.; Lee, M. C.; Xiong, G.; Zhang, W.; Yang, R.; Cieplak, P.; Luo, R.; Lee, T.; Caldwell, J.; Wang, J.; Kollman, P. J. Comput. Chem. 2003, 24, 1999. (38) Wang, J.; Wolf, R. M.; Caldwell, J. W.; Kollman, P. A.; Case, D. A. J. Comput. Chem. 2004, 25, 1157. (39) Wang, J.; Wang, W.; Kollman, P. A.; Case, D. A. J. Mol. Graph. Model. 2006, 25, 247. (40) Case, D. A.; Darden, T. A.; Cheatham, T. E., III; Simmerling, C. L.; Wang, J.; Duke, R. E.; Luo, R.; Crowley, M.; Walker, R. C.; Zhang, W.; Merz, K. M.; Wang, B.; Hayik, S.; Roitberg, A.; Seabra, G.; Kolossva´ry, I.; Wong, K. F.; Paesani, F.; Vanicek, J.; Wu, X.; Brozell, S. R.; Steinbrecher, T.; Gohlke, H.; Yang, L.; Tan, C.; Mongan, J.; Hornak, V.; Cui, G.; Mathews, D. H.; Seetin, M. G.; Sagui, C.; Babin, V.; Kollman, P. A. AMBER 10; University of California, San Francisco, 2008. (41) Jorgensen, W. L.; Chandrasekhar, J.; Madura, J. D.; Impey, R. W.; Klein, M. L. J. Chem. Phys. 1983, 79, 926.

10444

J. Phys. Chem. B, Vol. 113, No. 30, 2009

(42) Ryckaert, J. P.; Ciccotti, G.; Berendsen, H. J. C. J. Comput. Phys. 1977, 23, 327. (43) Essmann, U.; Perera, L.; Berkowitz, M. L.; Darden, T.; Lee, H.; Pedersen, L. G. J. Chem. Phys. 1995, 103, 8577. (44) Park, S.; Khalili-Araghi, F.; Tajkhorshid, E.; Schulten, K. J. Chem. Phys. 2003, 119, 3559. (45) Park, S.; Schulten, K. J. Chem. Phys. 2004, 120, 5946. (46) Deng, Y.; Roux, B. J. Phys. Chem. B 2009, 113, 2234. (47) Jensen, M. O.; Park, S.; Tajkhorshid, E.; Schulten, K. Proc. Natl. Acad. Sci. U.S.A. 2002, 99, 6731. (48) Amaro, R.; Tajkhorshid, E.; Luthey-Schulten, Z. Proc. Natl. Acad. Sci. U.S.A. 2003, 100, 7599. (49) Zhang, D.; Gullingsrud, J.; McCammon, J. A. J. Am. Chem. Soc. 2006, 128, 3019.

Shen et al. (50) Hwang, H.; Schatz, G. C.; Ratner, M. A. J. Phys. Chem. B 2006, 110, 26448. (51) Vashisth, H.; Abrams, C. F. Biophys. J. 2008, 95, 4193. (52) Ytreberg, F. M.; Zuckerman, D. M. J. Comput. Chem. 2004, 25, 1749. (53) Ekena, K.; Weis, K. E.; Katzenellenbogen, J. A.; Katzenellenbogen, B. S. J. Biol. Chem. 1996, 271, 20053. (54) Celik, L.; Lund, J. D.; Schiott, B. Biochemistry 2007, 46, 1743. (55) Nettles, K. W.; Bruning, J. B.; Gil, G.; Nowak, J.; Sharma, S. K.; Hahm, J. B.; Kulp, K.; Hochberg, R. B.; Zhou, H.; Katzenellenbogen, J. A.; Katzenellenbogen, B. S.; Kim, Y.; Joachmiak, A.; Greene, G. L. Nat. Chem. Biol. 2008, 4, 241.

JP903785H