Predicting The Toxicity of Emerging Pollutants via Molecular Docking: Interaction of Polyhalogenated Carbazoles with Human Androgen Receptor

Polyhalogenated carbazoles (PHCZs) are persistent chemical pollutants increasingly detected in different environmental matrices. They are structurally and chemically similar to other persistent organic pollutants (POPs), which are capable of disrupting the endocrine systems. However, PHCZs such as 3-chloro-9H-carbazole (3-CCZ), 3-bromo-9H-carbazole (3-BCZ), and 36-dibromo-9H-carbazole (36-BCZ) are rarely discussed in the context of their adverse effects on human health. Using molecular docking to investigate the potential toxicity of these PHCZs with the human androgen receptor (AR), this study finds that 36-dibromo-9H-carbazole and 3-bromo-9H-carbazole are potential AR antagonists, with the former being more toxic than the latter. This finding is on account of the presence of both Asn705 and Thr877 in the hydrophobic interaction of 36-BCZ, while only Thr877 is found in the hydrophobic interaction of 3-BCZ. Hence, PHCZs with higher bromine substitutions are more likely to be endocrine disruptors. Moreover, their binding sites with the human androgen receptor are similar to that of the androgen (agonist). Therefore, this study suggests that PHCZs may readily penetrate and disrupt the human androgen receptor (AR), providing the groundwork for future research studies and experimental validation on the molecular docking employed.


INTRODUCTION
Polyhalogenated carbazoles (PHCZs) are persistent and bio-accumulative emerging pollutants which occur in various environmental media (Yue et al., 2020). Their structural similarity to other well-known pollutants such as polychlorinated biphenyls (PCBs), polychlorinated dibenzofurans (PCDFs), and dioxins reflects their potential adverse effects on the environment and to human health, particularly in endocrine systems (Yue et al., 2020). The United States Environmental Protection Agency (USEPA) reported that endocrine-disrupting chemicals (EDCs) alter both the endocrine systems and the homeostasis systems, which are indispensable in the optimal functioning of living organisms (Lauretta et al., 2019;Kiyama et al., 2015). These EDCs are grouped according to their origin, but PHCZs are not formally included because of the limited mechanisms of their toxic effects as potential EDCs as well as their controversial distribution globally. Despite the lack of information on the global distributions, occurrence, and sources of these organic pollutants, the presence of the PHCZ pollutants in the environment continues to rise because of industrial wastes and wastewater discharges (Lauretta et al., 2019).
Common research targets of numerous studies are 3-chloro-9H-carbazole (3-CCZ), 3-bromo-9Hcarbazole (3-BCZ), and 3,6-dibromo-9H-carbazole (36-BCZ) due to their wide detection in the environment. For instance, the study by Mumbo et al. (2017) showed that these three PHCZ compounds, along with 3,6-dichloro-9H-carbazole (36-CCZ) slowly degrade in the aquatic environment as their maximum wavelength absorptions are in the near-ultraviolet region (295,296,299, and 301 nm) (Kiyama et al., 2015;Yan et al., 2018). Moreover, 3-CCZ is reported to be among the persistent PHCZs in soil samples due to its halogenation effects. It stays in the soil for more than 120 days, exceeding the European Chemicals Agency (ECHA) benchmark guidelines (Kiyama et al., 2015). A number of published articles relate the high structural similarity of these three PHCZs and dioxin-like compounds (DLCs) as indication that the former can induce the adverse physiological effects of the latter (Ma et al., 2019). For instance, 3-BCZ and 36-BCZ are two of the nine PHCZs that can be docked in the same binding site of dioxins (Ma et al., 2019). Ji et al. (2019) further reported that both 3-CCZ and 3-BCZ have the same level of potency as that of the dioxins.
Meanwhile, the human androgen receptor (AR) serves as the target receptor for this study. Human AR regulates the actions of androgens, which have important role in the development and function of the reproductive, musculoskeletal, cardiovascular, immune, and neural systems (Davey and Grossmann, 2016). It is a soluble protein from the nuclear receptor (NR) super family, which is responsible for the different regulations of the physiological functions in the body (Xu et al., 2013). It binds to androgens (agonist) with strong affinity to initiate the male reproductive development and growth (Tan et al., 2015). However, once the AR malfunctions upon binding with an antagonist, androgen insensitivity syndrome (AIS) occurs (Tan et al., 2015). AIS is a disorder in the external and internal genitalia in which a person with this condition has sexual development abnormality. Other related endocrine disorders include decreased sperm count, infertility and diabetes mellitus (Wahl and Smieško, 2018).
Several predictive models for ecological risk assessments involves comprehensive studies on the mechanistic interactions of these pollutants to their respective target sites (Escher and Hermens, 2002). Practical use of this approach can offer a possibility to reduce animal studies in the laboratory and eliminate tedious experiments that are often contentious and costly (Villeneuve and Garcia-Reyero, 2011). With the advent and development of several computational programs, researchers can veer away from the relatively costly use of the conventional testing paradigm. Hence, the use of such programs can address the complicated risk analysis of pollutants that usually comes with copious toxicological data. The National Research Council proposed computational methods such as molecular docking to aid in evaluation and assessment of toxic compounds (Huang and Zou, 2010). Molecular docking is used to predict the binding modes and affinities of ligands to its target macromolecule or protein (Huang and Zou, 2010). It serves as an efficient technique to predict the effects of the interactions between contaminants and associated target receptors (Walker and McEldowney, 2013) by identifying the most favorable binding mode(s) of the ligand to its target (Schleinkofer et al., 2006). With a relatively easy and robust computational approach, such as molecular docking, the binding modes and affinities of the PHCZs to the AR are predicted (Escher and Hermens, 2002). Direct binding at the target binding pockets is an essential mechanism to determine whether emerging pollutants such as PHCZs are possible endocrine-disruptors or not. Molecular docking analysis can therefore guide the researchers in evaluating the different chemical toxicity endpoints of the environmental Volume 33, Number 2, July-December 2022 • KIMIKA contaminants, typically employed in cell-based bioassays, reporter gene assays, in vivo/ex vivo evaluation assays and omic-based investigations.
The complexity of the ecological and toxicological risks of PHCZs led the present study to investigate the different interactions of the 3 PHCZs, such as 3-chloro-9H-carbazole (3-CCZ), 3bromo-9H-carbazole(3-BCZ), and 3,6 -dibromo-9H-carbazole(3,6-BCZ) to the human Androgen Receptor (AR) via molecular docking. PHCZs which have origins from anthropogenic byproducts are known for their persistent bioaccumulation in the environment. This explains why PHCZs have received much research attention. However, despite the publication of several studies providing evidence of their occurrence and distributions in the environment (Guo et al., 2017), their detrimental effects on living organisms are yet to be explored extensively, thus limiting the availability of the information pertinent to the mechanism of action of PHCZs (Yue et al., 2020). To date, there is no comprehensive study focused on the mechanistic, biochemical interactions of PHCZs with AR yet. The present work aims to contribute to the empirical research on the PHCZs by: 1) identifying the main interactions between different classes of PHCZs and AR, 2) studying the differences in the mode of interactions of these compounds with the protein receptor, and 3) deducing general insights relevant to the assessment of the potential binding sites for interactions with AR as influenced by the type and nature of halogen substitutions in PHCZs. Molecular models depicting the relationship of the chemical structures of these PHCZs with their respective binding interactions with AR were also explored. However, for this study, there was no attempt to characterize the absorption, distribution, metabolism, and/or excretion of the PHCZs in silico as these are intended to be determined experimentally via molecular (DNA) testing, as part of the future validation studies. The main computational software/applications used to carry out these modeling studies were Autodock Tools, IBM SPSS Statistics 25, and Biovia Discovery Studio 2020.

METHODS
All software programs and web servers used in this study are accessible and can be downloaded for free. Most of the parameters employed are adapted from different related literatures (Villeneuve and Garcia-Reyero, 2011;Kamerlin et al., 2020;Huang and Zhou, 2010;Walker and McEldowney, 2013;Schleinkofer et al., 2006;Guo et al., 2017;Singh et al., 2017), but some modifications were made depending on which suit best with the objectives of this study and on what can be accommodated by the computational programs employed in this study.

Materials and Software.
The crystal form of AR in PDB format was downloaded at the Protein Data Bank (PDB code: 2PNU), while the 3D conformers of 3-CCZ, 3-BCZ, and 36-BCZ in SDF were obtained from PubChem. The following software programs/web server were utilized to ensemble and characterize the docked conformations of each ligand; Swiss PDB Viewer, Avogadro, Pymol, AutodockTools, Biovia Discovery Studio, IBM SPSS Statistics 25.

Pocket Detection.
A pocket detection server CASTp 3.0 (Computed Atlas of Surface Topography of proteins) was employed to visualize and list the receptor's (PDB ID: 2PNU) pocket sites, which contain the participating amino acid residues (Villeneuve and Garcia-Reyero, 2011). This webserver uses the solvent accessible, molecular surface model and a standard default probe radius of 1.4̇ to compute the volumes and areas of these pockets (Tian et al., 2018). It was able to identify a total of 27 pockets, which represent the regions in the AR that can be accessible to different ligands and substrates for binding. However, in this study, only the first three pockets with the largest volumes were considered in this study, since they are most likely the active sites of the protein to which the ligand can interact with (Kamerlin et al., 2020). Table 1 provides the list of the participating amino acid residues in these pockets alongside their respective exact volumes and areas. For clarification, active sites are defined as the cavities in the protein which serve as potential binding sites for the ligands as these could contain the most number of amino acid residues which can possibly interact with the AR. As shown in Table 1, Pocket 1 is considered as the most active site of the human AR as it exhibits the largest volume and surface area, which constitutes most of the amino acid residues which can be crucial to the binding. According to literatures, it also contains the functionally important residues such as Asn705 and Thr877 in the androgen binding sites from which agonists and antagonists significantly interact with (Tamura et al.,2008). Macromolecule Preparation. The ligands that form complex with the AR were removed in Pymol, while Swiss PDB viewer was used to optimize the AR. Energy minimization computations was done in vacuo GROMOS96 43B1. This forcefield repairs the distorted geometries of the AR to achieve its most stable conformation since its initial structure was intended for its original bound ligands. An energy minimized report was then displayed upon completion of the optimization. The optimized AR was read in AutodockTools (ADT) using its graphical user interface (GUI). The water molecules were removed, while hydrogens and Gasteiger charges were added and assigned in the AR, respectively. The resulting file was saved in new PDB format for later docking.

Ligand Preparation.
The SDF format of the ligands was converted in PDB file using Open Babel. Geometry optimization of the molecular structure of the ligands was carried out using Avogadro. Since the ligands are organic compounds, MMFF94 was selected as the forcefield extension (Cerqueira et al., 2016), while the operating algorithm used was the steepest descent to reduce the net forces on the atoms (Jaidhan et al., 2014). All calculations were set to default as this is a typical workflow for small input ligands such as the PHCZs (Sikakwe et al., 2020). Energy minimization was achieved when dE has reached zero (Cerqueira et al., 2016). Then, each of the optimized ligands was converted to PDBQT format using Open Babel and was saved for later docking.

Molecular Blind Docking.
AutodockTools program was used to dock the PHCZs and the AR. The docking's working environment was created by copying the autodock4 and autogrid4 extension files from the Scripps Research Institute folder. The prepared files of the ligand and the receptor were placed in the same working environment. Then, the link address of the folder was copied and pasted into the startup directory of ADT. The PDB file receptor was loaded, followed by the ligand. Blind docking was carried out by adjusting the grid box at the center of the ligand, covering the entire area of the macromolecule thereby allowing the ligand to rotate freely. The Autogrid4 extension file was run together with the grid parameter file (GPF) to generate the receptor map files of the atoms. For the docking parameter file (DPF), a genetic algorithm (GA) was utilized as the search parameter to generate docked conformations (Wang et al., 2003). Seventy (70) independent dockings per ligand were carried out to produce ten (10) Table 2 summarizes the input values needed for GA depending on the number of torsions of the ligand. To save the docking parameter file (DPF), Lamarckian GA (4.2) was selected for the sampling complex conformation formed by the ligand and AR. The autodock4 extension file and the DPF were run to generate the docking log file (DLG).

Table 2. General guidelines for Energy Evaluations of Rigid Ligands and Proteins
Number of Torsions GA_num_evals GA_num_generations *0 25000 to 250 000 27 000 1-10 250 000 to 25 000 000 27 000 >10 >25 000 000 27 000 *torsion defined for the ligands in this study Autodock then provides the docking scores in kcal/mol, which served as the predicted binding energies of the complex. It can be found in the second column of the DLG file once the program has finished the computation. Autodock uses a systematic search method that explores the ligands' binding modes and accounts for any energy changes during docking or complex formation (Pantsar and Poso, 2018). Intermolecular energies such as Van der Waals, hydrogen bonding, desolvation, and electrostatic energies all contributed to the total docking score of each docked conformation (Schleinkofer et al., 2006 ). The docking results were presented in a histogram using Origin® Pro to visualize the three most populated clusters per PHCZ, that is, those containing the highest possible number of conformations (Xu et al., 2013).

Hierarchical Clustering Analysis (HCA) and Principal Component Analysis (PCA).
The docking scores of the 10 x 3 bound ligand poses per target were used as the weighing factor for principal component analysis (PCA) and hierarchical cluster analysis (HCA). HCA clusters each Result No. based on the similarities of docking scores by applying agglomerative schedule hierarchical clustering using Ward's method and Euclidean distance measure (Sikakwe et al., 2020). Ward's method is commonly used for quantitative variables, which analyzes the variance of the clusters by merging smaller clusters but are further apart. To simplify the large set of the data, only the main cluster with the smallest level of the Euclidean distance formed is selected, and its corresponding subclusters are subjected to further analysis. The HCA of this study did not consider the subclusters formed at a much lower distance if they belong to a higher cluster. Meanwhile, PCA used a correlation matrix to reduce the large dataset of the docking scores clustered in each Result No. by converting them into principal components (PC). Only the principal component (PC) with eigenvalues greater than 1 was retained as it contains the most number of important information in the dataset (Sikakwe et al., 2020). Both analyses were performed using IBM SPSS Statistic 25. The conformation with the lowest docking score contained in the significant Result Nos. from HCA and PCA was loaded in Biovia Discovery Studio 2020 to enable the characterization of the ligand-AR interactions.
Statistical Analysis. The significant differences among the binding energies and RMSD values of the 3 PHCZs were tested using one-way Analysis of Variance (ANOVA) in IBM SPSS Statistics 25 software. The test is said to be significant if p-value < 0.05. Validation Protocol. The known AR antagonists, such as prochloraz and vinclozolin, were individually docked with the AR. The results were then compared with results from published studies (Singh et al., 2017;Tian et al., 2018) to validate the present work. Meanwhile, docking of the PHCZs with androgen, a known agonist, was also performed. This validation was done following the foregoing procedures. Meanwhile, details on the overall molecular docking analysis employed in this study is shown in Figure 1.

RESULTS AND DISCUSSION
Validation Results. To test the robustness of the docking protocol employed in this study, exogenous ligands, such as vinclozolin and prochloraz, were individually docked with human androgen receptor (AR) for 70 independent runs. The docking results in Figure 2 show that prochloraz produced a relatively higher number of conformations compared with vinclozolin, attributable to the presence of more rotatable bonds on the structure of prochloraz, as illustrated in Figure 3. This finding is consistent with the readings of Autodock that prochloraz and vinclozolin consist of 8 and 2 rotatable bonds, respectively. According to Schleinkofer et al. (2006), searching for the complexes' possible conformations increases exponentially due to the highly branched torsions or more rotatable bonds in a ligand. Hence, the findings provide plausible evidence that docking results are time-variant, suggesting that extended calculation is observed as the input ligand becomes larger. This distinct variation is consistent with the rigid body docking principle, that is, the search for search for conformational space for complexes whose ligands with ≤ 6 rotatable bonds, is significantly faster without compromising accuracy, thereby, producing fewer conformations per cluster.
Further, the docking protocol employed in this study was able to establish hydrophobic interactions between AR and two amino acid residues, namely, Thr877 and Asn705 (Table 3). Such interactions demonstrated the ability of these ligands as AR antagonists. With reference to the findings of the Comparative Molecular Field Analysis (CoMFA), the structural diversity of AR antagonists is addressed by relating their mechanistic behavior through the Near 10 ̇ Polar Interaction rule (Tamura et al., 2008). This rule hypothetically states that the distance for the interactions of the antagonists to the receptor should be less than 10 Ȧ to prevent any H-bond interactions with Thr877 and Asn705 (Tamura et al., 2008). The observed potential binding interactions are in agreement with this rule as well as with the published reports, suggesting formation of predominantly hydrophobic interactions. Accordingly, this imply that both    Docking of PHCZs to the Human AR. Foremost, the computational approach employed in this study considers the protein receptor as rigid. As such, this study only considered the flexibility of the PHCZs by taking on different possible conformations as each binds with the protein receptor. Blind docking was done since the potential binding sites for the PHCZs in the human AR were not predetermined and that PHCZs are not originally bound in the human AR, therefore, their target pockets and potential binding sites are also not known. Blind docking will maximize the search for the binding sites within the receptor. If the PHCZs were made to bind with just a specific predetermined binding sites as indicated with previous studies, this study will be limited to assume that it will be the best binding site and the only binding site possible, neglecting that, binding could be influenced by possible ligand-specific (in this case PHCZs)-effects, which this study primarily sought. Thus, blind docking was performed using a docking box to determine the potential binding sites of the PHCZs. The coordinates of the docking box are not limited to certain KIMIKA • Volume 33, Number 2, July-December 2022 values as long as it is larger than the default docking box, which is enough to cover the entire AR allowing the PHCZs to move freely (Ghersi and Sanchez, 2009). In every docking, the parameters of the docking box vary to maximize the search for possible binding sites. Hence, the conformation or the orientation of each PHCZs also varies for each generated binding site. For instance, the PHCZs were individually blind docked with AR, and results show that 3-BCZ contains the highest number of conformation ( Figure 4). The possible conformations of 3-BCZ in each Result No. are relatively similar to that of the 3-CCZ and 36-BCZ. This finding is due to the highly rigid structure of all these PHCZs, with their rings completely fused, as shown in Figure 5. However, the complex with the highest number of conformations does not automatically imply that it has the strongest affinity, among others. This is because not all the predicted conformations of the PHCZs are totally bound in the AR when loaded in Autodock. Similarly, complexes with a lower number of conformations do not directly signify weaker binding. To address these docking inconsistencies, a total evaluation of the 10 x 3 bound PHCZs ligand poses are considered in the subsequent discussions.  The hierarchical clustering analysis (HCA) performed on the 10 x 3 bound PHCZ complexes is presented in three dendrograms, which display the number of clusters formed. Each cluster consists of Result Nos that are distinctively grouped into separate subclusters based on the similarity of their docking scores. The vertical and horizontal axes of the dendrograms represent the Euclidean distances measurement and the Result Nos from blind docking, respectively. The height of the vertical branches in dendrograms indicates the degree of similarity. The smaller the level of the Euclidean distance (vertical axis of dendrogram) formed, the greater the similarity within the groups merged in a particular cluster (Sikakwe et al., 2020). Furthermore, this study interprets HCA by focusing on the main clusters, which showed the highest degree of similarity. Subclusters fused at a lower level of Euclidean distance belonging to a higher main cluster were not considered. The statistical measures that provide support to the dendrogram analysis are also reported, which include mean, standard deviation and variance. These parameters are commonly used by statisticians to measure the variability within the given data set (Marsal, 1987). A value of zero variance indicates no variability in the data (R. Wilcox, 2003). Meanwhile, according to Ramirez and Cox (2012), the rule of thumb for the maximum acceptable range of the variance calculated is approximately four times the standard deviation. Hence, as a guide in the analysis, the data with the lowest variance was considered which was based on how close its value to zero. The calculated variance outside or greater than the range rule of thumb will be considered irrelevant.
HCA of 3-CCZ. Figure 6 shows the HCA results of 3-CCZ, which consisted of the three main clusters. cluster 1 contains the highest number of Result Nos such as 3, 9, 10, 5, 6 and 1, while cluster 2 comprises Results 4 and 8. In cluster 3, Results 2 and 7 were merged, both of which contain the least number of conformations with 26 and 29, respectively. It can be observed that clusters 1 and 2 are linked slightly higher than cluster 3, which makes the latter the most similar. Cluster 3 has binding energies closely clustered around the mean while those in clusters 1 and 2 were more dispersed over a wider range of values due to higher standard deviation and variance (Table 4). HCA of 3-BCZ.  HCA of 3-BCZ. Meanwhile, Figure 7 shows the HCA results of 3-BCZ, which comprised of four main clusters. Clusters 1 and 2 are equally formed at a lower Euclidean distance compared to clusters 3 and 4. Cluster 1 consists of Result Nos 5, 6, and 10, while cluster 2 has Result Nos 2 and 4. In cluster 3, the Result Nos formed are 1 and 7 and lastly, cluster 4 constitutes Result Nos 8, 9, and 3. Furthermore, the figure shows that cluster 2 has the largest similarity among the other main clusters since its Euclidean distance is relatively lower than clusters 3 and 4, and it rapidly connects its subclusters from the bottom of the dendrogram than cluster 1. Incidentally, Cluster 2 also contains the Result No. with the lowest standard deviation and variance (   Figure 8 displays three main clusters. Cluster 1 consists of Result Nos. 4, 10 and 3. Cluster 2 has the highest number of Result Nos formed such as 7, 9,5,6,2, and 8, while only Result No.1 comprised Cluster 3. Clusters 1 and 2 were fused at an equal Euclidean distance in which cluster 1 shows the closest link at the bottom of the dendrogram, which depicts greater similarity between its subclusters. Meanwhile, cluster 2 fused its subclusters at a much higher Euclidean distances and cluster 3 is clearly observed as being completely separated from all others as it took half of the dendrogram to be formed at a higher distance. The variance and standard deviation of cluster 3 indicate that its binding energies are totally spread out. Hence, the distribution of its docking scores is substantially different compared to other clusters. On one hand, cluster 1 was selected since it does not contain the Result No. with the highest SD and variance (Table 6).  Therefore, this study suggests that the number of conformations in the Result Nos does not significantly contribute to clustering as some vary greatly despite being fused in one main cluster. Moreover, it is important to note that no other parameters can elaborately explain the HCA of this study as the Result Nos were not measured at different scales. Rather, the statistical analysis of the Result Nos. simply aids the selection of the clusters for further analysis. All values of variance calculated were within the range rule of thumb. Hence, clustering the Result Nos is solely based on the homogeneity of dockings scores they contain. Overall, clusters 3, 2, and 1 of 3-CCZ, 3-BCZ, and 36-BCZ, respectively, showed the highest degree of similarity and were further characterized.
Principal Component Analysis (PCA). In this study, PCA was carried out to reduce the large dataset of docking scores into a smaller number of principal components (PCs). Based on the Kaiser Index (Mahmoudi et al., 2021), only the PCs with eigenvalues greater than or equal to 1 are retained as they describe the major trends in the data (Sikakwe et al., 2020;Le et al., 2020). PCA classifies which among the ten (10) Result Nos of each PHCZ is strongly correlated to the PC (Jaidhan et al., 2014). Table 7 reports the principal component (PC), eigenvalues, and percentage variance for each PHCZ. One PC was identified per PHCZ, which represents the total number of possible sources of variation in the data. In other words, PCs contain the greatest number of sources of information pertinent to docking scores. For instance, the PCs of 3-CCZ, 3-BCZ, and 36-BCZ explained 95%, 94% and 96% of the total variance in docking scores, respectively. The Result Nos 10, 5 and 9 of 3-CCZ, 3-BCZ and 36-BCZ, respectively correlated most strongly to each of the identified PC. This suggests that the possible conformations of each PHCZs will tend to have docking scores contained in the most strongly correlated Result No. Hence, only the docking scores were treated as describing variables that could influence the PC extracted.  Table 8 summarizes the Result Nos which showed the highest degree of similarity in HCA and strongest correlation to PC from PCA. Each value in column 3 of Table 8 represents the conformation with the lowest binding energy contained in each Result No, which is considered as the best-docked conformation. In this study, each conformation will be referred to as the Potential Binding Mode (PBM). Overall, PBM2, PBM5 and PBM3 of 3-CCZ, 3-BCZ and 36-BCZ, respectively showed the lowest binding energy among the PBMs of the PHCZs (Table 8, columns 2 and 3). This means that the PHCZs at these PBMs have a longer residence time or association with the AR's binding site, which may be long enough to elicit any effects. These findings are consistent with the generally accepted principle that low binding energies calculated, imply stronger and preferential affinity of the ligands with the target receptor (Sikakwe et al., 2020).  Table 8 were characterized using Biovia Discovery Studio 2020. It is important to note that the amino acids listed in Table 9 directly interacted with the N-H moiety of the PHCZs and not with another amino acid. Furthermore, as previously discussed, this study treated the AR as a rigid body; thus, the PHCZs were the only molecules expected to change its conformation, while finding its way to bind with the AR.

Characterization of the Potential Binding Modes (PBMs). The interaction mechanisms of the selected PBMs from
For instance, Figure 9a displays the interacting residues of the PBM2 of 3-CCZ, which include Gly708, Arg752,Phe876,Leu707,Met780,Leu873,and Met745. Specifically, the π-electron cloud of the phenyl ring of 3-CCZ formed π-alkyl and alkyl interactions with the electron group of Leu707 and Met745, while the lone pair of the chlorine atom formed similar interactions with KIMIKA • Volume 33, Number 2, July-December 2022 Phe876, Met780, and Leu873. The van der Waals interactions were also observed with Gly708 and Arg752. Furthermore, the N-H moiety of 3-CCZ and Gly708 may have acted as a hydrogen donor and formed H-bond with Leu704 with distance of ~2.00 Ȧ and ~1.87Ȧ, respectively.
In Figure 9b, nine interacting amino acid residues were observed from the docking of 3-BCZ with the AR, five of which, namely, Arg752, Leu707, Gly708, Met745, and Thr877 were involved in van der Waals interaction. The π-electron cloud of the phenyl and the five-membered N-containing ring of 3-BCZ interacted with Phe764 forming a somewhat T-shaped conformation. At the same time, its Br atom formed π-alkyl and alkyl interaction with Phe876 and Leu873. One hydrogen bond interaction was formed with the N-H moiety of 3-BCZ with Leu704 with H-bonding distance of 2.28 Å (Table 9).
Meanwhile, Figure 9c shows that Gly708 formed hydrogen bonding interaction with Leu704 (1.85 A ) ̇i n lieu of N-H moiety of 36-BCZ. The PBM of 36-BCZ further demonstrated van der Waals interaction with Arg752, Met745,Leu707,Gly708,Asn705,and Thr877. Other important interactions involved are π -sulfur bond with Met787, the π -alkyl and alkyl interaction with Met749, Val746, Phe876, Phe764, and π-π T-shaped interaction with Phe764. Overall, the interacting amino acid residues involved in the PBMs of the PHCZs are bound in the right-handed alpha-helix structure of the AR ( Table 9). The common interacting residues involved are Arg752, Gly708,Phe876,Leu707,Met745,and Leu704, which are all present in Pocket 1 as listed in Table 1. This implies that the three PHCZs can be docked in the first pocket of AR ( Figure  9d), which in this study was tagged as its most active site. It is also observed from Figure 9d that the PHCZs partly overlap with one another signifying a possible ligand-ligand interaction. At each docking, the orientation/conformation of the PHCZ may change as a result of energy minimization, a process which allows the formation of the most stable conformation for the ligand during the docking process. From the mechanistic standpoint, variations in the observed conformations of the PHCZs could possibly lead to the alteration of the protein structure, causing the protein to become dysfunctional (Ma et al., 2019). Dysfunctionality might then lead to negative effects, which may be related to the toxicity of the ligands under investigation. It is emphasized that there could be several binding sites in the protein with different amino acid residues interacting with the ligands. In fact, because of the complexity associated with the toxicity assessment, these binding sites have not all been established experimentally. Thus, the results presented in this study attempt to showcase possible interactions that may or may not lead to toxicity in vivo. For clarification on the toxicity evaluation and mechanism exploration, there is a possibility that the binding of any contaminant with the protein receptor doesn't take place in one site only and that binding may cause activation of the protein function, to some extent. But it does not follow that the activation could immediately indicate potential toxicity. In some case, the receptors can be mediators of ligand toxicity and this is what this study is trying to sort out, that is, to come up with a simulation and determine whether the observed interactions could be ligand-specific human AR variations.
The docking results obtained in this study revealed that the binding potency of the three PHCZs with the AR is in the order of: 36-BCZ (-7.62 kcal/mol) > 3-BCZ (-7.18 kcal/mol) > 3-CCZ (-7.08 kcal/mol). This relative potency is influenced by the varying number and type of halogen substitutions on the base structure of the PHCZs. The presence of two electronegative bromine atoms on the phenyl rings of 36-BCZ may have resulted to strong, preferential affinity with the AR. Meanwhile, although the Cl atom of 3-CCZ is more electronegative than the Br atom of 3-BCZ, the former PHCZ came last in terms of affinity, which may indicate that the large size (steric effect) of bromine may have induced hydrophobic interactions, demonstrating preferential interaction with bromine substituted PHCZ over the chlorine substituted PHCZ.
Risk Assessment of PHCZs in complex with the AR. The toxic mechanisms and potential effects of the PHCZs with the AR are best understood by evaluating the type and number of important amino acid residues involved in the interaction. Though the presence of hydrogen bonds and hydrophobic interactions are ubiquitous in any receptor-ligand binding however, these may not be completely disregarded in determining the efficiency of the effective interactions via simulation studies. Overall, 3-BCZ and 36-BCZ showed potential interaction with the AR, with the latter being the most pronounced. This observation could be attributed to the hydrophobic interaction of Thr877 and Asn705 with 36-BCZ. On the other hand, 3-BCZ only interacted with Thr877. As reported, the preferential interaction of these amino acid residues with PHCZs, could elicit adverse effects both on the standpoint of environment and human health, thereby suggesting the potential of PHCZs as endocrine disrupting compounds (Tian et al., 2018;Xu et al., 2013). Moreover, the potential ability of PHCZs as AR antagonists may be enhanced if substituted with bromine atoms as well as with increasing number of bromine atom substitutions.
To support the findings that the PHCZs under study are potential AR antagonists by hydrophobically interacting with Thr877 and Asn705, androgen (a known agonist) was docked KIMIKA • Volume 33, Number 2, July-December 2022 with the human AR following the same docking protocol discussed in the methodology section. The positive control, androgen, was also found to bind in Pocket 1 of the androgen receptor (AR), in almost similar manner as observed with the binding of PHCZs. Most of the interacting residues in the binding mode of androgen can also be found in the PHCZs (Figure 9). This suggests that PHCZs can immediately bind with the AR and might readily induce toxic effects such as inhibition, mimicking the response of the natural hormone, androgen.
This study did not cover omics-based investigation and/or in vivo or ex vivo assays, which would allow quite extensive experimental examination of the human AR-PHCZs interaction to indicate agonist/antagonistic activities or enhancement or inhibition of the ligand-dependent physiological function of the human AR. Human AR regardless whether bound to antagonist or to agonist, can elicit structural transformation which could lead to direct binding with functional components and performing a specific function or to interact with signaling pathway. For instance, in the study by Brinchmann et al., 2019, binding with contaminants such as polyaromatic hydrocarbons (PAHs) can trigger calcium responses in human endothelial cells via receptor-dependent nongenomic signaling, which are relevant to the toxicity of the said contaminant.
To highlight the importance of this study, docking results were correlated with the partition coefficient or the n-octanol/water partition ratio usually expressed as log ( ) (Hodges et al., 2019). The partition coefficient is a useful parameter in determining the environmental fate and biological effects of chemical pollutants. However, as this in silico study did not include the determination of pharmacokinetic parameters, the partition coefficient was not derived in this study. Instead, the potential ability of these PHCZs to diffuse through the cell membrane and target the androgen receptor was ascertained based on the partition coefficient of structurally related compounds such as polychlorinated biphenyls (PCBs) since both compounds are planar halogenated hydrocarbons (Pantsar and Poso, 2018). As per review of related literatures, PCBs with greater number of chlorine atoms have demonstrated higher tendency to diffuse through the cell membrane as they are readily soluble in fats (Mahmoudi et al., 2021). They alter and disrupt the membrane structure causing cell death (Gao et al., 2005), and can migrate to fatty tissues and bio-accumulate in the food chain (Singam et al., 2019). They also interfere with the transcriptional activation of the androgen receptor (Singam et al., 2019). The octanol-water partition coefficient [log ( )] of the PCBs is expected to increase with the increasing molecular weight of the PCBs (more chlorine atoms) (Liang et al., 1998).
In view of the foregoing discussion about the behavior of PCBs, the relative lipophilicity of the PHCZs may also be deemed influenced by the number of halogens in the structure. As such the relative lipophilicity of the PHCZs would most likely follow the order of: 36-BCZ (325 g/mol)>3-BCZ (246.1g/mol)> 3CCZ (201.65g/mol), with 36-BCZ having the most number of halogen substitutions. Moreover, since PHCZs are halogenated polycyclic aromatic hydrocarbons, most of them may also be classified as stable and generally have hydrophobic character with > 4 and stable (Pantsar and Poso, 2018). These literature-based inferences may shed light on the potential antagonistic behavior of 3-BCZ and 36-BCZ. However, it is emphasized that the results indicated herein are empirical models only and therefore the antagonistic activity of the PHCZs must be verified through experimental work such as molecular (DNA) testing or toxicity testing, to corroborate with the endocrine-disrupting ability of PHCZs as demonstrated by the docking results. As such, the calculation of the partition coefficient of the PHCZs to AR and other pharmacokinetic parameters will be determined during the experimental validation and shall constitute the next scope of work to be performed in line with this study.