Structural Characterization of the Binding Interactions of Various Endogenous Estrogen Metabolites with Human Estrogen Receptor α and β Subtypes: A Molecular Modeling Study

In the present study, we used the molecular docking approach to study the binding interactions of various derivatives of 17β-estradiol (E2) with human estrogen receptor (ER) α and β. First, we determined the suitability of the molecular docking method to correctly predict the binding modes and interactions of two representative agonists (E2 and diethylstilbesterol) in the ligand binding domain (LBD) of human ERα. We showed that the docked structures of E2 and diethylstilbesterol in the ERα LBD were almost exactly the same as the known crystal structures of ERα in complex with these two estrogens. Using the same docking approach, we then characterized the binding interactions of 27 structurally similar E2 derivatives with the LBDs of human ERα and ERβ. While the binding modes of these E2 derivatives are very similar to that of E2, there are distinct subtle differences, and these small differences contribute importantly to their differential binding affinities for ERs. In the case of A-ring estrogen derivatives, there is a strong inverse relationship between the length of the hydrogen bonds formed with ERs and their binding affinity. We found that a better correlation between the computed binding energy values and the experimentally determined logRBA values could be achieved for various A-ring derivatives by re-adjusting the relative weights of the van der Waals interaction energy and the Coulomb interaction energy in computing the overall binding energy values.


Introduction
The endogenous estrogens are a group of vitally important female sex hormones with diverse biological functions. Disruption of their actions contributes to the pathogenesis of a number of disease states in humans, including disruption of reproductive functions [1][2][3], development of breast cancer [4][5][6][7][8], and many other abnormal conditions [9][10]. Undoubtedly, it is of considerable scientific interest to develop the ability to predict the estrogenic potency and efficacy of a given chemical. Compared to the widely used methods such as receptor binding assay, receptor activity assay and crystallography, computational molecular modeling methods have the potential advantages of low cost, high speed, and high throughput. Progress in this area of research will help develop the ability to predict and/or identify new environmental estrogens among millions of natural and synthetic compounds that humans are potentially exposed to so that the general public can be informed and thus avoid unwitting exposure to these chemicals.
There are mainly two categories of computational methods commonly used in studying the ligand-receptor interactions. One is the more traditional QSAR (Quantitative Structure-Activity Relationship) approach, and the other is the molecular docking approach. The QSAR approach relies solely on the structures of ligands to predict the three-dimensional properties of the receptor binding site as well as the potential binding affinity of a given unknown ligand. In comparison, molecular docking analysis is based on the three-dimensional structures of both the receptor binding pocket and the ligand to predict its preferred binding conformations in the receptor's binding site. The relative binding affinity of a ligand can then be estimated by computing the binding energy values (DG or DE). These two approaches have their own advantages and disadvantages. Whereas the QSAR model usually can achieve a good prediction when the unknown test compounds are structurally very similar to the training set of compounds used, the molecular docking analysis theoretically is more versatile in predicting the binding conformations of structurally-diverse compounds.
In an earlier study [11], we have compared the relative binding affinity (RBA) of some fifty estrogen derivatives (mostly the endogenous metabolites of 17b-estradiol, E 2 ) for human ERa and ERb. These estrogen derivatives, while sharing the same core structure as E 2 with only one or two small functional group(s) added (their structures shown in Figure 1), have vastly different RBAs for human ERa and ERb, ranging from having a similar binding affinity as E 2 to having little or no binding affinity at all (see Table 1). In the present study, these E 2 derivatives were used as model compounds to study their binding interactions with human ERa and ERb using the molecular docking approach.
There are two main reasons that we chose to use these structurally similar estrogen derivatives as model compounds. First, although the crystallographic structure of human ERa in complex with E 2 is known, the interactions of the receptor with various E 2 metabolites are actually not known. Studying the interactions of these compounds with human ERs may shed light on how small modification in the E 2 structure can drastically disrupt its binding interactions with the receptors. These results, in Figure 1. The chemical structures of 17b-estradiol (E 2 ) and the 27 E 2 derivatives. The number of each carbon is labeled next to the atom in the E 2 structure. The names of the E 2 derivatives were shown in the upper left corner of each frame. The RBA values of the E 2 derivatives for ERa and ERb (data from Table 1 as % of RBA of E 2 ) were shown in the lower right corner of each frame. The numbers were rounded to the nearest integer due to space constraint. doi:10.1371/journal.pone.0074615.g001 turn, may help us better understand the detailed structural characteristics of the binding interaction of an estrogen with various amino acid residues in the LBDs of human ERs. Second, we have successfully used the QSAR approach in an earlier study [11] to accurately predict the binding affinity of these structurally similar estrogen derivatives for human ERs. In the present study, therefore, we sought to use a molecular docking approach to predict the binding affinities of these estrogen derivatives for human ERa and ERb, so that the predicted abilities of these two modeling approaches could be compared.
The research strategies used in this study are as follows: We first evaluated whether the molecular docking method can be used to reliably and correctly predict the binding interaction of a ligand with the ligand binding domain (LBD) of human ERa. Next, by using the validated molecular docking method, we then set to determine the binding interactions of a total of 27 representative E 2 metabolites/derivatives (structures shown in Figure 1) with the LBDs of both human ERa and ERb. Lastly, we also explored ways to achieve a better prediction of the binding affinity by modifying the calculation method for the binding energy (DE) values.

Building of the 3-D Structures of E 2 Derivatives
For the ligands used in this study, their chemical structures and RBA to ERs are shown in Figure 1. As mentioned above, these estrogens share the same core steroid structure as E 2 ; in fact, most of them are naturally occurring metabolites of E 2 formed in humans [12].
Building of the ligand structures as well as energy minimization were performed using the InsightII modeling program (Version 2005, Accelrys Inc. San Diego, CA) installed in the Red Hat Enterprise Linux WS4.0 (Red Hat Inc. Raleigh, NC) operating system on a Dell Precision 690 workstation. The 3-D structures of these E 2 derivatives were built with the Builder module in InsightII Hydrogen bond lengths were quantified by measuring distances between the hydrogen atoms of the 3-hydroxyl group of estrogen derivatives and the O e of ERa-E353 or ERb-E305, between the oxygen atoms of the 3-hydroxyl group of estrogen derivatives and the H g of ERa-R394 or ERb-R346 and between the hydrogen atoms of 17hydroxyl group of estrogen derivatives and N d of ERa-H524 or ERb-H475. For the D-ring derivatives, two hydrogen bond lengths were listed. The first is formed by hydrogen atoms of 17-hydroxyl groups and the second is formed by hydrogen atoms of 16-hydroxyl groups. Relative binding affinity (RBA) was also listed for comparison. doi:10.1371/journal.pone.0074615.t001 and minimized with the Discover module. Energy minimization was carried out with the Polak and Ribiere conjuate gradients method in the Discover module in InsightII until the final convergence criterion reached the 0.001 kcal/molÅ . Consistent valence force field (CVFF) was used for energy minimization.

Selection of the Human ERa and ERb 3-D Structures as Templates
The use of molecular docking approach to study the binding interaction of a ligand with its receptor protein is usually based on the known 3-D structure of the receptor protein. As listed in Table 2, presently about a dozen crystal structures of the ligand binding domains (LBDs) of the human ERa and ERb are available in the PDB database [13][14][15][16][17][18][19]22]. While some of the receptor proteins are in complex with agonists (such as E 2 and DES), others are in complex with antagonists (such as raloxifene and tamoxifen). We compared the ERa structures in complex with three ERa agonists, namely, E 2 (PDB codes 1ERE and 1A52), DES (PDB code 3ERD), and genistein (GEN) (PDB code 1X7R) by superimposing these structures on each other. The Root Mean Square Distance (RMSD), a parameter that is commonly used to reflect the average distance between the same atoms in different structures, was found to be only approximately 0.5 Å for these ERa structures when they were complexed with different agonists. By superimposing the ERb structures in complex with two different ERb agonists, E 2 (PDB codes 3OLS) and ERB041 (PDB code 1X7B), their RMSD was found to be 0.6 Å . This indicates that the receptor structures resolved by different research groups under rather different experimental conditions have very similar overall structures. Therefore, in this study, the x-ray crystal structure of human ERa in complex with E 2 (PDB code: 1ERE) [14] and the x-ray crystal structure of human ERb in complex with a synthetic ERb agonist ERB-041 (PDB code: 1X7B) [18] were used as templates for docking various estrogen derivatives. Hydrogen atoms were added to the proteins, and the structures were minimized to relax the strains imposed by the addition.

Molecular Docking and Simulation Studies
Molecular docking were performed using the InsightII modeling program (Version 2005, Accelrys Inc. San Diego, CA) installed in the Red Hat Enterprise Linux WS4.0 (Red Hat Inc. Raleigh, NC) operating system on a Dell Precision 690 workstation. The flexible docking procedures were carried out using the Simulated Annealing Docking method in the Affinity module of InsightII. The binding pocket in ERa was defined to include all residues within the 7-Å reach of the original ligand E 2 , which include M343, L346, A350, E353, L387, L391, R394, L402, F404, I424, H524, and L525.
A combination of Monte Carlo and Simulated Annealing methods was used to explore all possible conformations of a ligand molecule in the binding pocket. One hundred conformations were obtained and the one with the lowest potential energy was chosen for further minimization. The backbone of the protein was fixed during the docking procedures. Using similar protocols as Table 2. Summary of current 3D structures of ERs in complex with various ligands (listed according to the chronological order of the publications).

ER isoform Ligands PDB code References Comments
ERa-LBD E 2 1ERE [14] The first steroid receptor structure with an agonist bound inside the LBD.
ERa-LBD Raloxifene (RAL) 1ERR [14] The first ER structure with a SERM bound inside the LBD.
ERa-LBD DES/peptide 3ERD [16] The first ER structure with an agonist and a co-activator peptide bound.
ERa-LBD 4-hydroxyl tamoxifene (OHT) 3ERT [16] In this structure the ligand in the LBD induces a similar orientation of H12 as seen in 1ERR.
ERb-LBD RAL 1QKN [17] The first ERb structure with an antagonist bound inside the LBD.
ERb-LBD GEN 1X7J [19] The ERb structure with a partial agonist bound inside the LBD. GEN has 40-fold higher binding affinity for ERb than ERa.
ERb-LBD E 2 3OLS [22] An ERb structure with an agonist bound inside the LBD.
doi:10.1371/journal.pone.0074615.t002 Figure 2. The overlay of the ligand-binding domains (LBDs) of ERa and ERb. The protein structures were shown in cartoon and colored green and magenta for ERa and ERb, respectively. E 2 molecules were shown in stick and colored blue and red in ERa and ERb LBD, respectively. a-Helixes and b-sheets in the ER LBDs are labeled according to references [14,16]. Helix 2 structures are missing in both X-ray structures. doi:10.1371/journal.pone.0074615.g002 described above for ERa, flexible docking of various ligands into the LBD of human ERb was also carried out. Detailed docking parameters are as follows. When we set up the docking run, the ligand binding site was defined, and the ''solvation_grid'' was set to ''on''. Alpha carbons of the receptor were held fixed during the docking run. The Affinity docking method uses a combination of Monte Carlo (MC) and Simulated Annealing (SA) approaches. In the initial phase, an MC search generated 100 unique ligand binding orientations. The ''quar-tic_vdw_no_coul'' was selected as the ''nonbond_method'' because randomly placing the ligand in the binding pocket could potentially lead to very severe divergences in the Coulombic and VDW energies. The ''vdw cutoff'' was set to 5, the ''scale_vdw'' was set to 0.1, the ''scale_HB'' was set to 1, and the ''scale_tether'' was set to 1. In addition, the ''TempOrERange'' was set to 200, the ''MxrChange'' was set to 180, the ''EnerToler'' was set to 1e+6, the ''RMS_Toler'' was set to 1. Minimization step was set to 10,000.
After the initial structures had been generated by the initial MC phase, these structures were then read in as an ''.arc file'' for further refinement by SA. The SA phase consisted of 50 cycles of 100 fs molecular dynamics (MD) simulations. The temperature is re-scaled at each cycle, with the first cycle running at 500 K and the final cycle at 300 K. Following the SA phase, the structure is minimized by 10000 steps of conjugate gradients. The ''Cell_multipole'' was selected as the ''Nonbond_method''. The ''diel_const'' was set to 1, the ''scale_vdw'' was set to 1, the ''scale_coul'' was set to 1, the ''scale_HB'' was set to 0.1, and the ''scale_tether'' was set to 1. The ''diel_const'' was set to 1 and the ''dis_dep_dipole'' was set to off.

Calculation of Binding Energy Value
According to the InsightII program, the total interaction energy (DE binding ) between the receptor protein and the ligand is usually calculated using the Evaluate Intermolecular Energy function of the Superimposed structures of docking result and the crystal structure of ERa LBD in complex with E 2 . The known crystal structure of ERa LBD in complex with E 2 (PDB code: 3ERE) was colored in yellow with E 2 colored in green. The docked E 2 was colored in orange and ERa LBD was colored in magenta. The green dashes indicated the hydrogen bonds formed. All the structures were shown in ball and stick. The atoms involved with hydrogen bond formation were colored according to the atom type, i.e. white for hydrogen, red for oxygen and blue for nitrogen. Hydrogen atoms in other amino acids were not shown. doi:10.1371/journal.pone.0074615.g003  Docking module. The DE binding value includes two components, namely, the van der Waals (VDW) interaction energy (DE VDW ) and the Coulomb interaction energy (i.e., electrostatic interaction energy, DE Coulomb ). Hydrogen bond energy is included as part of the Coulomb interaction energy. The following equation 1 is commonly used to calculate the total interaction energy DE binding : Theoretically, DE binding can be used to reflect the experimentally determined logRBA, as expressed in equation 2 below: Here a higher value of DE would reflect a lower binding affinity. A higher degree of correlation between DE binding and logRBA would indicate that the computational docking model has a higher degree of accuracy to predict a compound's relative binding affinity for the receptor.

Validation of the Molecular Docking Method
Consistent with an earlier report [19], when the structures of the ERa and ERb LBDs are superimposed to each other (shown in Figure 2), the basic folds of their LBDs are very similar to each other, and the overall positions of the binding pockets, as well as the ligand conformation, are nearly identical. Based on this information, next we chose to use the ERa LBD as a representative example for validation of the computational methods used in this study.
The 3-D structure of the ligand binding domain (LBD) of human ERa in complex with DES, a non-steroidal estrogen, has been previously determined (PDB code: 3ERD). To determine whether the computational docking approach can produce the same binding mode for DES as observed in the crystallographic study, we performed a cross-docking validation by using the crystal structure of the human ERa LBD in complex with E 2 (PDB code: 1ERE) as a template receptor for the docking of DES. The docking model of the ERa-DES complex was then compared with the known crystal structure of the ERa-DES complex (PDB code: 3ERD) by superimposing these two structures.
As shown in Figure 3A, the docked structure of ERa-DES complex has a very similar overall structure as the known 3-D structure of this complex. The overall binding conformation of the modeled DES ligand is almost identical to its experimentally derived structure with the structures showing that the two aromatic rings of DES have the same steric orientation, and form four hydrogen bonds with residues E353, R394 and H524 in the binding pocket. These three amino acids have less than 0.4 Å RMSD between the docked structure and the crystal structure. In addition, several other important residues in the binding pocket, namely, A350, L346, M343, L525, I424, L402, L391, F404 and L346, form almost identical hydrophobic interactions with DES in both modeled and crystal structures.
Similarly, we used the known crystal structure of ERa LBD that was complexed with DES (PDB code: 3ERD) as a template and then docked E 2 into the LBD. The docking results were compared with the known crystal structure of the ERa-E 2 complex (PDB code: 1ERE). As shown in Figure 3B, the docking method produced a nearly identical structure of the ERa-E 2 complex compared to the known 3-D structure of this complex. Three identical hydrogen bonds are identified in both docked and known crystal structures, which are formed between the two hydroxyl groups of E 2 and the amino acid residues E353, R394 and H524 in the binding pocket. Also, all hydrophobic amino acids in the binding pocket have nearly the same orientation and positioning in these two structures. Together, these two examples show that the docking method used can correctly predict the binding mode of a ligand in the LBDs of human ERs with a high degree of reliability and accuracy. with ERa LBD determined by molecular docking. The green dashes indicate the hydrogen bonds formed. All the structures are shown in ball and stick. The amino acids were colored according to the atom type, i.e. green for carbon, red for oxygen, blue for nitrogen and white for hydrogen. Among the amino acids in the binding site, only E353, R394 and H524 are shown in this figure. E 2 was colored in white. The ligands were shown in the following different colors: in panel A, 6a-OH-E 2 (yellow), 6b-OH-E 2 (orange), 6-keto-E 2 (pink), 6-dehydro-E 2 (red), 7-dehydro-E 2 (magenta), 9(11)-dehydro-E 2 (light blue), 11a-OH-E 2 (purple) and 11 b-OH-E 2 (green); in panel B, E 1 (magenta) estriol (16a-OH-E 2 ) (yellow), 16b-OH-E 2 (orange), 16-keto-E 2 (pink), 17a-OH-E 2 (red), 15a-OH-E 3 (dark blue), 16a-OH-E 1 (light blue), 16-keto-E 1 (purple), 16a-OH-E 2 -17a (brown), 16b-OH-E 2 -17a (grey). doi:10.1371/journal.pone.0074615.g006

Prediction of the Binding Mode of E 2 Derivatives in the LBDs of Human ERa and ERb
By using the same docking approach as described above, next we docked a total of 27 E 2 derivatives into the LBDs of human ERa and ERb. These 27 analogs were selected from a pool of some 50 estrogen derivatives tested in our earlier study [11]. They were divided into three groups according to their structural relationship with E 2 (see Table 1): A-ring derivatives (9 compounds including E 2 itself), B/C-ring derivatives (8 compounds), and D-ring derivatives (10 compounds). These 27 ligands all share the same core structure as E 2 with only one or two small functional group(s) being added to the steroid core (see Figure 1). The purpose of the docking study was to determine how these subtle modifications in the ligand structures alter the binding conformations and binding energy values (which reflect the binding affinity) of these ligands.
A-ring derivatives. The overall 3-D binding modes of these A-ring derivatives were very similar to that of E 2 . Two hydrogen bonds are formed between their 3-hydroxyl groups and the amino acid residues E353 and R394 in ERa or E305 and R346 in ERb. A third hydrogen bond was formed between the 17-hydroxyl group of E 2 and ERa-H524 or ERb-H475. In the ligand-bound structures that were modeled, other hydrophobic amino acids in the binding pocket had a similar conformation to E 2 . The docking results of various A-ring derivatives of E 2 with ERa LBD are shown in Figure 4. For most of the A-ring derivatives, there were only subtle differences between their binding modes in comparison with the binding mode of E 2 . Apparently, these small differences in the binding modes contribute importantly to their differential binding affinities for the ERs. One of the most notable differences in the binding of these A-ring estrogen derivatives was that the hydrogen bonds formed with their 3-hydroxyl group were of slightly different length. It is hypothesized that the difference in the hydrogen bond distance, which would affect the hydrogen bond strength, is an important determinant of the binding affinity of a given ligand. To test this hypothesis, we computed for comparison the distances of the following three hydrogen bonds formed between the 3hydroxyl group of various estrogens: the hydrogen bond between the H atoms of the 3-hydroxyl group and the O e of ERa-E353 or ERb-E305, the hydrogen bond between the O atom of the 3hydroxyl group and the H g of ERa-R394 or ERb-R346, and the hydrogen bond between the H atom of the 17-hydroxyl group of estrogens and the N d of ERa-H524 or ERb-H475. The values are summarized in Table 1. To decipher the contribution of each hydrogen bond to the overall binding affinity, we determined the degree of correlation between the hydrogen bond lengths with the experimentally determined logRBA values. Given that the hydrogen bond interactions are generally thought to be electrostatic interactions in nature and that the electrostatic potential is usually in an inverse first-order relationship with the distance between the two atoms, the inverse first-order curve regression was Figure 8. Relationship between the correlation coefficient r and x/y for the estrogen derivatives in Table 1. For each x/y value, the total binding energy was calculated according to equation (3) and the correlation coefficient r value is calculated by correlating the total binding energy for each chemical with its logRBA. The x/y value is shown in log scale. doi:10.1371/journal.pone.0074615.g008 thus used here to correlate the experimentally determined logRBA values with the hydrogen bond length.
For the A-ring derivatives, two hydrogen bonds are formed between the 3-hydroxyl groups of these estrogen derivatives and the amino acid residues E353 and R394 in ERa or E305 and R346 in ERb. As shown in Figure 5, for ERa, the inverse correlation for the distance of the hydrogen bond with R394 is better than that for E353, suggesting that the hydrogen bond formed between the 3-hydroxyl group of various estrogens and R394 plays a more important role in the binding interaction than the hydrogen bond formed with E353. For ERb, the correlation for the length of the hydrogen bond formed with E305 is stronger than that for R346, again suggesting that these two hydrogen bonds contribute differentially to the binding affinity. The substantial correlation between the experimentally determined logRBA values and the hydrogen bond lengths indicates the vital importance of the hydrogen bonds with the 3-hydroxyl group of E 2 in determining their overall binding affinity of a ligand. This correlation also suggests that computing the lengths of the hydrogen bonds formed between the 3-hydroxyl group of a steroidal estrogen and ERs can be used as an important predictor of the binding affinity of a given A-ring derivative of E 2 . The hydrogen bonds formed between the 17hydroxyl group of A-ring derivatives of E 2 with H524 in ERa or H475 in ERb were of similar lengths and less important to determine their binding affinities.
Earlier studies suggested that 4-OH-E 2 has a markedly lower dissociation constant for the ERs compared to E 2 , which may mean that 4-OH-E 2 can remain in the ERa LBD longer [12,20]. Our computational docking model showed that 4-OH-E 2 could form an ideal hydrogen bond between its 3-hydroxyl group and E353 and R394 in the same way as seen with E 2 . In addition, 4-OH-E 2 can also form a hydrogen bond between its 4-hydroxyl group and ERa-L387. In the case of ERb, a similar hydroxyl group is also formed with L339. This interaction is believed to account for its slower dissociation with the ERs. In this context, it is also of note that our docking results show that 2-OH-E 2 can also form an additional hydrogen bond between its 2-hydroxyl group and ERa-E353 or ERb-E305. Interestingly, despite the presence of additional hydrogen bond, 2-OH-E 2 actually has a lower binding affinity for both ERa and ERb compared to E 2 . Our docking model of 2-OH-E 2 provides unique insights. The presence of an additional 2-hydroxyl group causes a unfavorable interaction with the hydrogen bond formation between its 3-hydroxyl group and ERa-E353 or ERb-E305, which likely is the main cause that reduces its overall binding affinity for the ERs.
Among the A-ring derivatives tested, 2-MeO-E 2 and 4-MeO-E 2 caused the most drastic shift of the ligand position relative to E 2 due to their bulkier methoxy substituents, which shifts the ligand position and results in less favorable van der Waals interactions and also interferes with the formation of hydrogen bonds between their 3hydroxyl groups and ERa-E353 and R394 or ERb-E305 and R346, ultimately resulting in drastically reduced binding affinities.
B/C-ring derivatives. The docking results of various B/Cring derivatives with the LBDs of ERa and ERb are shown in Figure 6. Similar to A-ring derivatives, we noticed that their binding modes are very close to that of E 2 . Two hydrogen bonds are formed between the 3-hydroxyl group and residues E353 and R394 in ERa or E305 and R346 in ERb. The lengths of their hydrogen bonds are similar to those formed with E 2 . No hydrogen bonds were observed with either the 6-and 11-hydroxyl groups of the B/C-ring derivatives in the modeled structures.
For B/C-ring derivatives, the correlation coefficient r values between their hydrogen bond lengths and logRBA values are 0.1, 0.2 and 0.4 for the hydrogen bonds with E353, R394 and H524 of ERa, respectively ( Figure 5). For ERb, the r values are 0.3, 0.1 and 0.4 for the hydrogen bonds with E305, R346 and H475, respectively. The poor correlations for the B/C-ring derivatives indicate that the binding interactions contributed by the 3-and 17b-hydrogen bonds are not the determining forces of the binding affinity of the B/C-ring derivatives. As discussed below, the repulsive forces play a dominant role in determining the binding affinity of B-ring derivatives.
D-ring derivatives. For some of the D-ring derivatives with a hydroxyl group at the C-16 position, one more hydrogen bond can be formed with H524, which may contribute to their relatively higher binding affinity for the ERs as compared to other estrogen metabolites ( Figure 6). The correlation coefficient r 2 values for the D-ring derivatives are 0.13, 0.06 and 0.04 for the hydrogen bonds observed with E353, R394 and H524 of ERa, respectively (hydrogen bond lengths were shown in Table 1; correlation data were not shown). For ERb, the r 2 values are 0.00 (no correlation at all), 0.00 and 0.08 for the hydrogen bonds with E305, R346 and H475, respectively. For the D-ring derivatives that form two hydrogen bonds with histidine, only the shorter hydrogen bond was used in the correlation. Similarly to the B/C ring derivatives, the poor correlation for the D-ring derivatives suggested that the binding interactions contributed by the 3 and 17b hydrogen bonds are not the determining force of the binding affinity of the D-ring derivatives. The role of repulsive forces in determining the binding affinity of D-ring derivatives is discussed later.

Prediction of Binding Affinity of Estrogen Derivatives with ERa and ERb LBDs Based on Binding Energy Calculation
According to the docking models developed in this study for each of the 27 E 2 derivatives, we have also computed their binding energy (DE) values in complex with ERa or ERb (listed in Table 1). The correlation between the computed DE binding values and the experimentally determined logRBA values for various E 2 derivatives is shown in Figure 7. For both ERa and ERb, the correlation is highest for A-ring derivatives, with r values 0.87 and 0.78, respectively. In contrast, there was no meaningful correlation seen between the computed values and experimental values for the B/C-ring derivatives or D-ring derivatives.
The poor correlations between the calculated binding energy values and the experimental values for B/C-ring derivatives or D-ring derivatives are, in part, due to the inherent problems associated with the currently-used molecular docking method [21], namely, the relative inaccuracy of the molecular mechanics approach for computing the interaction energy levels. Because of the hydrophilic substitutions (i.e., the hydroxyl group) in their C-6 or C-11 position, both of which are surrounded by hydrophobic amino acid residues, these hydroxyl groups are expected to reduce the chances for their binding interactions with the receptors because of the presence of strong repulsive forces. This low correlation also reveals a major deficiency of the currently used computational program, which cannot properly weigh the role of the repulsive forces that may completely diminish the chances for effective binding interactions.
In an attempt to improve the modeling program to better access the repulsive force, we sought to modify the relative weight of VDW interaction energy (DE VDW ) and Coulomb interaction energy (DE Coulomb ) in computing the final interaction energy (DE binding ), by converting the commonly used equation 1 to the following equation 3: where x and y are assigned parameters to adjust the relative weights of DE VDW and DE Coulomb in the total binding energy. Using the above equation, by altering the x/y ratios, we can compute the corresponding r value for the relationship between the computed DE binding and experimentally-determined logRBA. The relationships between x/y value and r value for A-ring, B/C-ring and D-ring derivatives are shown in Figure 8. The optimal r value was achieved with different x/y values for Aring and B/C-ring derivatives. For A-ring derivatives, in the case of ERa, when x/y is 2, the optimal r value is 0.9. For ERb, the curve is basically flat, which means that VDW and Coulomb interaction energy are equally important to determine the RBA of a ligand. For B/C-ring derivatives with ERb, the optimal r is 0.8 when x is 0, which indicates that Coulomb energy is more important in determining RBA. Interestingly, when x/y is 1, the correlation is the worst. However, for B/Cand D-ring derivatives, the prediction is still very poor for both ERs. It is apparent that the inability of the docking program to accurately calculate the binding energy value for the B/Cand D-ring derivatives is not simply due to the inappropriate balance of the electrostatic interaction and VDW interaction. Further studies are warranted to better resolve this problem.
Conclusions and practical implicationsThe present study sought to use the molecular docking approach to characterize the interaction of E 2 derivatives with human ERa and ERb. First, we tested the suitability of the molecular docking method for predicting the correct binding mode of a ligand inside the binding pockets of human ERa and ERb. Using DES and E 2 as examples, we demonstrated that the docked structures are almost exactly the same as the known crystal structures of ERa LBD in complex with either of these two ligands. Using the same docking approach, we also docked 27 structurally-similar E 2 derivatives into the LBDs of human ERa and ERb. While their binding modes are very similar to that of E 2 , there are notable subtle differences. These small differences are believed to contribute importantly to their different ER binding affinity. In the case of the A-ring estrogen derivatives analyzed, there is a strong inverse relationship between the length of the hydrogen bonds formed with ERs and the binding affinity. In the study, we optimized our approach's ability to predict the relative binding affinity of a given A-ring derivative of E 2 through re-adjusting the relative contribution of the van der Waals (VDW) interaction energy and Coulomb interaction energy in computing the overall binding energy (DE) value.
However, for B/C and D-ring estrogen derivatives, the correlations between calculated binding energy and experimentally determined binding affinity is very low and re-adjustment of the relative contribution of the van der Waals (VDW) interaction energy and the Coulomb interaction energy did not improve their correlations. Notably, some recent studies have also noted that while the docking programs can be successfully used to predict the ligand poses, they appear to lack the ability to accurately predict the ligand binding affinity [23,24]. The inability of the docking programs to predict the binding affinity likely is mainly due to their lack of adequate consideration of the repulsive forces between the hydrophilic substitutions at the C-6 and C-11 positions with the hydrophobic residues around these regions in the docking approach. In addition, it is also rather problematic to account for the desolvation of the hydrophilic groups on the ER ligands as well as in the receptor binding pocket, which may also be a partial contributing factor. Methodological improvements in these particular areas may greatly enhance the docking and scoring approaches, and ultimately, the accuracy in calculating the binding energy values between a ligand and its binding protein.

Author Contributions
Conceived and designed the experiments: PW BTZ. Performed the experiments: PW BTZ. Analyzed the data: PW BTZ. Contributed reagents/materials/analysis tools: PW CM BTZ. Wrote the paper: PW CM BTZ.