A detailed theoretical exploration on the THR-β binding affinities and antioxidant activity of some halogenated bisphenols

Abstract Natural halogenated phenolic compounds are unique bioactive structures which share features and physicochemical properties with thyroid hormones, who are essential regulators of neurological development and metabolism processes. Also, these structures can be used as natural antioxidants to minimize the diseases related to oxidative stress. In this work, the binding affinity of 32 natural and synthetic halogenated bisphenols were investigated on thyroid hormone receptor-β (THR-β) using the molecular docking, MM/GBSA, molecular dynamics, and a rigorous three-layer ONIOM ((M06-2X/6-31G*:PM6:AMBER) calculation. The desirable potency is observed for binding of selected compounds to THR-β. The Met313, Asn331, and His435 are the main interacting residues in the binding cavity which involved in the hydrogen and halogen bond interactions with the ligands. The most potent candidate on binding to the active site of THR-β is presented with respect to the results of mentioned calculations. Moreover, the antioxidant activity of compounds has been investigated using the quantum mechanical calculations. The electrostatic potential surfaces illustrate well the antioxidant capacity of compounds. The halogen substituents increase the antioxidant activity of the most stable conformers. The position and number of OH groups are crucial factors which affect the activity, whereas two adjacent hydroxyl groups enhance the antioxidant activity of selected compounds. Communicated by Ramaswamy H. Sarma


Introduction
Several natural halogenated phenolic compounds were extracted from red, brown, and green algae, where the bromophenols were identified as the most potent bioactive compound by several studies (Flodin & Whitfield, 1999;Gordon, 1999;Liu et al., 2011;Xu et al., 2004) . They exhibited a wide spectrum of pharmacological activities (Akbaba et al., 2013;Balaydı n et al., 2010;Duan et al., 2007;Kurata et al., 1997;€ Oztaskı n et al., 2017;Shi et al., 2009;Wang et al., 2005;Xu et al., 2003). Several studies reported the diseases related to oxidative stress such as cancer, coronary heart disease, hypertension, obesity, diabetes, cataract, and inflammatory diseases are best minimized by the use of natural antioxidants (Ahmad et al., 2014;Ortega-Dom ınguez et al., 2017;Ovalle-Magallanes et al., 2017;Zhong et al., 2015). More than 30 bromophenols from marine algae were reported to have antioxidant activity using in vitro methods (Duan et al., 2007;Hye et al., 2007;Ke et al., 2008;Wanyi et al., 2010). Although there are some discussions on the structure-activity relationship (SAR) of these compounds, recent studies revealed halogenated phenolic compounds as potential candidates in the prevention of diseases related to free radical attack (Ming et al., 2011). Some halogenated phenolic compounds from algae such as halogenated bisphenols structurally resemble the thyroid hormones (Erin et al., 2018;Vanessa et al., 2012). Hence, they have potential of binding to different nuclear hormone receptors and may affect a variety of physiological functions (Thangaraj et al., 2019). Thyroid hormones are essential regulators of numerous physiological processes, from neurological development to metabolism (Mondal & Mugesh, 2017). They are produced from thyroid gland and released in the blood stream containing variable number of iodine substituents in their structures, as shown in Scheme 1 (Thangaraj et al., 2019). Thyroid hormones exert their genomic actions by binding to the two homologous nuclear receptors, thyroid hormone receptors a (THR-a) and b (THR-b), which are members of the nuclear receptor superfamily of ligand-activated transcription factors (Shulin et al., 2013;Thangaraj et al., 2019). THRs are important targets in drug discovery and environmental toxicity researches (Fei et al., 2010;Matthieu et al., 2003).
The positive metabolic effects of thyroid hormones are mediated by THR-b isoform, the predominant liver, kidney, and brain THR, while THR-a is found mostly in the heart, brain, and skeletal muscles (Romero et al., 2020). Both THRs consist of three domains including an N-terminal domain (NTD), a DNA binding domain (DBD), and a C-terminal ligand binding domain (LBD) (Bleicher et al., 2008;Joharapurkar et al., 2012). The ligand binding pocket of THRs closely matches the thyroid hormones, so they have high specificity and selectivity against the THRs. Previous studies revealed that the compounds which deviate from the molecular and physicochemical properties of thyroid hormones are excluded from binding site (Akinori et al., 2003;Ishihara et al., 2003). However, halogenated phenolic compounds which share features and physicochemical properties with thyroid hormones may bind to THRs. Several studies examined the affinity of halogenated compounds such as polybrominated diphenyl ethers and tetrabromobisphenol to THRs for mammals, and suggested a strong preference for a hydroxylated diphenyl core structure similar to that of thyroid hormones (Kitamura et al., 2008;Leonetti et al., 2016;Ren et al., 2013;Shigeyuki et al., 2005;Stapleton et al., 2012), while the position of hydroxyl groups and halogenation pattern affect the binding affinities (Kollitz et al., 2018).
Due to the biological importance of bisphenols reported in the last decade, the current study aims to investigate the antioxidant capacity of some naturally occurring halogenated bisphenols and their synthetic derivatives. Moreover, their binding affinity against THR-b were considered using the computational approaches.
Molecular docking predicts the affinity of compounds to bind a biological target using the scoring functions (Meng et al., 2011;Stapleton et al., 2012). It determines the mode of action of compounds in the LBD of THR-b through the various interactions. The active compounds were also identified by the molecular mechanics/generalized born surface area (MM/GBSA) calculations which use the molecular mechanics to calculate the binding free energies (Genheden & Ryde, 2015). In addition, the molecular dynamics (MD) simulations can provide a worthwhile insight into the activity of compounds by probing their dynamics and binding mechanisms. They are powerful tools for understanding the driving forces underlying molecular recognition, accelerating drug discovery, and guiding molecular design (Buch et al., 2011;Faghih et al., 2015;Qi et al., 2019). Finally, the exact binding features and interaction details were investigated using theoretical methods based on the ONIOM (our own Nlayered integrated molecular orbital and molecular mechanics) methodology (Vreven et al., 2006). This is an important technique to study very large systems by the division of an entire system into several regions that can be studied at different theoretical levels (Chung et al., 2015). Therefore, the most interesting portion of system is allowed to consider via a high resolution quantum mechanics computations, which is not essential for the entire system.

Computational methods
Thirty-two halogenated bisphenol derivatives used in this work are shown in Scheme 2. The most stable conformers were obtained via two-dimensional relaxed scan around C1-C2 and C1-C2 0 rotatable bonds (see Scheme 2), from 0 to 180 by 10 intervals using the B3LYP method (Becke, 1992;Lee et al., 1988) in conjunction with the 6-31 G(d,p) basis set.
All calculations were performed using the Gaussian 09 program package . The global minima on the potential energy surfaces (PES) were optimized at the B3LYP/6-311þþG(d,p) level before further calculations. In addition to the most stable neutral conformers, the related anions, radicals, and cation radicals were also fully optimized at the above mentioned level of theory. The frequency calculations were performed to estimate the thermodynamic properties, and characterize the stationary points, which verified those as the ground state structures with no imaginary frequencies.
Three common chemical pathways for the scavenging of free radicals by antioxidants are proton coupled-electron transfer (PC-ET/HAT), sequential electron transfer proton transfer (SETPT), and sequential proton loss electron transfer (SPLET) (Leopoldini et al., 2011;Senthil Kumar & Kumaresan, 2012;Wright et al., 2001). The reaction environment can affect the balance between these mechanisms (Olszowy, 2019). These are quantified by some antioxidant descriptors, such as the bond dissociation enthalpy (BDE), the proton affinity (PA), the electron transfer enthalpy (ETE), the ionization energy (IE), and the proton dissociation enthalpy (PDE), which are defined in Equations (1)-(5) for the ArOH antioxidant  The X-ray crystal structure of THR-b in complex with the thyroid hormone T 3 was extracted from the RCSB protein data bank (PDB) with entry 3GWS code (Nascimento et al., 2006). The optimized structures were applied to molecular docking using the molecular operating environment (MOE) (Software Available from Chemical Computing Group Inc., 1010 Sherbrooke Street West, Suite 910, Montreal, Canada H3A 2R7) and the Maestro Schrodinger (Maestro, 2017) program packages.
The correction of the structural issues of crystallographic data, the optimization of the hydrogen bond network (protonate 3D), and the energy minimization of crystal structures were performed using QuickPrep application in the MOE program package (Software Available from Chemical Computing Group Inc., 1010 Sherbrooke Street West, Suite 910, Montreal, Canada H3A 2R7), from the chemical computing group, prior to docking. The MOE-Alpha site finder was used to generate the binding site. The compounds were docked into ligand binding cavity of THR-b using MOE-Dock, with the Triangle Matcher method used as the search protocol and the AMBER99 molecular mechanics force field for interactions, and further refined by the London dG rescoring function. The best pose of each ligand-receptor complex corresponds to the highest negative energy scoring in kcal mol À1 (Behazin & Ebrahimi, 2018). The interactions of ligands with the key residues take also into account to choose the best pose of each ligand.
In addition, the Maestro Schrodinger package version 11.1.012 (Maestro, 2017) was used for docking calculations. The preparation of crystal structure 3GWS was carried out by the Protein Preparation Wizard tool included in the Maestro suite (Maestro, 2017). This tool automatically sets the bond orders, atomic charges, and the orientation of misoriented groups such as amide groups of Asn and Gln. The crystal structures were protonated and unwanted bound molecules removed from structure. Then, minimization was performed with the OPLS-2005 force field with the aim to relax the atomic coordinates until the geometric convergence (0.01 Å RMSD). The position of ligand in the corresponding crystal structure was used to define the binding site, and the receptor grid generator workflow was used to generate the grid box in the region of 20 Å around the co-crystallized ligand. The docking calculations were performed using the Glide extra precision (XP) mode with the OPLS-2005 force field. There were no protonation states of ligands in the physiological pH of 7.4. For each ligand, 10 docking poses were calculated and other options were left at their default values. The topmost pose from each docking run was included in the analysis. Moreover, Prime-MM/GBSA approach (Schr€ odinger, 2017) was used to calculate the binding free energies of each ligand-receptor complex with an applied force field of OPLS-2005 in Maestro Schrodinger suite. The binding free energy of ligand-receptor complex was calculated using the following equation, The best pose of docking calculations for each ligand-receptor complex was employed as starting model for the ONIOM calculations. The three-layer ONIOM method implemented in the Gaussian 09 program package ) was used to calculate the binding energies of complexes. The high, medium, and low theoretical levels were applied for the three-layer ONIOM calculations, which are respectively considered by the high-level quantum mechanics (HQ), low-level quantum mechanics (LQ), and molecular mechanics (MM) methods. The further details of the threelayer ONIOM methodology are presented in the literatures (Chung et al., 2015;Shahraki & Ebrahimi, 2019). In this work, the ligand atoms were considered in the high-level layer, and the atoms of residues within a sphere with a radius of 10 Å around the ligand were included in the medium-level layer. The medium layer should sufficiently be large to consider all the possible non-bonding interactions between the ligand and the receptor. The low-level layer constitutes other residues of receptor, as shown in Figure 1.
The hybrid meta DFT method M06-2X/6-31G(d), the PM6 semiempirical QM method, and AMBER molecular mechanics were respectively used as the high, medium, and low theoretical levels. The M06-2X functional provides accurate descriptions for both hydrogen bond and dispersive interactions between the ligands and the binding site of receptor (Hohenstein et al., 2008). Besides, polarization functions on heavy atoms are essential for these noncovalent interactions. The ligands were fully optimized, and the atomic positions of the medium-level and low-level layers were frozen during the ONIOM optimization due to avoid unrealistic expansions of the receptor and to keep it close to the X-ray structure. The hydrogen atoms served as link atoms to saturate the dangling bonds between two layers.
The binding interaction energy between the ligand and protein was calculated using the following formula (Kitisripanya et al., 2011): where E complex , E ligand , and E protein are the energy of the ligand-receptor complex, ligand, and protein after ONIOM optimization, respectively.
The MD simulations were also performed on the docked complexes, as the initial coordinates, via the GPU-accelerated GROMACS 2020 (Kohnke et al., 2020). The topology and coordinate files of protein were generated using the pdb2gmx program taking parameters from the AMBER99SB-LIDN force field (Lindorff-Larsen et al., 2010). The antechamber program of Amber Tools 14 (Case et al., 2015) was applied to assign the partial atomic charges of ligands based on the AM1-BCC method (Jakalian et al., 2002). The van der Waals and bonded parameters of ligand were taken from the general amber force field (GAFF) . The full Amber topology and coordinate files of ligands were created using the parmchk and tleap programs implemented in the Amber Tools package (Case et al., 2015). The AMBER format files of ligands were converted to the GROMACS format using the acpype python script (Sousa da Silva & Vranken, 2012). Then, the coordinate and topology files of protein and ligands were merged to obtain the final starting structure and topology file for each complex. Systems were solvated in a cubic periodic box with a side length of 15 Å 3 by the addition of TIP3P water molecules (Jorgensen et al., 1983), and they were neutralized by adding counterions as required. The energy minimization of systems were performed with the steepest descent algorithm for 50,000 steps. Then, an NVT equilibration was performed to heat the systems to 300 K during a 50 ps at a constant volume, and the pressure was equilibrated to 1 atm during a 500 ps NPT simulation at a constant 300 K temperature employing V-rescale thermostat and Berendsen barostat. A cut-off radius of 10 Å was employed for van der Waals interactions, and the longrange electrostatics was treated using a particle mesh Ewald (PME) scheme. The final MD production run was carried out in 120 ns with a 0.2 fs time step. Visualization of protein-ligand complexes and MD trajectories were carried out with the VMD software package (Humphrey et al., 1996).

Results and discussion
The intended compounds have been classified into eight groups as shown in Scheme 2. There are four compounds in each group with X ¼ Br, Cl, F, and H. There are four OH substituents in groups A-F, and two OHs in groups G and H. The substituents locate at the same positions of two rings in each compound (Scheme 2). There are two rotatable bonds, C1-C2 and C1-C2 0 (Scheme 2), in bisphenol and its derivatives. The potential energy surfaces (PESs) were obtained by two-dimensional scans around two dihedral angles, C3C2C1C2 0 (u1) and C2C1C2 0 C3 0 (u2), to determine the global minimum structures. The PESs of compounds are shown in Figure S1 in Supplementary Information (SI). All compounds follow similar patterns in each group. The global minima structures on PESs were further optimized at the B3LYP/6-311þþG(d,p) level to obtain the most stable conformers for latter calculations. The most important geometrical parameters of the optimized structures are shown in Figure S2 in SI. The non-planar structures are identified as the most stable conformer for each compound.

Antioxidant activity
The antioxidant activities of compounds were estimated on the base of three mechanisms: the proton coupled-electron transfer (PC-ET) or H atom transfer (HAT), the sequential proton loss electron transfer (SPLET), and the sequential electron transfer proton transfer (SETPT).
The HAT mechanism is on the base of BDE value, which is a measure of ability to lose a hydrogen atom and forming the radical form of compound. An OAH bond is dissociated easier when the related BDE value is lower. The BDE values of OAH bonds calculated at the B3LYP/6-311þþG(d,p) theoretical level are given in Table S1, and the lowest BDEs of compounds are graphed in Figure 2(a).
According to Figure 2(a), the same sequence A < B < F < E < D is observed for BDEs of compounds in each group, which is followed by C < G % H, G % H < C, H < G % C, and G % H < C for compounds with Br, Cl, F, and H substituents, respectively. Therefore, the lowest BDE values correspond to the compounds of group A (1-4) which have two adjacent OH groups in each ring (Scheme 2). The BDE values of OH substituent located at the position 2 (p2, see Scheme 2) are almost 10 kcal mol À1 lower than those located at the position 1 (p1) for the compounds of group A (Table S1). In other words, homolytic cleavage of OH bond in p2 is more possible than that of p1 in transferring hydrogen atom to a free radical. Hence, the highest antioxidant activity is related to the OH in p2 of group A.
As can be seen in Figure 2(a), the lowest BDE corresponds to the compound 4 (X ¼ H) in group A, where the trend of BDE values is 4 (X ¼ H) < 1 (X ¼ Br) < 3 (X ¼ F) < 2 (X ¼ Cl) in this group. The BDE values increase slightly, by 1-2 kcal mol À1 , from group A to B and from group B to F (see Table  S1 and Figure 2a). Since the OH groups in p1 and p2 are equally affected by halogen substituents in the compound of groups E and F, there are slight differences in the BDE values of these compounds (Table S1).
The BDE values of OH at p1 are lower than p2 for compounds of groups C and D, although the differences are negligible in group D (Table S1). The OH group is located between two halogen substituents at p1 in these groups, which is similar to the OH position in the compounds of groups G and H. Therefore, the BDE values of compounds in group C are approximately close to those in groups D, G, and H (Figure 2a), while are almost 10 kcal mol À1 larger than those in group A (Table S1).
There is only one OH group on each ring, between two halogen substituents, in the compounds of groups G and H that is located (Scheme 2). The halogen substituents locate at the same positions in groups G and H, and the only difference is the bridged group, which is C(CH 3 ) 2 in group H instead of CH 2 in group G (Scheme 2). According to Table S1 and Figure 2a, the negligible differences between the BDE values of groups G and H demonstrate that the antioxidant activities are not significantly affected by the methyl groups locate at the bridge atom. The BDE values of groups G and H follow almost the same orders H < Cl < F < Br and H < F < Cl < Br, respectively.
In the SPLET mechanism, the first and most crucial step is the heterolytic OH bond cleavage that leads to the formation of phenolate anion, where the controlling parameter is the PA. The PA values calculated at the B3LYP/6-311þþG(d,p)  , 5, 9, 13, 17, 21, 25, and 29 calculated at the B3LYP/6-311þþG(d,p) level and mapped onto the isosurface of the electron density (0.001 electrons per Å 3 ). Blue and red regions represent positive and negative potential areas, respectively. level are summarized in Table S1, and the lowest PA value of each compounds are graphed in Figure 2(b). The anion formation is easier and the antioxidant activity is higher for the compounds with lower PAs. As can be seen in Table S1 and Figure 2(b), the PA values are lower in the compounds of groups A and B in comparison with other groups. Thus, that series of compounds are expected to be better antioxidants, which is in agreement with the results obtained from the HAT mechanism. The PA values of compounds in group C are almost the same as those of groups D-H, where the lowest difference is observed for the F substituted compounds (Figure 2(b)). The lowest/highest values of PA correspond to the compounds of group A/F, where the differences are almost 10 kcal mol À1 (Table S1). The trend of PA values is Br < Cl < F < H for all groups (see Figure 2(b)). Therefore, the OAH cleavage could be facilitated in the presence of halogen substituents, where the Br substituent has the highest influence.
According to Table S1, the differences in the PA values are almost 1.50 kcal mol À1 for the compounds of groups G and H. Therefore, the CH 3 groups located on the bridge atom do not affect the antioxidant property significantly.
The second step in the SPLET mechanism is the electron transfer from anion that forms the radical species. The ETE values calculated using Equation (3) are reported in Table S1 and graphed in Figure 2(c). The lowest/highest ETE values correspond to the compounds of group F/G. The trend in the ETE values of all groups is H < F < Cl < Br (Figure 2(c)). The opposite is true for the PA values, because electron transfer is more difficult for more sable anions.
In the SETPT mechanism, the electron donation of compounds is assessed by the ionization enthalpy (IE) according to the first step of mechanism. The single electron transfer to oxidant molecules is easier from a compound with lower IE. The IE values calculated at the B3LYP/6-311þþG(d,p) level are gathered in Table S1, and the values are also graphed in Figure 2(d). The order of IEs is H < Br < Cl < F in all groups, where the lowest values correspond to group F (Figure 2(d)). The electron transfer is also involved in the SPLET mechanism and is described by the ETE parameter.
The SETPT mechanism is followed by the proton transfer from the produced cation radicals. The proton dissociation enthalpies (PDE) calculated at the mentioned level are reported in Table S1, and the lowest PDE values are also graphed in Figure 2(e). The larger PDE values correspond to groups D and F. The lowest PDE value is observed in the presence of F substituent in each group, where the trend is F < Cl < Br < H.

ESP analysis
The molecular electrostatic potential (ESP) (Johnson et al., 1993) is a useful parameter to understand the electrophilic and nucleophilic sites on the compounds. In order to compare the electronic distribution on the compounds, the ESP surfaces were obtained using GaussView 5.0 program . Figure 3 depicts the ESP surfaces of compounds 1, 5, 9, 13, 17, 21, 25, and 29 that belong to groups A-H, respectively. The ESP surfaces of other compounds are presented in Figure S3. As can be seen in Figure 3, the blue and red regions represent the positive (electron-deficient region) and negative (electron-rich region) electrostatic potential areas, respectively. The positive and negative potential areas are observed around the H and O atoms, respectively. Vicinity of two OH groups in the compounds of groups A and B decreases the negative potential around the O atom at p2 (Figure 3). Hence, deprotonation from p2 is easier than p1, which makes the antioxidant activity of p2 larger than p1 in groups A and B.
According to Figure 3, the large positive area, which is observed between two rings of compounds 9, 17, 25, and 29 (in groups C, E, G, and H), increases the electron affinity of compounds. The positive potential (blue color) is slightly lower between two rings of compounds 9 and 17 in comparison to 25 and 29. In addition, the related red areas around the OH group are larger than those of compounds 1 and 5 (groups A and B). Hence, deprotonation of compounds is harder in groups C, E, G, and H in comparison to groups A and B.
The positioning of OH between two halogen substituents increases in the negative potential around the O atom of p1 in compound 13 (group D). Thus, PA is larger at p1 in comparison with p2 in this group. As can be seen in Figure 3, deprotonation of p2 in the group D (13) is easier than those of both p1 and p2 in group C (9), which makes larger the antioxidant activities of group D in comparison to C. The negative potential calculated around the OH at p1 is larger than that at p2 for 21 in the group F. Approximately similar ESP surfaces are observed for the compounds of groups F (21), A (1), and B (5). Accordingly, the similar antioxidant activity is expected for the compounds of those groups.
There is no halogen substituent in compounds 4, 8,12,16,20,24,28, and 32 of groups A-H, so the red area is distributed on a larger region of ESP surface ( Figure S3), which increases the PA of that relative to other compounds. Consequently, deprotonation of these are more difficult in comparison with the others. The results obtained from the ESP surfaces are in accord with the trend of antioxidant activities. Therefore, the ESP surfaces illustrate well the antioxidant capacity of compounds.

Binding affinities to thyroid hormone receptor-b
(THR-b)

Molecular docking
The results obtained from docking calculations on 32 bisphenol derivatives (Scheme 2) using the MOE software (Software Available from Chemical Computing Group Inc., 1010 Sherbrooke Street West, Suite 910, Montreal, Canada H3A 2R759) are gathered in Table 1. The binding energies of  reference ligands, T 0 -T 4 , are also given to make a reasonable comparison between binding affinities. As can be seen in Table 1, the binding affinities of all compounds, with the exception of 4, 8, 16, 28, and 29, are approximately higher than those of reference ligands. The score values are in the range of À10.91 to À15.41 kcal mol À1 that sequentially correspond to compounds 28 (group G) and 2 (group A). Therefore, the intended compounds possess the tolerable ability to bind the active site of receptor. However, no ligand presented free binding energy lower than the reference compound T 2 .
According to Table 1, the lowest scores in each group of ligands corresponds to the compound without any halogen substituents. Therefore, halogen bond interaction is important in interaction between ligands and the binding site residues. The highest scores correspond to the compounds 2 (À15.41), 14 (À15.35), 17 (À15.29), 18 (À15.23), and 9 (À15.04), where the data in the parentheses are the scores in kcal mol À1 , which verify them as the proper ligands to bind with the active site. The docking procedures were repeated to avoid the margin of error; and results gathered in Table S2. The change in the results after the repetitions is negligible. Figure 4 shows the 2D and 3D view of interactions between ligand 2 and the binding site residues. There is at least a halogen bond interaction in the ligand-active site pair, as shown in Figure 4. The topmost pose of 2 is stabilized by two halogen bond interactions between binding site residues Ala234 and Thr329 and compound 2 with distances of 3.28 and 2.90 Å, respectively (Figure 4). Figure S4 shows the 2D view of interactions of ligands 9, 14, 17, and 18 with the binding site residues. The residues Asn331 and Gly344 interact with compound 9 via a halogen bond and a hydrogen bond interaction with distances of 3.11 and 2.76 Å. In addition, Leu330 contributes to a pi-H interaction with ligand 9. One halogen bond interaction is observed between residue Thr239 and ligand 14 with a distance of 2.67 Å, as illustrated in Figure S4.
According to Figure S4, compounds 17 and 18 are stabilized in the binding site via hydrogen bond interactions with the residues Met313 and Asn331, respectively, and the distances of nearest atoms are sequentially equal to 2.57 and 2.82 Å. A pi-H interaction is also observed between residue Leu330 and compounds 17 and 18 in the binding site ( Figures S4).
Docking calculations have also been performed using the Maestro Schrodinger package (Maestro, 2017) to confirm the results of the MOE docking. The docking scores obtained from XP (extra precession) mode in Schrodinger software are reported in Table 1. For the reference ligands T 0 -T 4 , docking scores are in the range of À8.84 to À12.29 kcal mol À1 , while that is changed from À9.05 to À13.64 kcal mol À1 for intended compounds. Therefore, all selected ligands are suitable to bind the active site of receptor. As can be seen in Table 1, the XP score values of all compounds are higher than that of reference ligands T 0 and T 2 . In addition, the affinities of compounds 1-7 and 17-19 to bind the active site are larger than the reference ligand T 3 based on XP scores (Table 1). Moreover, the docking scores of compounds 1-3 and 17-19 are also higher than that of T 4 . Similar to the MOE results, the lowest Glide XP score in each group of compounds, A-H, corresponds to the compound without any halogen substituent except for groups D and H. Figure 5 shows the binding poses of compounds 1-3 and 17-19 that are visualized in 3D pattern using Maestro suite. Compounds are consistently binding to the residues Phe272, Ile275, Met310, Met313, Leu330, Asn331, Gly344, Leu346, and His435 of active site via hydrogen and halogen bond interactions. The residues Phe269,Phe272,Ile275,Ile276,Ala279,Met310,Met313,Ala317,Leu330,Leu341,Leu346,Ile353,Met442,Phe451,and Phe455 are also contributed to bind with ligands via the hydrophobic interactions.
A visual representation of the conformational rotations of protein structure can be given by the Ramachandran plot. Therefore, these plots were provided to validate the docking models, as shown in Figure 6. The dihedral angles of amino acid residues appear as crosses in the plot. The black, dark grey, grey, and light grey regions represent the highly preferred conformations. The less favorable regions are poorly populated due to steric repulsion of adjacent atoms. All torsion angles of protein are approximately within the allowed regions of the plot (see Figure 6). Therefore, the docking models was validated by the Ramachandran plots.

ONIOM calculations
The binding energies (DE) of some ligand-receptor complexes were calculated by the quantum mechanics calculations via the ONIOM methodology. The best four compounds based on the MM/GBSA binding free energies (4, 15, 17, and 18) were chosen for further investigation by a rigorous three-layer ONIOM method. The DE values are graphed in Figure 7. The trend in the DE values is 15 (À18.85) > 4 (À24.05) > 17 (À32.19) > 18 (À43.24 kcal mol À1 ). The more negative value corresponds to 18, in which two OH groups are surrounded by two Cl substituents in each ring (Scheme 2).
The non-bonded interactions of 4 and 15 in the binding site are illustrated in Figure 8. Three hydrogen bond interactions are observed between 4 and the residues Asn331 (2) and His435 in the binding site. The distances of nearest atoms involved in these interactions are equal to 2.20, 1.95, and 2.89 Å, respectively. Also, the hydrophobic nature of binding site cause hydrophobic interactions between 4 and Ile276, Ala279, Met310, Met313, Ala317, Leu330, and Leu346 residues.
The residues Met313, Asn331, and His435 participate in both hydrogen and halogen bond interactions with 15 in the binding site of receptor (Figure 8). The hydrogen and halogen bond lengths are in the range of 1.72-2.33 and 2.28-2.87 Å, respectively. Moreover, the residues Ile276, Ala279, Met310, Met313, Ala317, Leu330, and Leu346 bind with 15 via the hydrophobic interactions (Figure 8).
The orientation of 17 in the binding site was suitable to interact with the Met313, Asn331, and His435 residues via the four hydrogen bond interactions, as shown in Figure 9. Also, the residue Gly344 participate in a halogen bond interaction with ligand. Moreover, several hydrophobic interactions are observed between ligand and residues Phe272, Ile275, Ile276, Ala279, Met310, Met313, Ala317, Leu346, Ile353, and Met442 in the active site.

MD simulations
The 3GWS protein and the two best docked complexes, 3GWS-17 and 3GWS-18, were subjected to MD simulations with time periods of 120 ns. The stability of systems was examined by further analysis of MD trajectories using root mean square deviation (RMSD), root mean square fluctuation (RMSF), and hydrogen bond analysis. The RMSD of protein Ca was analyzed to evaluate the stability and conformational changes of protein during the course of simulations for 3GWS-17 and 3GWS-18 complexes and the ligand-free 3GWS protein. All systems show relatively stable RMSD fluctuations after 50 ns simulation, as shown in Figure 10(a). Therefore, those are equilibrated after this time. There are slight differences in the RMSD values for complexes and ligand-free THR-b protein, which indicate small conformational changes of THR-b in these complexes. However, the structural changes of protein are not negligible upon binding of ligands 17 1nd 18. After equilibrating of systems, the RMSD values for 3GWS-18 complex are lower than those of 3GWS-17 that verify the more stability of 3GWS-18 complex.
The RMSF analysis of residues can be used to evaluate the structural movement and flexibility of THR-b upon binding of ligands 17 and 18. Figure 10 The hydrogen bonds, the most important directional interactions in receptor-ligand complexes, cause the rigidity of receptor structure and specificity of intermolecular interactions. Figure 10(c) and (d) depict the number of hydrogen bonds formed between receptor and compounds 17 and 18. It is observed throughout the simulation time that the hydrogen bond network of both complexes are stable.

Conclusions
The antioxidant activity of 32 natural and synthetic halogenated bisphenol derivatives have been studied by DFT calculations through hydrogen atom transfer, sequential electron transfer proton transfer, and sequential proton loss electron transfer mechanisms. The halogen substituents enhanced the antioxidant activity of compounds, where the lowest activity corresponds to the compounds without any halogen substituents. The antioxidant activities of compounds are influenced by the distance between the position of OH groups and halogen substituents. Also, the number of OH groups is an important factor that affects the antioxidant capacity of compounds; according to all three mechanisms, the lowest activity corresponds to the structures with only one OH group in each ring. Furthermore, two adjacent hydroxyl groups increase the antioxidant activity of compounds. The vicinity of two OH groups decreases the negative potential around the O atom and simplifies the deprotonation. Moreover, the antioxidant activity is not significantly affected by the methyl groups locate at the bridge atom. The results obtained from the ESP surfaces are in agreement with the antioxidant activity of compounds. Therefore, the ESP surfaces illustrate well the antioxidant capacity of compounds.
Furthermore, the binding affinity of compounds against thyroid hormone receptor-b was considered using the computational approaches including molecular docking, MM/ GBSA, and ONIOM calculations. According to the MOE docking, compounds 2, 9, 14, 17, and 18 are the most suitable ligands to interact with the main residues Asn233, Ala234, Phe272, Ala279, Met310, Met313, Arg316, Arg320, Thr329, Leu330, Asn331, Leu341, Gly344, and Leu346 of the binding site of THR-b via the hydrogen bond, halogen bond, and pi-H interactions. The halogen bond is an important interaction that binds the ligands with the active site of THR-b.
In addition, the three-layer ONIOM calculations were performed on the best four compounds 4, 15, 17, and 18 detected based on binding free energies. Those results show that the affinities of compounds 17 and 18 are higher for binding to THR-b. Compounds bind with the active site via the residues Phe272, Ile275, Ile276, Met313, Asn331, Gly344, and His 435. The hydrophobic nature of binding site causes several hydrophobic interactions with the ligands, where the residues Ile276, Ala279, Met310, Met313, Ala317, Leu330, and Leu346 are common to take part in those interactions.
The affinity of compounds 17 and 18 to bind to the THRb cavity was also verified by 120 ns MD simulation. The results indicated the stable hydrogen bond network between ligands and residues in the binding site across the simulation time, the tight binding of ligands with low fluctuations of binding site residues, and the RMSD values lower than 0.5 in equilibration zone.
The results of docking obtained from both softwares, MM/ GBSA free energies, MD simulations, and more rigorous ONIOM calculations verified compounds 17 and 18 as the most potent candidates to bind with the ligand-active site of THR-b.

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

Funding
We thank the University of Sistan and Baluchestan for financial support and Computational Quantum Chemistry Laboratory for computational facilities.