Importance of Receptor Conformations in Docking Calculation-Based

This is an open access article published under an ACS AuthorChoice License, which permits copying and redistribution of the article or any adaptations for non-commercial purposes.

Article Cite This: ACS Omega 2019, 4, 6620−6629

http://pubs.acs.org/journal/acsodf

Importance of Receptor Conformations in Docking CalculationBased Risk Assessment for Endocrine Disruptors against Estrogen Receptor α Hitoshi Kesamaru,† Keitaro Suyama,‡ and Takeru Nose*,†,‡ †

Department of Chemistry, Faculty and Graduate School of Science, and ‡Faculty of Arts and Science, Kyushu University, Fukuoka 819-0395, Japan

Downloaded via 85.202.194.213 on April 10, 2019 at 19:41:12 (UTC). See https://pubs.acs.org/sharingguidelines for options on how to legitimately share published articles.

S Supporting Information *

ABSTRACT: Employment of appropriate receptor conformations as templates is essential for appropriate identification of the latent receptor binding ability of chemicals using in silico docking calculations. In this study, we performed docking calculations using a number of agonistand antagonist-bound conformations of 83 estrogen receptor (ER) αligand-binding domains as templates to clarify the type of receptor conformations required for reasonable identification of endocrine disruptors. Our results showed that 17β-estradiol and diethylstilbestrol (ERα agonists) bound preferentially to the agonist conformations, whereas 4-hydroxytamoxifen and raloxifene (ERα antagonists) bound selectively to the antagonist conformations. We also observed that bisphenol A analogues, which are partial agonists, bound more moderately and preferentially to the agonist conformations as compared with the antagonist conformations. Additionally, the docking calculations were able to estimate biological agonist or antagonist activity of chemicals based on the receptor conformation selectivity. Furthermore, structural analyses of the ligand-binding domains and docking calculation utilizing C-terminaltruncated receptors indicated that the C-terminal regions of these domains were capable of discriminating agonists from nonagonists. These results suggest that both agonist- and antagonist-binding conformations of receptors are necessary to predict the binding affinity and biological activity of chemicals for docking calculation-based risk assessment. Furthermore, this in silico method can be beneficial for drug discovery because it is useful for rapid searching of ligands for receptors and preventing the side effects caused by unfavorable receptor binding.

1. INTRODUCTION Nuclear receptors (NRs) are receptor proteins involved in regulating a wide range of physiological functions. Following agonist ligand binding, NRs change their conformation to the activated form (Figure 1), while antagonist ligand binding results in a conformational change to an inactive state. These conformational changes allow NRs to differentiate chemical agonists and antagonists.1−4 Under physiological conditions, NRs are susceptible to endocrine disruption by environmental chemicals,5 including the so-called endocrine-disrupting chemicals (EDCs).6 These chemicals can bind to human NRs and cause inappropriate activation or inhibition of gene transcription, resulting in adverse biological and physiological effects involving development, reproduction, cell differentiation, nervous system function, and immune responses in mammals and especially in fetuses and infants.7−9 A representative example of an EDC is bisphenol A (BPA), which is used to manufacture polycarbonate plastics, epoxy resins, and other products. Exposure to BPA in mice and rats leads to early onset of sexual maturation in females,10−12 increased prostate size in male offspring,13−15 and altered immune function.16−18 © 2019 American Chemical Society

Therefore, to avoid health damage by BPA, several bisphenol analogues known as “next-generation bisphenols” have been developed and used as alternatives to BPA.19 However, some BPA analogues, such as bisphenol AF (BPAF),20 have been receiving much attention because these analogues also bind to NRs and could induce endocrine disruption similar to BPA. Furthermore, it is highly likely that latent chemicals will emerge as EDCs similar to BPA; therefore, all NRs are regarded as candidate targets of these yet undetermined EDCs.21 Based on this likelihood, it is important to develop reliable assays capable of predicting the risk of NR-related endocrine disruption associated with various chemicals. The recent emphasis on implementing strict ethical standards in research has increased the demand to switch from in vivo to in vitro or in silico methods for risk assessment of chemical substances. In particular, in silico methods for risk assessment have been developed in recent years. To date, many reliable and large-scale procedures for estimating the risk Received: January 7, 2019 Accepted: March 29, 2019 Published: April 10, 2019 6620

DOI: 10.1021/acsomega.9b00050 ACS Omega 2019, 4, 6620−6629

ACS Omega

Article

ligand-binding domain (LBD) complexes and three crystallographic conformations of the protein: the apo structure and a protein structure with a bound agonist or antagonist ligand.28 We have also reported previously on an agonist/antagonist differential-docking screening method using a total of four ERα-LBDs, including two agonist-bound conformations and two antagonist-bound conformations.26 In order to further advance these studies and to clarify the optimal method for selecting the receptor structures in docking calculation, we tested 9 chemicals29−33 (including ERα agonists, antagonists, and nonbinder) (Table 1) against crystal structures representing 83 human ERα-LBDs from which the docked ligand (if present) was removed and categorized them into 2 groups: agonist-bound conformations (ABC) and nonagonist-bound conformations (NABC). By statistically comparing the results of these calculations, we investigated whether docking calculations were capable of accurately predicting ERα chemical binding activity. Moreover, using modified receptor structures in which H11 and H12 were truncated, we also investigated how H11 and H12 positioning affects the docking calculation.

Figure 1. Different spatial positioning of H12 on ERα-LBDs determined by X-ray crystallography. H12s are colored in magenta (ABC; 1GWR), yellow [NABC (apo-conformation); 1A52], and green [NABC (antagonist-bound conformation); 3ERT].

of chemicals have been reported, such as quantitative structure−activity relationship and virtual screening approaches utilizing in silico methods that do not require the actual chemical compounds.22−24 Among the in silico techniques, docking calculations25,26 represent a possible first-line screening method for EDCs because of their rapidity and safety (as they eliminate the need for handling dangerous chemicals). Docking calculations require the X-ray crystallographic 3D structure of the receptor as a template. However, the quality and choice of the structures impact the results substantially. In most cases, the receptor structure used as the docking template has been treated as rigid (i.e., all atoms of the protein structure are fixed to save calculation time and increase the number of test compounds) during docking calculation using common and reliable procedures such as AutoDock, DOCK, and Autodock Vina.27 Under this rigid condition, if the receptor structure is not appropriately selected, binding affinity of the chemical will not be estimated appropriately because the conformations of each residue on the receptor structure are incapable of adapting to all binding ligands. Therefore, it is necessary to select suitable template structures to enable accurate evaluation of binding affinities between receptor and ligands via docking calculations.26 Nevertheless, it is not always possible to select an appropriate receptor structure for in silico studies. Thus, docking calculation studies using only one receptor structure may not be appropriate to evaluate the binding affinities of various chemicals. In this context, statistical analyses of the binding energy values obtained by docking calculations using multiple receptor structures might be an effective method to reduce prediction errors. In this study, we performed docking calculations using estrogen receptor alpha (ERα) as the NR model to study whether the physiological activity of ER ligands could be predicted with high accuracy by statistical analyses of the binding energy values. ERα is a member of the steroid receptor family and a representative target of EDCs.3 Furthermore, numerous agonist-bound and nonagonist-bound crystal structures of ERα/ligand complexes have been reported and deposited in the Protein Data Bank (PDB; www.rcsb.org); ERα is suitable for use in such statistical analysis of EDC binding. In silico docking studies using different conformations of ERα have been reported previously. For example, Celik et al. utilized three newly found quasistable structures of ERα and

2. RESULTS AND DISCUSSION 2.1. Analysis and Classification of Receptor Structures. Crystal structures of the ER were classified by the positioning of the H12 helix. The LBDs of NRs share a highly homologous structure, including 12 α-helices (H1 to H12) and one or more underlying β-sheets. The NR protein conformation is influenced by the positioning of the H12 helix, which depends on the biological characteristics of the ligand (agonist or antagonist). Within the LBD, a ligand-binding pocket (LBP) mainly composed of multiple hydrophobic residues and fewer hydrophilic residues is responsible for ligand binding. Upon ligand binding, the LBD attains an altered conformation of H12 depending upon the structure of the ligand bound in the LBP, with these alterations observed in H12 that is present in the C-terminal region (Figure 2a). H12 conformation is primarily categorized as representing either an activated conformation or a nonactivated conformation. When an agonist binds to wild-type ERα, the amide hydrogen atom of Leu540 on H12 forms a hydrogen bond with the oxygen atom of the Asp351 side chain on H3.34 Additionally, in the case of the Y537S mutant of ERα, it was reported that Asp351 formed a hydrogen bond with Ser537, forcing it to form an activation conformation.35 Therefore, ABC and NABC were classified based on whether the distance between Leu540 and Asp351 for the wild-type ERα and Ser537 and Asp351 for the Y537S mutant ERα was within a hydrogen bond distance of ∼3.5 Å or not. By applying these criteria to 83 ERα-LBDs (Table S1 in the Supporting Information), we assigned 56 ERα-LBDs to the ABC group and 27 ERα-LBDs to the NABC group. The ratio of ERα-LBD structures in the ABC group to those in the NABC group was almost 2:1. Therefore, we assumed that analysis of all LBDs using docking calculations would result in a higher identification rate of high-affinity chemical agonists as compared with chemicals lacking agonist-binding capability. Thus, it was concerned that calculation results by using ABC and NABC should be evaluated separately. We also further observed each of the ERα-LBPs in terms of LBP volume, as it was reported that differences in LBP volume could affect ligand selectivity in vitro.2,36,37 Although it is not always possible to identify residues within the loop region 6621


ACS Omega

Article

Table 1. Chemical Structure of Compounds Utilized in This Study

between H11 and H12 in NABC structures that form a part of the pocket wall of the LBP, it was estimated that the average LBP volume in structures from the ABC (55 structures) and NABC (10 structures) groups were 632 ± 67 and 707 ± 89 Å3,

respectively. Thus, structures from the NABC group had larger LBPs than those from the ABC group. This difference in LBP volume might provide a reason for differences in ligand selectivity observed in in silico docking calculations. 6622


ACS Omega

Article

group, the Kdw/Kdnw ratio was 0.97 in the presence of E2 and 1.5 in the presence of (Z)-4-OHT. For structures in the NABC group, the Kdw/Kdnw ratio was 0.99 for E2 and 1.5 for (Z)-4OHT, indicating highly similar ratios between the two groups. E2 showed relatively high affinity for structures in the ABC group, whereas (Z)-4-OHT showed high affinity for structures in the NABC group. These results suggested that the presence of water molecules in and around the LBD did not significantly influence the results of docking calculations performed on structures in either the ABC or NABC group. Thus, water molecules in the LBP crystal structures were preserved in the subsequent analysis. 2.3. Comparison of Ligand Binding Energy between ABC and NABC ERα-LBD Structures. The binding affinities of nine chemicals, including agonists, antagonists, and partial agonists, to ERα-LBD structures were investigated by docking calculation using the structure classification mentioned above (Table 4). Typical ERα agonists [E2 and diethylstilbestrol (DES)] and antagonists [(Z)-4-OHT and raloxifene (RAL)] were considered as positive controls in the docking calculation. In contrast, BPA, as an EDC, and bisphenol C (BPC), BPAF, and 2,2-bis(p-hydroxyphenyl)-1,1,1-trichloroethane (HPTE), as EDC candidates for ERα, were further assessed for their receptor binding strength and activation potency. Using the average calculated ΔG against each test chemical listed in Table 4, we obtained the dissociation constant (Kd) values associated with their receptor binding activity against the target structures by using the abovementioned theoretical formula i (Table 5). To confirm the reproducibility of the docking calculation, the calculated average ΔG values were plotted against log Kd values previously determined using wet assays, and Pearson’s product−moment correlation (r) was calculated.38 A strong correlation (r = 0.737) was observed between the calculated average ΔG values and Kd values determined by wet assays (Figure S2 in the Supporting Information). Among the tested chemicals, biphenyl, a nonbinder for ERα, showed significantly weaker binding affinity than the other chemicals (Kd = 10 μM) as expected.33 Therefore, the following discussion is related to the eight chemicals, except for biphenyl, that were expected to bind to ERα. The statistical difference between ΔG values from the ABC group and those from the NABC group was determined by Welch’s t-test. The ΔG values of agonist ligand binding for the ABC group were found to be statistically smaller than those observed in the NABC group, whereas the ΔG of antagonist binding for the NABC group was smaller than that observed in the ABC group. Based on these values, Kd ratios, that is the Kdnabc/Kdabc, were calculated for each chemical to investigate the selectivity of the conformations that indicated which receptor group preferentially bound to the chemical (Table 5). Usually, E2 (a potent mammalian estrogenic steroid hormone) binds as an agonist to induce an activated conformation in ERα. Our results showed that, on average, E2 bound to ERα structures in the ABC group with 1.7-fold greater affinity as compared with structures in the NABC group. Although E2 and DES share similar molecular features, including molecular volumes (E2 = 234.51 Å3; DES = 245.02 Å3), two characteristic hydroxyl groups, and ERα agonistic activity, DES exhibited only 1.2-fold greater affinity for ERα structures in the ABC group than it did for those in the NABC group. In contrast with the antagonists, we observed no significant differences in agonist binding between ABC and NABC groups (Table 5). Additionally, the calculated DES binding affinity for

Figure 2. Representative position of H12 in ABC. (a) H12 is indicated by an arrow. (b) Hydrogen bonds formed between Leu540 and Asp351 in the ABC of wild-type ERα-LBD. (c) Hydrogen bonds formed between Ser537 and Asp351 in the Y537S mutant.

2.2. Effects of Water Molecules on LBP/Ligand Docking Calculations. To investigate the effects of water molecules, the docking calculations were performed using the same LBDs in the presence or absence of water molecules, as they play roles as spacers or mediators of hydrogen bonds between a ligand and residues of the LBP. There was no change in the calculation condition other than the presence or absence of water molecules (Table 2). The dissociation Table 2. Calculated Binding Energy of E2 and (Z)-4-OHT to ERα in the Presence or Absence of Water Moleculesa ΔG (kcal/mol) water preserved ABC NABC

−10.2 −10.8 −9.89 −11.5

E2 (Z)-4-OHT E2 (Z)-4-OHT

± ± ± ±

water removed −10.2 −11.0 −9.89 −11.8

0.1 0.1 0.07 0.1

± ± ± ±

0.1 0.1 0.06 0.1

Values are mean ± SE.

a

constant (Kd) values associated with receptor binding activity of chemicals against the target structures were calculated by using the theoretical formula i ΔG = RT ln Kd

(i)

Ratios of the two calculated Kd values [Kdw (waterpreserved) and Kdnw (water-removed)] are shown in Table 3, and in the calculations, 17β-estradiol (E2) and (Z)-4-OHT were used as representative ligands. For structures in the ABC Table 3. Calculated Dissociation Constant of the Binding of E2 and (Z)-4-OHT to ERα in the Presence or Absence of Water Moleculesa,b,c Kd (nM) water preserved ABC NABC

E2 (Z)-4-OHT E2 (Z)-4-OHT

32.7 11.8 56.3 3.58

± ± ± ±

3.3 1.2 3.9 0.4

water removed

Kdnw/Kdw

± ± ± ±

0.97 1.5 0.99 1.5

33.8 8.02 57.1 2.45

3.4 0.80 3.4 0.25

Values are mean ± SE. bKdw: dissociation constant value of the binding of a chemical to ERα-LBDs with water molecules calculated from ΔG. cKdnow: dissociation constant value of the binding of a chemical to ERα-LBDs without water molecules calculated from ΔG. a

6623


ACS Omega

Article

Table 4. Calculated Binding Energy between Test Chemical and ERαa ΔG (kcal/mol) all

ABC

mean ± SE

chemicals

−10.1 −8.80 −11.1 −11.5 −7.95 −8.31 −7.36 −8.58 −6.48

E2b DES (Z)-4-OHTb RALb BPAb BPCb BPAFc HPTEb biphenyl

± ± ± ± ± ± ± ± ±

mean ± SE

SD

0.1 0.04 0.1 0.2 0.03 0.03 0.03 0.04 0.04

−10.2 −8.84 −10.8 −10.7 −8.01 −8.38 −7.41 −8.66 −6.55

0.46 0.36 0.91 1.9 0.26 0.31 0.28 0.33 0.42

± ± ± ± ± ± ± ± ±

0.1 0.07 0.2 0.3 0.04 0.05 0.04 0.06 0.08

NABC SD 0.47 0.37 0.94 1.8 0.23 0.28 0.23 0.32 0.46

mean ± SE −9.90 −8.73 −11.5 −13.2 −7.83 −8.17 −7.24 −8.42 −6.33

± ± ± ± ± ± ± ± ±

0.07 0.06 0.1 0.2 0.06 0.06 0.06 0.06 0.05

SD 0.37 0.33 0.57 0.78 0.29 0.32 0.33 0.30 0.24

Mean: means of calculated ΔG value. SE: standard error. SD: standard deviation. bP < 0.05 in Welch’s t-test between ΔG values of ABC and NABC. c0.05 < P < 0.1 in Welch’s t-test between ΔG values of ABC and NABC.

a

shows 160 000-fold lower induction of in vivo estrogenic activity against ethynyl-estradiol following oral administration.43 In this study, we estimated that the BPA binding affinity for all structures were 39-fold lower than that of E2. This estimated weak binding affinity agreed with the previous results cited above and also suggested that BPA is a very weak or partial ERα-agonist. According to Kd ratio analysis, BPA exhibited 1.4-fold higher affinity for the ABC group relative to those in the NABC group (Table 5), indicating that BPA is not a potent or complete ERα-agonist. Other bisphenol analogues, including BPC, BPAF, and HPTE, bound to structures in the ABC group with 1.4-, 1.3-, and 1.5-fold higher affinity relative to their binding to structures in the NABC group, respectively (Table 5). These results confirmed that the BPA analogues may act as ERα-agonists, in agreement with previous in vitro results.20,30,31 In this study, we used RAL and (Z)-4-OHT as representative antagonists in the chemicals tested in silico. RAL bound to ERα structures in the NABC group with 68-fold higher affinity relative to its binding to structures in the ABC group. We also observed that the antagonist (Z)-4-OHT bound to ERα structures in the NABC group with 3.3-fold higher affinity relative to its binding to structures in the ABC group. The observed NABC selectivity may be due to the large LBP volumes of the NABC structures, which may be able to properly accept relatively large antagonist ligands. Overall, our docking calculation studies predicted that ERα agonists preferentially bound to the ABC group of the LBD structure, whereas antagonist ligands showed stronger binding affinity to NABC than ABC structures. Because these results were consistent with previous reports,26 the necessity of proper selection of ABC/NABC for docking calculation was appropriately confirmed in this study. In addition, we suggest that the possibility of transcriptional activity of unknown chemicals might be predicted by investigating whether the test ligand could differentiate ABC and NABC based on the docking calculation. 2.4. Importance of C-Terminal Truncation on ERαLBD Docking Calculations. The C-terminal region of ERαLBDs demonstrated an apparent variety of conformations between structures (Figure 1). To investigate the effects of conformational variance on in silico docking calculations, we compared the root-mean-square deviation (RMSD) values of all full-length ERα-LBD structures (Table 6). The mean RMSD for all full-length ERα-LBD structures was 2.55 ± 0.03

Table 5. Calculated Dissociation Constant between Test Chemical and ERαa Kd (nM) chemicals E2b DES (Z)-4-OHTb RALb BPAb BPCb BPAFc HPTEb biphenyl

ABC

NABC

mean ± SE

mean ± SE

Kdnabc/Kdabc

± ± ± ± ± ± ± ± ±

1.7 1.2 0.31 0.014 1.4 1.4 1.3 1.5 1.4

33.6 333 12.2 14.5 1352 724 3721 452 15 879

± ± ± ± ± ± ± ± ±

3.4 23 2.4 4.3 54 36 149 27 1270

55.8 401 3.75 0.213 1832 1032 4958 677 23 016

3.9 24 0.37 0.043 110 62 297 41 1151

Mean: means of calculated ΔG value. SE: standard error. bP < 0.05 in Welch’s t-test between ΔG values of ABC and NABC. c0.05 < P < 0.1 in Welch’s t-test between ΔG values of ABC and NABC. a

structures in the ABC group (ΔG = −8.84 kcal/mol) was weaker than that for E2 (ΔG = −10.2 kcal/mol), despite the similar 50% effective concentrations (EC50; E2 = 32 pM; DES = 13 pM).39 The complex structure that was determined by Xray crystallography of DES and ERα-LBD (PDB: 3ERD) was clearly different from the calculated complexes using ABC in this study. Among the amino acid residues of ERα-LBDs, Met421 and His524, His524 of 3ERD changed its conformation to preferably interact with the two linearly oriented hydroxyl groups of DES (Figure S3). Thus, this discrepancy in the calculated binding affinity may be caused by the lack of a specific DES-accepting conformation of ERα-LBD. It was confirmed from the crystal structure of DES/ERα-LBD (PDB: 3ERD) that the conformation of H11 of ERα-LBD was similar to that found in NABC rather than that in ABC, although the position of H12 changed to the activated form found in ABC. The ERα-LBD can induce such a unique conformation when DES binding occurs. In the present study, the crystal structures of the ERα-LBD were simply classified into two groups (ABC and NABC) by the positioning of the H12 helix to apply the docking calculation. Therefore, in case of chemicals that induce a unique binding conformation of ERα-LBDs, such as DES, the accuracy of the activity prediction method proposed in this study might be limited by unconventionality. BPA is a known ERα ligand; however, its binding affinity is 500- to 15 000-fold lower than that of E2,14,40−42 and it also 6624


ACS Omega

Article

calculation of agonist ligands (E2 and BPA), ΔΔG values slightly increased compared to those with full-length structures. On the other hand, in the presence of antagonist ligands ((Z)4-OHT and RAL), ΔΔG with the ΔK529C variants significantly decreased compared to those with the full-length structures. Decreasing of ΔΔG values could indicate that the selectivity of antagonist ligands for the ΔK529C variants was lower than that for the full-length structures. Thus, H12 seems to be important for ligand selectivity particularly in case of antagonists. Additional docking calculations were performed using the ΔN519C series to determine the importance of H11, in addition to H12. We observed that ΔΔG with the ΔN519C variants from both the ABC and NABC groups showed low ligand selectivity for E2 and BPA (ΔΔG = 0.02 and 0.10 kcal/ mol, respectively). Therefore, H11 may be important for ligand selectivity in the case of agonists. By contrast, in the presence of (Z)-4-OHT and RAL, ligand selectivity was lower in both ΔK529C and ΔN519C than in the full-length LBD (Figure 4). The relative decrease in the ABC/NABC selectivity of ligands using the C-terminal-truncated LBDs in the docking calculation indicates that ligand selectivity is chiefly determined by the ERα C-terminal region (Figure 4). Therefore, the docking calculations were unable to estimate ligand selectivity in the absence of the region that included H11 and H12 of ERα-LBD (from N519 to the C-terminus). For further evaluation of the importance of H12 and H11 on ligand orientation in the docking calculation, docked conformations of E2, (Z)-4-OHT, and RAL were analyzed. For analysis, the complex structures between ligands and full length, ΔK529C, or ΔN519C categorized in ABC (1GWR) and NABC (3ERT and 1ERR) were utilized (Figure 5). By the deletion of both H11 and H12 of 1GWR, binding conformation of E2 was changed. In contrast, when H11 was preserved, conformation of E2 was almost identical to that observed with the full-length LBD. Therefore, the binding conformation of E2 may be mainly affected by H11. On the other hand, conformations of antagonist ligands, (Z)-4-OHT and RAL, exhibited only little change in the NABC regardless of existence of H11 and H12 (Figure 6). Furthermore, binding conformation of (Z)-4-OHT and RAL in the truncated LBDs belonging to ABC, 1GWRΔK529C and 1GWRΔN519C, were almost similar to those in the NABC group. Thus, the conformation of antagonist ligands in the LBDs is affected depending on the orientation or existence of H12. 2.5. Importance of Docking Calculation Using Multiple LBDs. In this study, we performed the docking calculation using 83 ERα-LBDs as templates, which were categorized into 2 groups, ABC and NABC. In a docking calculation, the use of only one receptor LBD that has been determined as a

Table 6. RMSD Analysis for Each Receptor Conformation Group of ERαa comparison subjects all conformations ABC vs ABC NABC vs NABC ABC vs NABC

RMSD (full-length LBD) (Å) 2.55 1.22 1.40 4.09

± ± ± ±

0.03 0.01 0.02 0.03

Values are mean ± SE.

a

Å, whereas the values for structures within the ABC and NABC groups were 1.22 ± 0.01 and 1.40 ± 0.02 Å, respectively. As expected, a large RMSD value (4.09 ± 0.03 Å) was derived from comparisons between the ABC and NABC groups. Furthermore, to consider the importance of the C-terminal region in docking calculation, we prepared two in silico Cterminal-curtailing mutants involving truncation at K529 (ΔK529C) and another at N519 (ΔN519C). As shown in Figure 3, N519 is located in the middle of H11, and K529 is

Figure 3. Spatial positioning of Asn519 and Lys529 on ERα-LBD (PDB: 1GWR). (a) Whole LBD structure and (b) enlarged view around H11 and H12 are shown. Asn519 and Lys529 are shown by stick models. Asn519 is located on H11, whereas Lys529 is located on the loop between H11 and H12.

located at the end of H11 in close proximity to the loop between H11 and H12 (Figure 3). A series of ΔK529C (lacking the H12 moiety) and ΔN519C (lacking the H11 and H12 moiety) mutants was prepared for all 83 ERα-LBDs present in both the ABC and NABC groups. To evaluate the effects of the C-terminal region of the LBDs with regard to ligand binding, docking calculations of E2, BPA, (Z)-4-OHT, and RAL were performed using ΔK529C series as templates. Using the obtained ΔG values (Figure 4), ΔΔG was defined as the difference in ΔG between ABC and NABC. When the ΔK529C variants were used as templates for docking

Figure 4. Comparison of binding energies of the ligands for the full-length LBD (A), ΔK529C (B), and ΔN519C (C). Ligands indicated are: E2, BPA, (Z)-4-OHT, and RAL. *P < 0.05 by Welch’s t-test between ΔG values of ABC and NABC. 6625


ACS Omega

Article

Figure 5. Predicted conformation of E2, (Z)-4-OHT, and RAL bound to the ABC. Predicted conformations of E2 with (a) full-length (magenta), (b) ΔK529C (cyan), and (c) ΔN519C (green) of 1GWR-LBD are shown. Ligand conformation from PDB is colored in gray (PDB: 1GWR). Predicted conformations of (Z)-4-OHT with (d) full-length (magenta), (e) ΔK529C (cyan), and (f) ΔN519C (green) of 1GWR-LBD are shown. Ligand conformation from PDB is colored in gray (PDB: 3ERT). Predicted conformations of RAL docked with (g) full-length (magenta), (h) ΔK529C (cyan), and (i) ΔN519C (green) 1GWR-LBD are shown. Ligand conformation from PDB is colored in gray (PDB: 1ERR).

conformation of the LBD. However, it also means that LBD structures determined as ligand complexes optimized their structures specifically for each bound ligand. Because such an LBD structure optimized for a specific ligand may not be suitable for docking of various ligands, it is necessary to use multiple structures for docking calculations to perform statistical analyses. It would also be valuable to know the complete structure of the solution state of the LBD as a docking template, but this is currently only presumed in in silico calculations. One of the future challenges is to establish a statistical docking calculation strategy using apo-LBDs with molten globule structures whose structure is difficult to determine. As shown in this study, using a number of LBDs, one can aid the appropriate estimation of receptor agonists in docking calculation with statistical analysis. Thus, to perform effective docking calculations, it is important to increase the number of receptor structures used as templates in the future.

Figure 6. Predicted conformation of (Z)-4-OHT and RAL bound to the NABC. (a) Predicted conformations of (Z)-4-OHT in full-length (magenta) ΔK529C (cyan) and ΔN519C (green) of 3ERT-LBD are shown. Ligand conformation from PDB is colored in gray. (b) Predicted conformations of RAL docked with full-length (magenta), ΔK529C (cyan), and ΔN519C (green) 1ERR are shown. Ligand conformation from PDB is colored in gray.

receptor/ligand complex can cause problems. For example, the LBD determined as a receptor/small agonist complex might not afford the binding of a large antagonist due to the volume shortage of the LBP. This false recognition makes it difficult to identify seed compounds in drug discovery and EDCs in the risk assessment of chemicals. Thus, it was considered that docking calculation using multiple templates would be an effective method to reduce errors in ligand activity prediction. Celik et al. demonstrated a docking calculation using multiple LBDs, including three structures constructed by molecular dynamics simulation: an agonist-bound, an antagonist-bound, and an apo-form LBD determined by X-ray crystallography.28 They described that docking into the quasistable conformations of the ERα-LBDs leads to computed binding affinities similar to or better than those calculated for the three X-ray structures. Based on this result, they emphasized that the new structures may be important for assessing the influence of EDCs on NRs. We also reported that the chemicals can be differentiated into either agonists or antagonists by docking calculations of the chemical for the LBD, including two agonist-bound and two antagonist-bound LBDs.26 Recently, a number of NR-LBD/ligand complex structures have been determined by X-ray crystallography. Regarding the structure of NR-LBDs, there is a report that a ligand (rosiglitazone)bound PPARγ-LBD is more stable than the apo-LBD.44 This result suggested that the ligand binding induced a stable

3. CONCLUSIONS In this study, we performed in silico ligand docking studies using multiple ERα-LBDs as templates to assess whether the physiological activity of ER ligands could be predicted with high accuracy by statistical analyses of the binding energy values to ABC and NABC. In particular, we focused on and tested bisphenol-based chemicals in the docking calculation and found weak agonistic activity of the bisphenol-based chemicals via ERα. Our results suggested that receptor binding affinity and potential bioactivity for various bisphenol analogues could be estimated properly by the calculation. Moreover, we revealed that the C-terminal regions of ERαLBDs involving H11, H12, and an intermediate loop structure were important for estimation of ligand affinity by the docking calculation. Our study highlights that for accurate discrimination of EDCs that exhibit receptor binding properties and presumed biological activity, it is necessary to analyze results of the docking calculation using different receptor conformations. Acquisition of both ABC and NABC structures of NRs will enable the development of in silico docking calculation-based risk assessments capable of EDC identification. Because the in silico method presented here can rapidly detect receptor 6626


ACS Omega

Article

maximum number of energy evaluations was 2.5 × 107, and the maximum number of generations was 1.0 × 106. Other parameters were represented by default values implemented by the program. In the calculation, the free energy change (ΔG) associated with each conformation was defined as the sum of free energy changes, van der Waals forces, electrostatic interactions, hydrogen bonding events, desolvation activity, and torsion energetics as designed by the developer. The dissociation constant (Kd) was calculated from ΔG using a theoretical formula i

agonists, it is useful not only for computational toxicology research but also for pharmaceutical studies aiming to discover novel drug candidates.

4. EXPERIMENTAL SECTION 4.1. Preparation of Receptor Structures for Docking Calculations. As templates for docking calculations, 83 ERαLBD structures reported by 2015 were prepared using crystal structures deposited in the PDB (see Table S1 in the Supporting Information). When preparing the 3D structures for use as docking templates, hydrogen atoms were added to 83 structures and optimized using Discovery Studio software (v4.0; Accelrys, San Diego, CA, USA). All calculations were performed using a Dell Precision T3610 workstation (Dell, Round Rock, TX, USA). In all cases, cocrystallized ligands and other small molecules were removed. With regard to water molecules in the crystals, both water-removed and waterpreserved LBDs were prepared for further analysis. The size of the LBP in ERα-LBDs and RMSD between receptor structures was evaluated using Discovery Studio (Accelrys). Receptor structures were classified as ABC or NABC (including both the antagonist-bound and apo conformations). Receptor structures that belong to ABC were categorized on the condition that the distance between the amide nitrogen of Leu540 and the oxygen atom on Asp351 side chain in wild-type ERα, or alternatively Ser357 and Asp351 in the Y537S mutant, was within hydrogen bond distance (∼3.5 Å). On the other hand, NABC was composed of receptor structures other than that classified as ABC. 4.2. In Silico Preparation of C-Terminal-Truncated ERα-LBDs. To assess the effects of deletions of the C-terminal regions of LBDs on ligand selectivity, C-terminal-truncated ERα-LBDs were prepared using Discovery Studio (Accelrys). A series of ΔK529C mutations, resulting in deletion of H12 (from K529 to the C-terminus), were computationally prepared for 83 ERα-LBD structures. Each prepared structure was energy-minimized by Discovery Studio (Accelrys). A series of ΔN519C mutants that deleted both H11 and H12 (from N519 to the C-terminus) were prepared in the same manner. 4.3. In Silico Ligand Preparation for Docking Calculation. The ligands tested were as follows: E2 and DES as agonists for ERα, (Z)-4-hydroxytamoxifen [(Z)-4OHT] and RAL as antagonists for ERα, BPA and BPC as partial agonists for ERα, BPAF and HPTE as EDCs act as agonists for ERα, and 1,1′-biphenyl (biphenyl) as a nonbinder for ERα. All chemical structures of the ligands were constructed using Accelrys Draw (v4.0; Accelrys), and energy-minimized using the minimization protocol of Discovery Studio (Accelrys). Ligand volume was evaluated using the analytical tools included in Discovery Studio (Accelrys). The list of chemical structures used in this study is provided in Table 1.29−33 4.4. Docking Calculations. The binding affinities of ligands were evaluated by docking calculations using AutoDock version 4.2 (Scripps Research Institute, San Diego, CA, USA) on a Dell Precision T3500 workstation (Dell). The grid maps representing the receptor molecule were generated by AutoGrid (Scripps Research Institute), with each grid centered at the LBP. The docking calculation area was set within 126 × 126 × 126 grids at 0.197 Å per grid cube to include the entire LBP. Docking calculations were performed 100 times using the algorithm Lamarckian GA.45 The run parameters used in this study were as follows: the number of GA runs was 100, the

ΔG = RT ln Kd

(i)

The structure resulting from 100 rounds of free energy minimization was adopted as the calculation result. The lowest binding energy value among the calculation results from the 100 runs was used for the analysis. A total of 83 ERα structures (Table S1 in the Supporting Information) and 9 ligands (Table 1) were used in this study. The calculation results were analyzed and compared according to their ΔG values between ABCs and NABCs. Data normalization was determined using a normal quantile−quantile plot (Figure S1 in the Supporting Information).46 Subsequently, comparison of the means and standard errors of the mean were performed using Welch’s ttest.47

■

ASSOCIATED CONTENT

* Supporting Information S

The Supporting Information is available free of charge on the ACS Publications website at DOI: 10.1021/acsomega.9b00050. 3D structures of ERα used in this study; results of the Q−Q plot for estimation of average binding energy ΔG; Pearson’s product−moment correlation between calculated ΔG and log Kd calculated from reported Kd; and comparison of the structures between calculated and Xray crystallographically determined DES-docked ERα complexes (PDF).

■

AUTHOR INFORMATION

Corresponding Author

*E-mail: [email protected]. Phone/Fax: +81-92-8026025. ORCID

Takeru Nose: 0000-0001-9771-7267 Author Contributions

The manuscript was written through contributions of all the authors. All the authors have given approval to the final version of the manuscript. Funding

This work was supported by JSPS KAKENHI grant number JP15H02827. Notes

The authors declare no competing financial interest.

■

ABBREVIATIONS (Z)-4-OHT, 4-hydoroxytamoxifen; ABC, agonist-bound conformations; AF-1/2, activation function 1/2; BPX, bisphenol X (where X is any of the following A, AF, C, and E); DES, diethylstilbestrol; E2, 17β-estradiol; EDCs, endocrine-disrupting chemicals; ER, estrogen receptor; H12, α-helix numbered 12; HPTE, 2,2-bis(p-hydroxyphenyl)-1,1,1-trichloroethane; 6627


ACS Omega

Article

production, proliferative responses of lymphoid cells, and TH1 and TH2 immune responses in mice. Br. J. Pharmacol. 2003, 138, 1271− 1276. (18) Yoshino, S.; Yamaki, K.; Li, X.; Sai, T.; Yanagisawa, R.; Takano, H.; Taneda, S.; Hayashi, H.; Mori, Y. Prenatal exposure to bisphenol A up-regulates immune responses, including T helper 1 and T helper 2 responses, in mice. Immunology 2004, 112, 489−495. (19) Sartain, C. V.; Hunt, P. A. An old culprit but a new story: bisphenol A and “NextGen” bisphenols. Fertil. Steril. 2016, 106, 820− 826. (20) Matsushima, A.; Liu, X.; Okada, H.; Shimohigashi, M.; Shimohigashi, Y. Bisphenol AF Is a Full Agonist for the Estrogen Receptor ERα but a Highly Specific Antagonist for ERβ. Environ. Health Perspect. 2010, 118, 1267−1272. (21) Tabb, M. M.; Blumberg, B. New modes of action for endocrinedisrupting chemicals. Mol. Endocrinol. 2006, 20, 475−482. (22) Vedani, A.; Dobler, M.; Smieško, M. VirtualToxLab - A platform for estimating the toxic potential of drugs, chemicals and natural products. Toxicol. Appl. Pharmacol. 2012, 261, 142−153. (23) Zhang, L.; Sedykh, A.; Tripathi, A.; Zhu, H.; Afantitis, A.; Mouchlis, V. D.; Melagraki, G.; Rusyn, I.; Tropsha, A. Identification of putative estrogen receptor-mediated endocrine disrupting chemicals using QSAR- and structure-based virtual screening approaches. Toxicol. Appl. Pharmacol. 2013, 272, 67−76. (24) Vedani, A.; Dobler, M.; Lill, M. A. Combining protein modeling and 6D-QSAR. simulating the binding of structurally diverse ligands to the estrogen receptor. J. Med. Chem. 2005, 48, 3700−3703. (25) Luty, B. A.; Wasserman, Z. R.; Stouten, P. F. W.; Hodge, C. N.; Zacharias, M.; McCammon, J. A. A molecular mechanics/Grid method for evaluation of ligand-receptor interactions. J. Comput. Chem. 1995, 16, 454−464. (26) Nose, T.; Tokunaga, T.; Shimohigashi, Y. Exploration of endocrine-disrupting chemicals on estrogen receptor α by the agonist/antagonist differential-docking screening (AADS) method: 4-(1-Adamantyl)phenol as a potent endocrine disruptor candidate. Toxicol. Lett. 2009, 191, 33−39. (27) Pagadala, N. S.; Syed, K.; Tuszynski, J. Software for molecular docking: a review. Biophys. Rev. 2017, 9, 91−102. (28) Celik, L.; Lund, J. D. D.; Schiøtt, B. Exploring Interactions of Endocrine-Disrupting Compounds with Different Conformations of the Human Estrogen Receptor α Ligand Binding Domain: A Molecular Docking Study. Chem. Res. Toxicol. 2008, 21, 2195−2206. (29) Kuiper, G. G. J. M.; Carlsson, B.; Grandien, K.; Enmark, E.; Häggblad, J.; Nilsson, S.; Gustafsson, J. Comparison of the ligand binding specificity and transcript tissue distribution of estrogen receptors α and β. Endocrinology 1997, 138, 863−870. (30) Sun, J.; Huang, Y. R.; Harrington, W. R.; Sheng, S.; Katzenellenbogen, J. A.; Katzenellenbogen, B. S. Antagonists selective for estrogen receptor α. Endocrinology 2002, 143, 941−947. (31) Delfosse, V.; Grimaldi, M.; Pons, J.-L.; Boulahtouf, A.; le Maire, A.; Cavailles, V.; Labesse, G.; Bourguet, W.; Balaguer, P. Structural and mechanistic insights into bisphenols action provide guidelines for risk assessment and discovery of bisphenol A substitutes. Proc. Natl. Acad. Sci. U.S.A. 2012, 109, 14930−14935. (32) Gaido, K. W.; Leonard, L. S.; Maness, S. C.; Hall, J. M.; McDonnell, D. P.; Saville, B.; Safe, S. Differential interaction of the methoxychlor metabolite 2,2-Bis-(p-Hydroxyphenyl)-1,1,1-trichloroethane with estrogen receptors α and β. Endocrinology 1999, 140, 5746−5753. (33) Paris, F.; Balaguer, P.; Térouanne, B.; Servant, N.; Lacoste, C.; Cravedi, J.-P.; Nicolas, J.-C.; Sultan, C. Phenylphenols, biphenols, bisphenol-A and 4-tert-octylphenol exhibit α and β estrogen activities and antiandrogen activity in reporter cell lines. Mol. Cell. Endocrinol. 2002, 193, 43−49. (34) Kim, J. H.; Lee, M. H.; Kim, B. J.; Kim, J. H.; Han, S. J.; Kim, H. Y.; Stallcup, M. R. Role of aspartate 351 in transactivation and active conformation of estrogen receptor α. J. Mol. Endocrinol. 2005, 35, 449−464.

IC50, 50% inhibition concentration; LBD, ligand-binding domain; LBP, ligand-binding pocket; NABC, nonagonistbound conformations; NR, nuclear receptor; RAL, raloxifene

■

REFERENCES

(1) Heldring, N.; Pike, A.; Andersson, S.; Matthews, J.; Cheng, G.; Hartman, J.; Tujague, M.; Ström, A.; Treuter, E.; Warner, M.; Gustafsson, J.-Å. Estrogen receptors: how do they signal and what are their targets. Physiol. Rev. 2007, 87, 905−931. (2) Brzozowski, A. M.; Pike, A. C. W.; Dauter, Z.; Hubbard, R. E.; Bonn, T.; Engström, O.; Ö hman, L.; Greene, G. L.; Gustafsson, J.-Å.; Carlquist, M. Molecular basis of agonism and antagonism in the oestrogen receptor. Nature 1997, 389, 753−758. (3) Germain, P.; Staels, B.; Dacquet, C.; Spedding, M.; Laudet, V. Overview of nomenclature of nuclear receptors. Pharmacol. Rev. 2006, 58, 685−704. (4) Olefsky, J. M. Nuclear receptor minireview series. J. Biol. Chem. 2001, 276, 36863−36864. (5) Rouiller-Fabre, V.; Guerquin, M. J.; N’Tumba-Byn, T.; Muczynski, V.; Moison, D.; Tourpin, S.; Messiaen, S. b.; Habert, R.; Livera, G. Nuclear receptors and endocrine disruptors in fetal and neonatal testes: a gapped landscape. Front. Endocrinol. 2015, 6, 58. (6) Janošek, J.; Hilscherová, K.; Bláha, L.; Holoubek, I. Environmental xenobiotics and nuclear receptorsInteractions, effects and in vitro assessment. Toxicol. In Vitro 2006, 20, 18−37. (7) WHO-UNEP. State of the Science of Endocrine Disrupting Chemicals (2012); Bergman, A., Heindel, J. J., Jobling, S., Kidd, K. A., Zoeller, R. T., Eds.; WHO Press: Geneva, Switzerland, 2013; pp 1− 260. (8) vom Saal, F. S.; Cooke, P. S.; Buchanan, D. L.; Palanza, P.; Thayer, K. A.; Nagel, S. C.; Parmigiani, S.; Welshons, W. V. A physiologically based approach to the study of bisphenol A and other estrogenic chemicals on the size of reproductive organs, daily sperm production, and behavior. Toxicol. Ind. Health 1998, 14, 239−260. (9) Markey, C. M.; Luque, E. H.; Munoz de Toro, M.; Sonnenschein, C.; Soto, A. M. In Utero Exposure to Bisphenol A Alters the Development and Tissue Organization of the Mouse Mammary Gland. Biol. Reprod. 2001, 65, 1215−1223. (10) Honma, S.; Suzuki, A.; Buchanan, D. L.; Katsu, Y.; Watanabe, H.; Iguchi, T. Low dose effect of in utero exposure to bisphenol A and diethylstilbestrol on female mouse reproduction. Reprod. Toxicol. 2002, 16, 117−122. (11) Howdeshell, K. L.; Hotchkiss, A. K.; Thayer, K. A.; Vandenbergh, J. G.; vom Saal, F. S. Exposure to bisphenol A advances puberty. Nature 1999, 401, 763−764. (12) Nikaido, Y.; Yoshizawa, K.; Danbara, N.; Tsujita-Kyutoku, M.; Yuri, T.; Uehara, N.; Tsubura, A. Effects of maternal xenoestrogen exposure on development of the reproductive tract and mammary gland in female CD-1 mouse offspring. Reprod. Toxicol. 2004, 18, 803−811. (13) Gupta, C. Reproductive malformation of the male offspring following maternal exposure to estrogenic chemicals. Proc. Soc. Exp. Biol. Med. 2000, 224, 61−68. (14) Nagel, S. C.; vom Saal, F. S.; Thayer, K. A.; Dhar, M. G.; Boechler, M.; Welshons, W. V. Relative binding affinity-serum modified access (RBA-SMA) assay predicts the relative in vivo bioactivity of the xenoestrogens bisphenol A and octylphenol. Environ. Health Perspect. 1997, 105, 70−76. (15) Timms, B. G.; Howdeshell, K. L.; Barton, L.; Bradley, S.; Richter, C. A.; vom Saal, F. S. Estrogenic chemicals in plastic and oral contraceptives disrupt development of the fetal mouse prostate and urethra. Proc. Natl. Acad. Sci. U.S.A. 2005, 102, 7014−7019. (16) Sawai, C.; Anderson, K.; Walser-Kuntz, D. Effect of bisphenol A on murine immune function: modulation of interferon-gamma, IgG2a, and disease symptoms in NZB X NZW F1 mice. Environ. Health Perspect. 2003, 111, 1883−1887. (17) Yoshino, S.; Yamaki, K.; Yanagisawa, R.; Takano, H.; Hayashi, H.; Mori, Y. Effects of bisphenol A on antigen-specific antibody 6628


ACS Omega

Article

(35) 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.; Joachimiak, A.; Greene, G. L. NFκB selectivity of estrogen receptor ligands revealed by comparative crystallographic analyses. Nat. Chem. Biol. 2008, 4, 241−247. (36) Pike, A. C. W.; Brzozowski, A. M.; Hubbard, R. E.; Bonn, T.; Thorsell, A. G.; Engström, O.; Ljunggren, J.; Gustafsson, J. Å.; Carlquist, M. Structure of the ligand-binding domain of oestrogen receptor beta in the presence of a partial agonist and a full antagonist. EMBO J. 1999, 18, 4608−4618. (37) Shiau, A. K.; Barstad, D.; Loria, P. M.; Cheng, L.; Kushner, P. J.; Agard, D. A.; Greene, G. L. The structural basis of estrogen receptor/coactivator recognition and the antagonism of this interaction by tamoxifen. Cell 1998, 95, 927−937. (38) Hazra, A.; Gogtay, N. Biostatistics series module 6: Correlation and linear regression. Indian J. Dermatol. 2016, 61, 593−601. (39) Folmar, L. C.; Hemmer, M. J.; Denslow, N. D.; Kroll, K.; Chen, J.; Cheek, A.; Richman, H.; Meredith, H.; Grau, E. G. A comparison of the estrogenic potencies of estradiol, ethynylestradiol, diethylstilbestrol, nonylphenol and methoxychlor in vivo and in vitro. Aquat. Toxicol. 2002, 60, 101−110. (40) Blair, R. M.; Fang, H.; Branham, W. S.; Hass, B. S.; Dial, S. L.; Moland, C. L.; Tong, W.; Shi, L.; Perkins, R.; Sheehan, D. M. The estrogen receptor relative binding affinities of 188 natural and xenochemicals: Structural diversity of ligands. Toxicol. Sci. 2000, 54, 138−153. (41) CERI. Report on Evaluation and Method Development for Hormone-like Effects of Exogenous Substances. 2000 Contract Investigation/Research on Environment-Compatible Technology Development on Behalf of the Ministry of Environment and Industry, 2001. (42) Sheeler, C. Q.; Dudley, M. W.; Khan, S. A. Environmental estrogens induce transcriptionally active estrogen receptor dimers in Yeast: Activity potentiated by the coactivator RIP140. Environ. Health Perspect. 2000, 108, 97−103. (43) OECD. Second Meeting of the Validation Management Group on Screening and Testing for Endocrine Disrupters (Mammalian Effects), 2000; pp 20−21. (44) Kojetin, D. J.; Burris, T. P. Small molecule modulation of nuclear receptor conformational dynamics: implications for function and drug discovery. Mol. Pharmacol. 2013, 83, 1−8. (45) Morris, G. M.; Goodsell, D. S.; Halliday, R. S.; Huey, R.; Hart, W. E.; Belew, R. K.; Olson, A. J. Automated docking using a Lamarckian genetic algorithm and an empirical binding free energy function. J. Comput. Chem. 1998, 19, 1639−1662. (46) Wilk, M. B.; Gnanadesikan, R. Probability Plotting Methods for the Analysis of Data. Biometrika 1968, 55, 1−17. (47) Welch, B. L. The Generalization of ’Student’s’ Problem When Several Different Population Varlances Are Involved. Biometrika 1947, 34, 28−35.

6629


Importance of Receptor Conformations in Docking Calculation-Based

Recommend Documents