Molecular docking and 3D-QSAR studies on the glucocorticoid receptor antagonistic activity of hydroxylated polychlorinated biphenyls

Abstract The glucocorticoid receptor (GR) antagonistic activities of hydroxylated polychlorinated biphenyls (HO-PCBs) were recently characterised. To further explore the interactions between HO-PCBs and the GR, and to elucidate structural characteristics that influence the GR antagonistic activity of HO-PCBs, molecular docking and three-dimensional quantitative structure–activity relationship (3D-QSAR) studies were performed. Comparative molecular similarity indices analysis (CoMSIA) was performed using both ligand- and receptor-based alignment schemes. Results generated from the receptor-based model were found to be more satisfactory, with q2 of 0.632 and r2 of 0.931 compared with those from the ligand-based model. Some internal validation strategies (e.g. cross-validation analysis, bootstrapping analysis and Y-randomisation) and an external validation method were used respectively to further assess the stability and predictive ability of the derived model. Graphical interpretation of the model provided some insights into the structural features that affected the GR antagonistic activity of HO-PCBs. Molecular docking studies revealed that some key residues were critical for ligand–receptor interactions by forming hydrogen bonds (Glu540) and hydrophobic interactions with ligands (Ile539, Val543 and Trp577). Although CoMSIA sometimes depends on the alignment of the molecules, the information provided is beneficial for predicting the GR antagonistic activities of HO-PCB homologues and is helpful for understanding the binding mechanisms of HO-PCBs to GR.


Introduction
Polychlorinated biphenyls (PCBs) are highly persistent synthetic chemicals and resistant to biodegradation. Even though the production of PCBs was prohibited in many countries beginning in the late 1970s, these compounds are still one of the most ubiquitous and ABSTRACT The glucocorticoid receptor (GR) antagonistic activities of hydroxylated polychlorinated biphenyls (HO-PCBs) were recently characterised. To further explore the interactions between HO-PCBs and the GR, and to elucidate structural characteristics that influence the GR antagonistic activity of HO-PCBs, molecular docking and three-dimensional quantitative structure-activity relationship (3D-QSAR) studies were performed. Comparative molecular similarity indices analysis (CoMSIA) was performed using both ligand-and receptor-based alignment schemes. Results generated from the receptor-based model were found to be more satisfactory, with q 2 of 0.632 and r 2 of 0.931 compared with those from the ligand-based model. Some internal validation strategies (e.g. cross-validation analysis, bootstrapping analysis and Y-randomisation) and an external validation method were used respectively to further assess the stability and predictive ability of the derived model. Graphical interpretation of the model provided some insights into the structural features that affected the GR antagonistic activity of HO-PCBs. Molecular docking studies revealed that some key residues were critical for ligand-receptor interactions by forming hydrogen bonds (Glu540) and hydrophobic interactions with ligands (Ile539, Val543 and Trp577). Although CoMSIA sometimes depends on the alignment of the molecules, the information provided is beneficial for predicting the GR antagonistic activities of HO-PCB homologues and is helpful for understanding the binding mechanisms of HO-PCBs to GR. persistent groups of endocrine disrupting chemicals (EDCs) in the environment worldwide [1].
When PCBs enter organisms, cytochromes P450 (CYPs) oxidise them to hydroxylated PCBs (HO-PCBs) [2]. Although HO-PCBs are readily excreted, these compounds can be bound to blood proteins such as transthyretin, which facilitates the passing of HO-PCBs to blood in the umbilical cord [3]. Over the past few years, several studies have reported that HO-PCBs are present in the blood and placental tissue of humans and wildlife. Gomara et al. [4] firstly reported OH-PCB concentrations in placenta samples (53-261 pg/g wet weight) from a population in Madrid, Spain. Bytingsvik et al. [5] estimated the plasma concentrations of OH-PCBs in polar bear (Ursus maritimus) mothers (80 ± 38 ng/g wet weight) and their 4-month old cubs (49 ± 21 ng/g). The results indicated that exposure to OH-PCBs was a potential health risk for polar bears and the levels of OH-PCBs were above levels associated with health effects in humans and wildlife.
There is increasing evidence that HO-PCBs have endocrine disrupting potencies by interfering with various hormone receptors such as estrogen receptor (ER) and thyroid receptor (TR) [6][7][8]. Recently, HO-PCBs have been demonstrated to have glucocorticoid receptor (GR) activity [9]. GR is a member of the steroid-thyroid-retinoid receptor superfamily. By binding with glucocorticoid, the GR can block or modulate the undesirable effects of glucocorticoid excess, and affects immune response, development and metabolism in target tissues [10]. Takeuchi et al. [9] performed a systematic in vitro reporter gene assay to characterise the GR agonistic and antagonistic activities of HO-PCBs. To show the ability of the presented method in identifying positive compounds from negative compounds, a total 100 of the available OH-PCBs were tested for GR agonistic/antagonistic activities. It was noteworthy that 30 OH-PCBs acted as GR antagonists, while none of them exhibited agonistic activity.
As yet, the knowledge about the toxicological mechanisms of HO-PCBs with various hormone receptors is relatively limited. Moreover, the number of OH-PCBs is huge (e.g. there is a total of 837 monohydroxylated PCBs) [11] and detailed toxicological information for each compound is difficult to obtain. Therefore, there is increasing interest in presenting the complex biological interactions between ligands and nuclear receptors by simplified computational methods (in silico) known as quantitative structure-activity relationship (QSAR) studies [12]. Some QSAR studies have been performed on HO-PCBs with estrogenicity [13], thyroid hormonal activity [11] and androgenicity [14]. However, there is a lack of QSAR studies on the GR-antagonistic activity of HO-PCBs. Furthermore, the information about the binding of HO-PCBs with GR is limited. As a useful method to improve our knowledge of the protein-ligand interactions at the atomic level, molecular docking has become an integral part of many structure-based computational simulations [15].
In the present study, molecular docking and three-dimensional QSAR (3D-QSAR) were adopted to explore the structural characteristics which affect the GR antagonistic activity of HO-PCBs. Comparative molecular similarity indices analysis (CoMSIA) was employed to generate a 3D-QSAR model and to reveal the regions where interactive fields may affect the activity. Two different alignment rules (i.e. ligand-and receptor-based alignments) were applied to derive the rational conformations for developing the 3D-QSAR model. The robustness and predictive ability were tested by several internal and external validation strategies. Molecular docking was applied to study the interactions between HO-PCBs and GR. The results could provide some useful information on the structural characteristics that affect the GR antagonistic activity of HO-PCBs and elucidate the protein-ligand interaction mechanism.

Data sets
The data on GR antagonistic activity for 30 HO-PCBs were taken from the literature [9] and,were expressed as 20% relative inhibitory concentration (RIC 20 , M). In this study, the antagonistic activities, used as dependent variables in the 3D-QSAR analyses, were expressed as pRIC 20 (-log RIC 20 ). The whole data set was divided with a ratio of 4:1 into a training set (containing 24 compounds) for 3D-QSAR model generation and a test set (containing 6 compounds) for model validation, respectively. The compounds of the test set were selected randomly based on the fact that they could appropriately represent the structural diversity as well as the distribution of activities. All the compounds and associated GR-antagonistic activities are listed in Table 1. To further evaluate the feasibility and predictive capacity of the generated Experimental activities were derived from the literature [9]. c the predicted activities were calculated from the comSiA model.

No.
Compound a Cl sites pRIC 20  3D-QSAR model, the GR antagonistic activities of eight additional HO-PCBs, which were not included in the data set, were predicted by the model obtained.
The 3D-structures of HO-PCBs were generated using the SYBYL 7.3 molecular modelling package (Tripos Inc., St Louis, MO, uSA). Energy minimisation was performed using the Tripos force field [16] with a distance-dependent dielectric function and the Powell conjugate gradient algorithm. Convergence criterion of 0.001 kcal/mol·Å and a maximum of 1000 iterations were used for the energy-minimising process. MMFF94 charges were assigned to each compound. The energy-minimised structure was used for molecular docking (allowed to be flexible) and 3D-QSAR studies.

Molecular docking
The molecular docking and 3D-QSAR study were performed using SYBYL. The Surflex-dock program [17] interfaced with SYBYL was employed for docking studies to explore the probable binding conformation of HO-PCBs in the binding pocket of the target protein. The crystal structure of human GR with antagonist conformation (PDB code: 3H52) used in the molecular docking was retrieved from the Research Collaboratory for Structural Bioinformatics (RCSB) Protein Data Bank [18]. Prior to docking, the original ligand and the water molecules of GR were removed, polar hydrogen atoms were added to protein, and Kollman-all charges were assigned to protein atoms. In this paper, the automatic-based docking mode was applied with parameters in default values. During the docking process, the ligand molecules were considered to be flexible and the receptor was regarded as being rigid. The active site was defined as a cube covered the binding site with an edge length of 12 Å and centered at the midpoint of the longest atom-atom distance in the respective co-crystallised ligand in the protein. All molecules in the data set were docked into the active site of the receptor, and 10 conformations were obtained for each ligand molecule. The conformation with the top score of each ligand was selected as the most likely bioactive conformation [19] and used directly for the following CoMSIA research.
To verify the accuracy of molecular docking, the original ligand was generated by SYBYL and re-docked into the binding pocket of the receptor [20]. The docking procedure was reliable if the root mean square deviation (RMSD) value between the best-scored conformation determined by docking and the experimental one was less than 2.0 Å [21,22].

Molecular alignment
Molecular alignment of compounds is the key procedure in the development of 3D-QSAR models [23]. In order to derive the optimal 3D-QSAR statistical model, two different alignment rules were adopted: ligand-based and receptor-based alignments. In the former approach, all compounds were aligned to the most potent compound 2'-OH-CB58 by the Align Database command in SYBYL. During the latter process, the bioactive conformations of all compounds derived from molecular docking were assigned MMFF94 charges, and aligned automatically together and used directly for CoMSIA research.

CoMSIA analysis
The CoMSIA models were developed for the aligned molecular data set. In CoMSIA, five molecular fields including steric, electrostatic, hydrophobic, hydrogen bond donor, and hydrogen bond acceptor fields were calculated. For the generation of CoMSIA molecular fields, a 3D cubic lattice was created with a grid spacing of 2 Å and extending 4 Å units beyond the aligned molecules in all directions. The default value of 0.3 was used as the attenuation factor (α).
The 3D-QSAR models were derived using partial least squares (PLS) regression, which was implemented in SYBYL. Five CoMSIA descriptors (i.e. energies of steric, electrostatic, hydrophobic interactions, and hydrogen bond donor/acceptor) were used as independent variables and pRIC 20 as the dependent variable. In the PLS regression, a leave-one-out (LOO) cross-validation analysis was performed to determine the optimum number of components (OnC) and cross-validation correlation coefficient (q 2 ). q 2 is a statistical parameter used to evaluate the predictive ability of the model and it is calculated by the following equation: where n is the number of compounds, and Y predicted , Y observed and Y mean are the predicted, observed and mean values of the target property (pRIC 20 ), respectively.
The OnC was adopted in the non-cross-validation analysis and the final PLS regression model for CoMSIA was generated. The non-cross-validation results were evaluated by several statistical parameters such as coefficient of determination (r 2 ), standard error of estimate (SEE) and F value. To speed up the PLS analysis and reduce the noise, column filtering was used at the default value of 2.0 kcal/mol. The CoMSIA model was considered reliable and acceptable if q 2 was greater than 0.50 and simultaneously r 2 was greater than 0.90 [24,25]. To refine the obtained models, region focusing was performed on conventional CoMSIA models [26,27], which involved giving additional weight to the lattice points in a given CoMSIA region to increase the contribution of those points in a further analysis. In region focusing, discriminant power values were used as weights with different weighing factors applied in addition to grid spacing to obtain more predictive models [28].

Model validation
Some internal validation strategies such as the leave-many-out (LMO) cross-validation analysis (leave 10% out and leave 20% out, each of 25 runs) [29], the bootstrapping analysis [30] for 100 runs and the Y-randomisation for 10 times [31] were used as rigorous tests to further assess the robustness and statistical significance of the derived model.
The external validation method, the most valuable validation method [32][33][34], was applied to assess the predictive ability of the obtained model. The overall predictive ability of the CoMSIA model was externally validated by predicting the activity of compounds in the independent test set (not included in the original training set) [35,36]. The predictive ability of the model was expressed by the predicted r 2 (r 2 pred ), which was calculated using Equation (2): (2) r 2 pred = 1 −

PRESS SD
where SD is the sum of squared deviations between the experimental activities of the compounds in the test set and the average activity of the compounds in the training set, and PRESS is the sum of squared deviations between the experimental and predicted activity values.

Docking studies
Before the docking analysis, the original ligand generated by SYBYL was re-docked into the binding pocket of the GR. A low RMSD 0.7107 Å between the top-scored conformation of original ligand and the experimental conformation elucidated that the docking parameters specified in the Surflex-dock program were reliable and can well recover the interaction modes for the observed natural GR complex. All compounds in the data set were then docked into the active site of the GR in the same way. To elucidate the interactions, attention was focused on the protein-ligand interactions of the most potent compound, 2'-OH-CB58. The top-scored conformation of 2'-OH-CB58 in the active site of the GR is displayed in Figure 1. It can be seen that 2'-OH-CB58 was anchored at the binding site via the hydrogen bond and hydrophobic interactions. Specifically, the hydroxyl group of 2'-OH-CB58 formed a hydrogen bond with the carbonyl oxygen atom of Glu540 (OH···O=C, 2.3 Å, 110.8°). Additionally, a benzene ring of 2'-OH-CB58 was accommodated in the hydrophobic pocket that was composed of residues Ile539, Val543, Ala573 and Trp577. These aforementioned residues (Glu540, Ile539, Val543, Ala573 and Trp577) were important in the binding pocket and were the main contributors to the ligand and receptor interactions.

CoMSIA statistical results
The statistical results of the PLS analysis using two alignment approaches (ligand-based and receptor-based) are summarised in Table 2. The original ligand and receptor based CoMSIA models were unsatisfied due to a low r 2 value below 0.5. Subsequently, these less significant models were optimised through the region focusing method with a discriminant power value of 1.0 and a grid spacing of 1.0 Å. The statistical results in Table 2 showed that the receptor-based CoMSIA model with region focusing gave a significant result. In CoMSIA analysis, bioactive conformations are desired for molecular alignment. Different alignment rules resulted in different CoMSIA models. Receptor-based alignment led to models with better statistics than those from the ligand-based approach, presumably because the alignment using conformation derived from molecular docking is more realistic [37,38]. Therefore, our attention focused mainly on this receptor-based model for further discussion. using five components, this model yields a cross-validated q 2 of 0.632, a non-cross-validated r 2 of 0.931, a SEE of 0.040 and an F value of 63.819, indicating a strong strength for the correlation between the experimental and predicted pRIC 20 values. The relative field contributions of the steric, electrostatic, hydrophobic, hydrogen bond donor and hydrogen bond acceptor fields were 13.3, 23.7, 25.6, 19.2 and 18.2%, respectively. This suggests that all the fields make considerable contributions, and that electrostatic and hydrophobic interactions play relatively more important roles in the GR antagonistic activity of HO-PCBs.

Analysis of model validation
In addition to the LOO cross-validation, the LMO cross-validation analysis using 10 (leave 10% out) and 5 (leave 20% out) groups repeated 25 times each was performed to investigate the stability of the CoMSIA model. The mean q 2 values of leave 10% out cross-validation and leave 20% out cross-validation were 0.619 and 0.635, respectively. The q 2 values of the LMO analysis were comparable with those of the LOO analysis, further confirming that the obtained CoMSIA model was stable and possessed reasonable internal predictive ability. The robustness and stability of the model were also verified by the bootstrapping validation method for 100 runs. The bootstrapping method produces new data sets by randomly choosing samples from the original data sets. This method yielded a r 2 boot value of 0.952 and a SEE boot value of 0.025, suggesting good internal consistency in the underlying data set. All these results supported the robustness and the statistical validity of the developed model. To eliminate the possibility of chance correlation and further test the robustness of the CoMSIA model, the Y vector (pRIC 20 values) was shuffled 10 times. The obtained q 2 and r 2 values were in the range of -0.368 to 0.135 and 0.118 to 0.467, respectively, indicating the results from the CoMSIA model were not caused by chance correlation.
External validation was also carried out. The reliability and the predictive ability of the model were determined by a set of six test compounds not included in the model generation. The predicted r 2 pred of 0.879 was achieved, suggesting that the model possessed good external predictive ability and high accommodating capacity. The predicted activities calculated using the CoMSIA model and the residual values of the training set and the test set are listed in Table 1. A graphic representation of the experimental versus predicted antagonistic activity for both the training and the test set is shown in Figure 2. Table SI1 in the supplementary material which is available via the multimedia link on the online article webpage shows the predicted pRIC20 of eight additional HO-PCBs using the obtained CoMSIA model. 6-OH-CB35, 2'-OH-CB68 and 4'-OH-CB79, all with no experimental antagonistic activity in reporter gene assay [9], had a lower pRIC20 relative to those compounds with antagonistic activity (Table 1), indicating a good predictive capacity of the developed CoMSIA model. In addition, 4'-OH-CB14, 6'-OH-CB101, 3'-OH-CB138 and 3'-OH-CB180 showed a considerable predicted pRIC20 (Table SI1). Therefore these compounds may have antagonistic activity and should be further tested by experimental assay.

CoMSIA contour maps
To further identify the important structural regions of the compounds which could explain differences in the biological activity, 3D colour-coded field contour maps were generated. The coefficients of contour maps were generated using the field type 'Stdev*Coeff' (the standard deviation and the coefficient) with default values of 80% favourable and 20% unfavourable contributions. The contour maps of CoMSIA fields superimposed on those of the most potent compound 2'-OH-CB58 (Figure 3).
In the CoMSIA steric contour map ( Figure 3A), the steric interactions are represented by green and yellow contours, where bulky groups near the green regions increase the antagonistic activity but become sterically unfavourable near the yellow regions. A large yellow contour covers the 2′ and 6′ positions of the benzene ring, implying their preference for the sterically moderate and less crowded substituents here. This can be supported by the fact that 6-OH-CB36 was more potent than 6'-OH-CB26. The activity discrepancies between 3'-OH-CB53 and 4'-OH-CB50 can also be explained by this contour.
The electrostatic contour map is depicted in Figure 3B, where the electrostatic interactions are represented by blue and red contours. The blue regions, which represent the electropositive groups near these regions, are favourable to the activity and the red regions indicate that the electronegative groups close to these regions may increase the activity. A blue region covers the 6′ position of benzene ring, indicating that the electropositive groups are favourable for enhancing activity and that the electronegative group in these positions will decrease the activity. 4'-OH-CB33 exhibited higher potency than its counterpart 4'-OH-CB30, verifying this conclusion. In addition, 2'-OH-CB58 was more potential than its counterpart 2'-OH-CB121.
The hydrophobic contour map of the CoMSIA model is presented in Figure 3C, where the yellow region was hydrophobically favourable for the activity, while the white region was hydrophilically favourable. A large white contour covers the 2′ and 3 positions of the benzene ring, implying that hydrophobic groups were favored at these positions for increasing the potency. That 6-OH-CB36 had higher potency than 6′-OH-CB26, 6-OH-CB70 was more potential than 2′-OH-CB61, and 6′-OH-CB112 had higher potency than 2′-OH-CB121 verified this conclusion. This result was consistent with the observations obtained from previous molecular docking analysis, as the 2′ position of the benzene ring was pointing toward the hydrophobic pocket, interacting with the Ile539, Val543, Ala573 and Trp577 regions. Additionally, a medium-sized white contour covering the 5 position indicates that occupancy by a hydrophilic group in this region would promote antagonistic activity. This may explain why 4'-OH-CB106 exhibits higher activity than 6-OH-CB83.
The CoMSIA hydrogen bond donor and hydrogen bond acceptor contour maps are shown in Figure 3D and Figure 3E, respectively. The hydrogen bond donor field is represented by cyan and purple contours, in which cyan contours indicate regions where the hydrogen bond donor group would be beneficial to the activity, while the purple contours represent areas where the hydrogen bond donor group would decrease the activity. The hydrogen bond acceptor field is shown in magenta and red contours, in which the magenta contours demonstrate regions where a hydrogen bond acceptor group would improve the activity, whereas the red contour reveal regions where the unfavourable acceptor group would decrease the activity. near 2 position of the benzene ring, a purple region and a red region were located. Both hydrogen bond donor and acceptor contours emerge in the same region, revealing that a hydroxyl group may take part in hydrogen bond interaction with the receptor. Therefore, the presence of a hydroxyl group as hydrogen bond donor or acceptor substituent may form hydrogen bonds with residues and influence the potential activity. At this position, careful selection of groups with electron-withdrawing or electron-donating ability therefore is required [39].

Conclusions
In this work, 3D-QSAR and molecular docking were performed to build models for the GR antagonist activity of HO-PCBs. The conformations of all compounds, obtained from energetic minimisation and molecular docking, were used for the CoMSIA calculations. The obtained optimal receptor-based CoMSIA model exhibited good robustness and predictive ability according to several statistical parameters. Based on the detailed graphical interpretation of contour maps, the derived model helped to point out some significant structural features that influenced the GR antagonistic activity of HO-PCBs. Docking studies revealed that the hydrogen bond interaction with Glu540 and hydrophobic interactions with residues in the hydrophobic pocket played important roles in the binding of ligand with receptor. The results obtained from 3D-QSAR and molecular docking are able to accurately predict the GR antagonistic activity of HO-PCBs, identify structural features influencing the GR antagonistic activity, and provide a better understanding of the interaction between HO-PCBs and the GR.

Disclosure statement
no potential conflict of interest was reported by the authors.

Funding
This work was supported by the national Science and Technology Support Program (no. 2012BAC02B00).