In silico investigations of anti-androgen activity of polychlorinated biphenyls

Chemosphere 92 (2013) 795–802

Contents lists available at SciVerse ScienceDirect

Chemosphere journal homepage: www.elsevier.com/locate/chemosphere

In silico investigations of anti-androgen activity of polychlorinated biphenyls Xiaolin Li a, Li Ye b, Xiaoxiang Wang a, Wei Shi a,⇑, Hongling Liu a, Xiangping Qian b,c, Yongliang Zhu b, Hongxia Yu a,⇑ a b c

State Key Laboratory of Pollution Control and Resource Reuse, School of the Environment, Nanjing University, Nanjing 210023, PR China Suzhou NeuPharma Co., Ltd., Suzhou 215123, PR China College of Pharmaceutical Sciences, Soochow University, Suzhou 215123, PR China

h i g h l i g h t s  The AR-antagonistic activity of PCBs was studied by 3D-QSAR.  The binding mode between PCBs and AR were investigated.  The combination of 3D-QSAR, molecular docking, and MD simulations was performed.

a r t i c l e

i n f o

Article history: Received 1 June 2012 Received in revised form 11 April 2013 Accepted 13 April 2013 Available online 9 May 2013 Keywords: Polychlorinated biphenyls Androgen receptor 3D-QSAR Molecular docking Molecular dynamics simulations

a b s t r a c t Polychlorinated biphenyls (PCBs) have attracted great concern as global environmental pollutants and representative endocrine disruptors. In this work, a molecular model study combining three-dimensional quantitative structure-activity relationship (3D-QSAR), molecular docking, and molecular dynamics (MD) simulations was performed to explore the structural requirement for the anti-androgen activities of PCBs and to reveal the binding mode between the PCBs and androgen receptor (AR). The best comparative molecular similarity indices analysis (CoMSIA) model, obtained from receptor-based alignment, shows leave-one-out cross-validated correlation coefficient (q2) of 0.665 and conventional correlation coefficient (R2) of 0.945. The developed model has a highly predictive ability in both internal and external validation. Furthermore, the interaction mechanisms of PCBs to AR were analyzed by molecular docking and MD simulation. Molecular docking indicated that all the PCBs in the data set docked in a hydrophobic pocket. The Binding free energies calculated by Molecular mechanics–Poisson Boltzmann surface area (MM–PBSA) not only exhibited a good correlation with the experimental activity, but also could explain the activity difference of the studied compounds. The binding free energy decomposition analysis indicates that the van der Waals interaction is the major driving force for the binding process. Ó 2013 Elsevier Ltd. All rights reserved.

1. Introduction Polychlorinated biphenyls (PCBs) are an important group of synthetic organic compounds, which include 209 congeners characterized by the number and position of the chlorine atoms on the biphenyl core. Despite their production and use were banned in most countries, PCBs remain an environmental problem due to their stability, mobility and ongoing leaking to the environment from old appliances and electrical equipments. Due to their stability and lipophilic nature, PCBs concentrate in lipids and adipose tissues and accumulate in the food chain (Bremle et al., 1995; Portigal et al., 2002). ⇑ Corresponding authors. Tel.: +86 25 89680617; fax: +86 25 89680356. E-mail addresses: [email protected] (W. Shi), [email protected] (H. Yu). 0045-6535/$ - see front matter Ó 2013 Elsevier Ltd. All rights reserved. http://dx.doi.org/10.1016/j.chemosphere.2013.04.022

PCBs are one of the most widespread and persistent groups of endocrine disrupting chemicals (EDCs) in the environment. A number of endocrine related disturbances linked to exposure for PCBs have received much attention in the last years (Decastro et al., 2006; Takigami et al., 2010). Most studies evaluating the endocrine disrupting activity of PCB have focused on the aryl hydrocarbon receptor (AhR) and the estrogen receptor (ER) as molecular targets (Besselink et al., 1998; Layton et al., 2002; Oh et al., 2007). However, relatively fewer studies have addressed the androgen disrupting activity of PCBs. Up to now, large-scale in vitro assessment of 209 PCB congeners remains a labor intensive, time-consuming and expensive operation. Besides, due to some purified PCB standards are difficult to be obtained, their true effects cannot be concluded in many studies. A more efficient and cost-effective alternative method would be to employ in silico molecular modeling approaches for the purpose of study their androgen receptor

796

X. Li et al. / Chemosphere 92 (2013) 795–802

(AR)-antagonistic activity. Hence, as a kind of well-known modeling method, quantitative structure-activity relationship (QSAR) is suggested to model and predict the AR-antagonistic activity of PCBs. QSAR has been successfully used to predict the toxicity of chemicals (Cheng et al., 2011), endocrine disrupting activities of polybrominated diphenyl ether (PBDE) and Hydroxylated PBDE (Li et al., 2010; Yang et al., 2011a), and AhR binding affinity of PCBs (Li et al., 2011). However, there is a lack of in silico studies on the AR-antagonistic activity of PCBs. Besides, the research focused on the molecular level understanding of the interaction mode of PCBs and the protein of AR remains scarce. Similar three-dimensional QSAR (3D-QSAR) models by comparative molecular similarity indices analysis (CoMSIA) method have been previously developed for prediction of AR-antagonistic activity of PBDEs and identification of the interactions of PBDEs with the AR (Yang et al., 2009), but the detailed information about the ligand–protein interactions does not involved. In the present study, a combined computational approach by 3D-QSAR, molecular docking, and molecular dynamics (MD) simulations was performed to reveal the key structural features influencing the AR-antagonistic activity of PCBs and investigate the detailed binding mode between PCBs and AR. To derive the rational conformation for developing 3D-QSAR models, we applied two different alignment rules; ligand-based and receptor-based alignments. The reliability and robustness of the developed best model were estimated by bootstrapping analysis and leavemany-out (LMO) cross-validation analysis. Moreover, the detailed interaction mechanisms of two PCBs with much different activity were analyzed on the basis of the results of molecular docking and MD simulation. The interaction features between both compounds and AR were determined by analyzing energy contributions to the binding free energy of the protein–ligand complexes. The derived results can be utilized to accurately predict the ARantagonistic activity of related analogues and provide some insights into the interaction mechanism of PCBs to AR.

Table 1 Structure and activity of the target compounds.

a

2. Materials and methods 2.1. Data set The structures of the PCBs and their biological activities were taken from the literature (Hamers et al., 2011). The AR-antagonistic potencies (IC50) of individual PCBs were determined in the ARCALUX bioassay (chemically activated luciferase gene expression). The AR-CALUX bioassay use reporter cell lines carrying a luciferase gene under the transcriptional control of response elements for activated receptors. Compound PCB-153 (IC50 = 8.5) was not selected for CoMSIA analysis to avoid the over fitting of the model (Liu et al., 2010). The total set of compounds (23 compounds) was divided into a training set (containing 19 compounds) for model generation and a test set (containing 4 compounds) for model validation. The test set compounds were selected by considering whether they can represent structural diversity and the distribution of biological activities of the training set appropriately. The structures and biological activities of PCBs are shown in Table 1. 2.2. Molecular modeling and docking The 3D-QSAR and molecular docking computations were performed by using Sybyl 7.3 molecular modeling software package (Tripos Inc, St. Louis, MO). The 3D structure of each compound was constructed using the Sketch Molecule module in Sybyl. The geometry of these compounds was subsequently optimized using Tripos force field (Clark et al., 1989) with the Gasteiger–Hückel

Compound

Substitution order

PCB-19a PCB-28 PCB31 PCB-47 PCB-51 PCB-52a PCB-53 PCB-74 PCB77 PCB-80 PCB-81 PCB-95 PCB-100 PCB-101 PCB-104 PCB-105 PCB-114 PCB-118 PCB-122 PCB-123 PCB-125 PCB-126 PCB-128 PCB-136a PCB-138 PCB-156 PCB-157 PCB-167 PCB-168 PCB-169 PCB-170a PCB-180 PCB-189 PCB-190 PCB-194 PCB-206 PCB-209

2, 20 , 6 2, 4, 40 2, 40 , 5 2, 20 , 4, 40 2, 20 , 4, 60 2, 20 , 5, 50 2, 20 , 5, 60 2, 40 , 40 , 5 3, 30 , 4, 40 3, 30 , 5, 50 3, 4, 40 , 5 2, 20 , 3, 5, 6 2, 20 , 4, 40 , 6 2, 20 , 4, 5, 50 2, 20 , 4, 6, 60 2, 3, 30 , 4, 40 2, 3, 4, 40 , 5 2, 30 , 4, 40 , 5 20 , 3, 30 , 4, 5 2, 30 , 4, 40 , 50 20 , 3, 4, 5, 60 3, 30 , 4, 40 5 2, 20 , 3, 30 , 4, 40 2, 20 , 3, 30 , 6, 60 2, 20 , 3, 4, 40 , 50 2, 3, 30 , 4, 40 , 5 2, 3, 30 , 4, 40 , 50 2, 30 , 4, 40 , 5, 50 2, 30 , 4, 40 , 50 , 6 3, 30 , 4, 40 , 5, 50 2, 20 , 3, 30 , 4, 40 , 5 2, 20 , 3, 4, 40 , 5, 50 2, 3, 30 , 4, 40 , 5, 50 2, 3, 30 , 4, 40 , 5, 6 2, 20 , 3, 30 , 4, 40 , 5, 50 2, 20 , 3, 30 , 4, 40 , 5, 50 , 6 2, 20 , 3, 30 , 4, 40 , 5, 50 , 6, 60

IC50 (lM)

Residual

exp.

pred.

0.38 0.76

0.75 0.71 1.65 0.76 0.89 2.7 1.2 2.0 1.19 1.5 0.95 1.4 0.6 2.5 1.8 0.99 1.26 0.51 0.33 1.12 0.22 0.64 0.41 1.00 0.98 1.16 0.65 1.19 0.18 1.29 0.71 2.3 1.1 1.5 1.63 1.53 1.61

0.60 0.50 3.2 1.3 2.1 1.6 1.4 1.0 2.7 1.9

0.47 0.39 0.25 0.53 0.48 0.71 1.0

0.16 1.1 2.1 1.3

0.37 0.05 0.16 0.39 0.5 0.1 0.1 0.1 0 0.4 0.2 0.1

0.04 0.06 0.03 0.11 0.07 0.29 0.02

0.02 0.39 0.2 0.2

Test set.

charges (Gasteiger and Marsili, 1980). Repeated minimizations were performed by Powell method with a maximum iteration of 1000 to reach a energy convergence gradient value of 0.001 kcal mol1 Å. The minimized structures were used as initial conformations for molecular docking. Molecular docking simulations of PCBs into the binding pocket of AR were performed using the Surflex-dock program interfaced with Sybyl in this study. The detailed algorithm for Surflex-Dock is presented in the literature (Jain, 2007). The crystal structure of human AR complexed with hydroxyflutamide (PDB code: 2AX6) used in Surflex-Dock was obtained from the Protein Data Bank (http://www.rcsb.org/pdb/). Prior to docking, receptor was prepared by extracted hydroxyflutamide and structural water molecules, added polar hydrogen atoms, and assigned Kollman-all charges to protein atoms. During the docking process, automaticbased docking mode was applied and parameters with default values were employed in all operations. In this study, 10 conformations were obtained though Surflex-dock for each ligand. The binding conformation assumed to represent the bioactive conformation of ligand was selected based on the following two criteria: (i) the conformation gained the highest total docking score, and (ii) the orientation of the conformation was similar to that of the natural ligand hydroxyflutamide.

2.3. Molecular alignment and CoMSIA analysis Molecular alignment of compounds is of vital importance for 3D-QSAR studies (Liu et al., 2012). In the current work, two

X. Li et al. / Chemosphere 92 (2013) 795–802

different alignment rules (ligand-based and receptor-based alignment) were adopted to derive valid and reliable 3D-QSAR models. For ligand-based alignment, the most potent compound in the data set (PCB-168) was selected as the template to fit the remaining compounds by using of Align Database command in Sybyl. For receptor-based alignment, the bioactive conformation of each chemical firstly derived from molecular docking was assigned MMFF94 charges and imported into a Sybyl molecular database for CoMSIA study without further energy minimization. The CoMSIA models were developed in the Sybyl with the default parameters. CoMSIA descriptor fields, including steric, electrostatic, hydrophobic, hydrogen bond-donor and hydrogen bond-acceptor fields, were calculated using sp3 carbon carrying a charge of +1.0 as the probe atom. The attenuation factor (a) was set to default value of 0.3. Since PCBs cannot form hydrogen bond with amino acid residues in the binding pocket of AR (see below), only the other three descriptor fields were used as independent variables in the following CoMSIA regression analysis.CoMSIA regression models were built by partial least squares (PLSs) regression (Wold et al., 1984). In the PLS regression analysis, a leave-oneout (LOO) cross-validation analysis was performed to determine the optimum number of components (ONCs) and cross-validated correlation coefficient (q2). Then non-cross-validation analysis was done with ONC and threshold column filtering of 2.0 kcal mol1 to generate the regression model for CoMSIA. CoMSIA models were considered reliable and acceptable if the q2 was greater than 0.50 and the correlation coefficient (R2) was greater than 0.90 (Diao et al., 2010; Wang et al., 2005). To refine the obtained models, region focusing (Gilbert et al., 2007; Hevener et al., 2008) was performed on conventional CoMSIA models. 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 (Yang et al., 2010). Subsequently, bootstrapping analysis for 100 runs (Yang et al., 2011b), leavemany-out (10 groups, LMO10; and five groups, LMO5) procedures repeated 25 times (Trossini et al., 2009), and Y-randomization performed 10 times (Mouchlis et al., 2012) were used as rigorous tests to further assess the robustness and statistical significance of the derived models. Finally, the external predictive ability of the models was evaluated by the test set. The AR-antagonistic activities of compounds in the test set were predicted using models produced by training set and the predictive R2 (denoted by R2pred ) was calculated. 2.4. Molecular dynamics (MD) simulations The MD simulations were carried out using AMBER 10 software package (Case et al., 2005). The docked structures of AR with two PCBs (highly active compound PCB-168 and lowly active compound PCB-52) systems were used as the initial structures for MD simulations. The force field parameters for the protein were generated by FF03 AMBER force field (Duan et al., 2003). For the ligand, the GAFF parameter assignments (Wang et al., 2004) were made by using Antechamber program (Wang et al., 2004), and the partial charges were assigned by using the AM1-BCC method (Jakalian et al., 2000), respectively. Each protein–ligand complex was electroneutralized by adding suitable sodium counterions. The complexes were subsequently solvated in a global TIP3P (Jorgensen et al., 1983) water molecules with a margin distance of 29 Å. The particle mesh ewald (PME) method (Essmann et al., 1995) was used to deal with the long-range electrostatic interactions. In all simulations, a cutoff distance of 12 Å was set for the non-bonded interactions (Sun et al., 2011). After the starting systems had been treated, in order to relax any steric conflicts, the systems were energy minimized before MD simulations by three steps: First, the hydrogen atoms, the sol-

797

vent, and the ligand and solvent were minimized in three consecutive rounds, each of which consisted of 10 000 steps procedure (1000 steps using the steepest descent method and 9000 steps using the conjugate gradient method), while the rests were restrained with a force of 10 kcal mol1 Å2. Second, with constraint to the solute (10 kcal mol1 Å2), both systems were gradually heated from 0 to 200 K by 25 000 steps to ensure that water molecules were fully optimized. Finally, each system was minimized with no restraint by 1000 steepest descent steps and 49 000 conjugate gradient steps. Then, MD simulations were performed as followings: First, the system was gradually heated from 0 to 300 K over 50 ps and maintained at 300 K followed by 200 ps of constant pressure equilibration. Second, a production run for 3 ns (3000 ps) was performed at 300 K with 1.0 atm pressure. The dynamics equilibration was monitored by check the stability of temperature and pressure of the system. Coordinate trajectories were recorded every 1 ps throughout all MD runs. A total of 200 snapshots of the simulated structures with counterions and water molecules stripped in the last 1 ns stable MD trajectory at 5 ps intervals were extracted to perform the following binding free energy calculations. 2.5. Binding free energy calculations and energy decomposition For each protein–ligand complex, the binding free energy (DGbind) was calculated by the Molecular mechanics–Poisson Boltzmann surface area (MM–PBSA) approach (Kollman et al., 2000) implemented in the Amber. The aforementioned 200 snapshots extracted from the last 1 ns of the MD trajectory were used for calculations. The total free energy of binding was composed of the following molecular species (complex, protein, and ligand) as follows:

DGbind ¼ Gcpx  ðGpro þ Glig Þ ¼ DEMM þ DGsol  T DS

ð1Þ

where DGcpx, DGpro, and DGlig represent the free energy of protein– ligand complex, protein, and ligand, respectively. DGbind was composed of three parts: the molecular mechanical (MM) free energy (DEMM), the solvation free energy (DGsol), and entropy contribution to the free energy (TDS). The DEMM was the sum of the internal energy of bonds, angles and torsions (DEval), electrostatic interaction energy (DEele), and van der Waals interaction energy (DEvdw):

DEMM ¼ DEval þ DEele þ DEvdw

ð2Þ

DGsol was divided into two parts, the polar solvation free energy (DGp), and the nonpolar solvation free energy (DGnp):

DGsol ¼ DGp þ DGnp

ð3Þ

where the DGp was determined by Poisson Boltzmann (PB) equation. The values of the dielectric constant for solute and solvent were set to 1 and 80, respectively. The DGnp was determined according to:

DGnp ¼ cSASA þ b

ð4Þ

where SASA was the solvent accessible surface area. The solvation parameters of c and were set to 0.00542 kcal mol1 Å2 and 0.92 kcal mol1, respectively. As the calculation of entropy contribution was time-consuming and several earlier literatures (Zhou et al., 2005; Zhao et al., 2008; Shen et al., 2012) have shown that entropy does not contribute much to the binding free energy to the same protein, so the entropy contribution was not calculated in this work. Furthermore, the binding free energy decomposition was performed to evaluate the contribution of each residue to the binding process. Given the Generalizd Born (GB) model is much faster than the PB model and allowing the decomposition of polar solvation

798

X. Li et al. / Chemosphere 92 (2013) 795–802

free energy into atomic contributions in a straightforward manner (Nimmanpipug et al., 2011; Yang et al., 2011c), the energy decomposition was performed using MM-GBSA method based on the same 200 snapshots as those adopted in the binding free energy calculations. 3. Results and discussion 3.1. Molecular docking Before docking analysis, it is necessary to verify the rationality of the docking parameters specified for docking method and docking accuracy (Liu et al., 2011). The most straightforward method to evaluate the performance of the molecular docking is redock the cocrystallized ligand into the binding pocket of receptor. The docking method is reasonable if the root mean square deviation (RMSD) value between the best-scored conformation and the experimental one is less than 2.0 Å (Wang et al., 2003; Shen et al., 2012). A low RMSD 0.77 Å for the heavy atoms of the top-scored conformation was obtained by redock hydroxyflutamide into the binding pocket of AR, indicating that the docking parameters specified in Surflexdock method were reliable and can recover well the interaction modes for observed natural AR complex. Then, all compounds in the data set were docked into the active site of AR by the same way. As shown in Fig. 1a, all 23 PCBs in the data set are docked into the binding pocket of AR. All compounds were well placed in the binding pocket in the same way, and share a similar binding mode. The binding mode of the most active compound PCB-168 docked into the AR is taken as an example and illustrated in Fig. 1b. Accordingly, in the binding pocket, PCB-168 did not form hydrogen bonds with any amino acid residues of AR. Interaction between the activation function 2 (AF2) region of AR and some co-activator is necessary for activating AR. Similar to PBDE, no androgenic activity was detected for the studied PCBs could not induce the AF2 region into the active conformation. As observed from the figure, PCB-168 interacted with the receptor AR via hydrophobic contacts and electrostatic interactions. The A benzene ring of PCB-168 was posi-

tioned in a basically hydrophobic pocket which was formed by Leu707, Met745, Val746, Met749, and Phe764. The B benzene ring formed hydrophobic contacts with residues Leu701, Leu873, Phe876, Ala877, and Leu880. These residues were the important residues in the binding pocket and were the main contributors to the ligand and receptor interactions. The interactions between a set of PBDEs and AR were analyzed by Yang et al. (2009) using molecular docking, and residues Leu707 and Met745 were found also provided favorable hydrophobic contacts for the binding of PBDEs and AR. However, there were also some residues (Leu 704 and Met 742), which only appeared in the docking results of PBDEs and AR. These differences may due to the ether bond of PBDEs make the benzene rings more susceptible to the impact of residues and easier to rotate and move. 3.2. 3D-QSAR analysis 3.2.1. CoMSIA statistical results In the present study, two alignment rules (ligand-based and receptor-based) were adopted to derive reliable models. The statistical parameters of the 3D-QSAR models obtained from CoMSIA for the AR-antagonistic activity of PCBs are summarized in Table 2. As shown in Table 2, both the ligand-based and receptor-based CoMSIA model exhibited fewer satisfactory results (q2 < 0.5, R2 < 0.9). Subsequently, these less significant models were optimized through the region focusing method with a discriminant power value of 1.0 and a grid spacing of 0.5 Å. The statistical results in Table 2 showed that the receptor-based CoMSIA model with region focusing had a better statistical result. Therefore, we focused on this CoMSIA model in the following discussion. In the CoMSIA model, the cross-validation PLS analysis for the training set generated a q2 of 0.665 with ONC of 6. The non-cross-validation PLS regression sequentially yielded a R2 of 0.945 with standard error of estimate (SEE) of 0.208. The field contributions of steric, electrostatic, and hydrophobic were 9.7%, 70.9%, and 19.4%, respectively. The results illuminated that electrostatic and hydrophobic contributions shared the large part to the activity.Furthermore, the robustness and stability of the model were verified by internal

Fig. 1. Docking modes between PCBs and AR: (a) all 23 compounds with AR-antagonistic activity determined by in vitro bioassays (Hamers et al., 2011), and (b) PCB-168 at the binding pocket.

799

X. Li et al. / Chemosphere 92 (2013) 795–802 Table 2 Statistical results of CoMSIA models.

CoMSIA (without region focusing) PLS statistics ONC 4 0.312 q2 r2 0.879 SEE 0.293 F 23.530 Field contribution% S 1.0 E 82.3 H 16.7 a

Receptor-based

Training set Test set

2.5

CoMSIA (with region focusing)

CoMSIA (without region focusing)

CoMSIA (with region focusing)

6 0.120 0.892 0.300 15.207

2 0.421 0.890 0.260 60.396

6 0.665 0.945 0.208 33.457

1.8 72.4 25.8

7.5 73.9 18.7

9.7 70.9 19.4

a

Abbreviations: S (steric), E (electrostatic), H (hydrophobic).

validation strategies such as bootstrapping analysis and LMO cross-validation. The bootstrapping analysis for 100 runs was performed and a bootstrapped R2boot value of 0.987 and a SEEboot value of 0.103 were obtained. LMO10 and LMO5 cross-validation procedures were carried out for 25 runs and the mean q2 value of LMO10 and LMO5 were 0.632 and 0.603, respectively. The q2 values of LMO analysis were comparable to those of LOO analysis, which further confirmed that the obtained CoMSIA model was stable and possessed reasonable internal predictive ability. The Y vector (IC50 values) was randomized for 10 times. The obtained q2 and r2 values were in the range of 0.940 to 0.319 and 0.233 to 0.481, respectively. The low values of q2 and r2 indicate that the results from the CoMSIA model were not due to chance correlation.The external validation method, as the most valuable validation method (Zhu et al., 2009; Hao et al., 2011b), was applied to assess predictive ability of the generated model. An external test set was used to validate the reliabilities of the built model. The R2pred value of test set was 0.879, indicating the model have good predictive ability. The predicted activity values for the training set and the test set by the optimal CoMSIA model are listed in Table 1. The correlations between the experimental and predicted activity values are depicted in Fig. 2. With the CoMSIA-based QSAR model, the unknown AR-antagonistic activities of fourteen selected PCBs which were the most commonly searched for in humans exposed to organochlorinated compounds according to the WHO classification (ATSDR, 2000) were predicted in Table 1.

3.2.2. CoMSIA contour maps The results of the CoMSIA model were viewed as 3D colorcoded contour maps. The contour maps of the StDevCoefficient plots for the steric, electrostatic, and hydrophobic fields from the CoMSIA model are shown in Fig. 3, using PCB-168 as a reference structure. The contribution of favorable and unfavorable levels was maintained as default value of 80% and 20%, respectively. The CoMSIA steric contour map is displayed in Fig. 3a, the green and yellow contours represented the steric favorable and unfavorable regions, respectively. As seen from this figure, a small greencolored contour was mapped near the 3-position of A ring. This suggested that bulkier groups were favored at this position, as illustrated by the fact that PCB-80 was more potent in activity than compound PCB-52. The higher potency of PCB-122 than PCB-104 was also such a case. A large green contour covered the 30 - and 40 -positions of B ring, indicated that these positions were prefer for larger substitutes to increase activity. This may explain the difference of experimental activity of PCB-168 and PCB-136 as well as PCB-126 and PCB-101. In addition, the medium sized region of yel-

IC50 (pred.)

Ligand-based

3.0

2.0 1.5 1.0 0.5 0.0 0.0

0.5

1.0

1.5

2.0

2.5

3.0

IC50 (exp.) Fig. 2. Plot of predicted IC50 values versus the corresponding experimental values for the training and test set compounds.

low contour located near 50 -positions of B ring indicated that bulkier group at this position may decrease the activity. For electrostatic field, Fig. 3b depicts the contour maps of favorable electropositive (blue1) and favorable electronegative (red) regions. Two blue regions cover the 2- and 4-positions of A ring represented that the electropositive groups are favorable for enhancing activity. These key positive-favorable properties were consistent with the experimental results, e.g. negative charged Cl atom in these positions will decrease the activity. The fact that PCB-122 had higher potency than PCB-101 and PCB-19 was more potent than PCB-28 verified this conclusion. Besides, a blue region was found near the 30 -positions of B ring, suggested that the presence of Cl atom herein was also harmful to the activity. This may explain why PCB-125 exhibit higher activity than PCB-126. In addition, two red contours located on the 20 - and 50 -positions of B ring demonstrated that electronegative substituents were favorable for the improvement of activity. The hydrophobic contour map of the CoMSIA model is presented in Fig. 3c, where the yellow regions was hydrophobic favorable for the activity, while the white region was hydrophilic favorable. A yellow region was observed close to the 40 -positions of B ring indicated that hydrophobic group here was good for increasing the potency, which was in agreement with the fact that PCB-74 was more potential than PCB-101. The activity discrepancies between PCB-28 and PCB-52 can also be explained by this contour. It is shown that this position is pointing toward the hydrophobic pocket, interacting with the Leu701, Leu873, Phe876, Ala877, and Leu880 regions. This result was consistent with that of molecular docking. Additionally, two white contours were mapped near the 2- positions of A ring and 50 -positions of B ring, implied that hydrophilic groups were favored at these positions. The white contour near the 2- positions of A ring was well consistent with the higher activity of PCB-122 than PCB-100, and higher activity of PCB-128 than PCB-138 was due to the white contour near the 50 - positions of B ring.

3.3. MD simulations In order to compare the difference between the most active and the least active compounds, 3.25 ns MD simulations were performed on the docking complexes of PCB-52-AR system and PCB168-AR system considering the effects of the flexibility of receptor and the water solvation. 1 For interpretation of color in Fig. 3, the reader is referred to the web version of this article.

800

X. Li et al. / Chemosphere 92 (2013) 795–802

Fig. 3. CoMSIA contour maps based on the most active compound PCB-168: (a) steric contour map, (b) electrostatic contour map, and (c) hydrophobic contour map.

Fig. 4. Superimposition of the average structure from the last 1 ns of the MD simulation (magenta) and the initial structure (cyan): (a) PCB-52-AR, and (b) PCB-168-AR.

3.3.1. Validation of the MD simulations Firstly, after about 150 ps MD, the RMSD of the backbone atoms of the systems with respect to the starting structure reached equilibrium (Supplementary material Fig. S1), The RMSD reached 0.22 and 0.18 Å for PCB-52-AR system and PCB-168-AR system, respectively, and did not alter meaningfully in the simulation process. The RMSD values indicating the structures of the developed model were stable and equilibrated during the MD simulations. Secondly, to compare the obtained MD-simulated binding modes and the docking modes of the both systems, a superimposition of the average structures of the 200 snapshots extracted from the last 1 ns of the MD trajectory and the initial docked structures of the both complexes is depicted in Fig. 4. In this picture, initial structure of the docked complex was represented by cyan stick, and the average structure of MD-simulated structure was represented by magenta stick, respectively. It can be seen that the average structure obtained from MD simulation was in good agreement

with the docked structure for both systems, suggesting the sampling method was reasonable and the docking model was valid. Finally, to further validate the reliability of the simulation, the binding affinities of both compounds with AR were calculated using the MM–PBSA method. The results include the overall binding free energy and all of the energy terms for both compounds are summarized in Table 3. As listed in Table 3, there was a large difference of the calculated binding free energy between PCB-52 (21.16 kcal mol1) and PCB-168 (29.61 kcal mol1), which illuminated the difference of binding affinity. The IC50 of PCB-168 was about 20 times greater than that of PCB-52. Consequently, the ranking of the calculated binding free energies were in accordance with the experimental order (Hao et al., 2011a). The consistency between the computational results and the experimental activity suggested that the computational models should be reliable.

X. Li et al. / Chemosphere 92 (2013) 795–802 Table 3 Binding free energy (kcal mol1) of both systems. Energy

PCB-52-AR

PCB-168-AR

DEele DEvdw DEval DEMM DG p DGnp DGsol DGbind

2.63 36.20 0 33.57 14.65 2.24 12.41 21.16

1.96 41.80 0 39.84 12.70 2.47 10.23 29.61

3.3.2. Analysis of the interaction mechanism According to the free energy calculation, PCB-168 had stronger binding affinity. As shown in from Table 3, the individual components of the binding free energies indicated that the van der Waals interactions were the largest contribution to the free energy The nonpolar salvation contribution terms, which correspond to the burial of solvent accessible surface area upon binding, contributed slightly favorably. The binding processes within the ligand–protein complexes were counteracted by unfavorable electrostatics interactions and polar solvation free energy terms. In the two studied complexes, the highly favorable nonpolar interaction contribution (the sum of DEvdw and DGnp) provided the dominant driving force for binding, this phenomenon probably owing to burial effects of ligand’s hydrophobic groups upon binding and the solute-water dispersion interaction. The polar interaction contribution (the sum of DEele and DGp) between ligand and receptor were unfavorable for both systems, which indicated that the large polar energy between the solvent (water molecules) and the ligands counteracted the binding. Comparing PCB-52-AR system and PCB-168-AR system, it showed that the nonpolar interaction contribution and polar interaction contribution of the PCB-168-AR system were 5.83 kcal mol1 and 2.62 kcal mol1 less than those of PCB-52-AR system, respectively. The difference can partly explain the reduced binding affinity of PCB-52 and can be correlated to activity of PCBs. PCB-52 displayed lower activity than PCB-168 because of different number and position of Cl atom between these two compounds. According to Table 3, we can conclude that this difference was likely to cause the reduction of van der Waals interaction energy. This result was in accordance with the aforementioned 3D-QSAR study which indicated that bulky substitutes at the 30 - and 40 -positions and hydrophobic group at the 40 -position are favorable for the anti-androgen activity of PCBs. To gain a detailed picture of ligand–protein interaction and identify key residues related with the binding process, the binding free energies of two systems were decomposed on per residues located within 6 Å of the ligand using MM-GBSA approach (Hao et al., 2011a). The electrostatic, van der Waals and total energy contribution of each residue to the PCB-52-AR system and PCB-168-AR system are listed in Supplementary material Table S1. The major favorable energy contributions (almost larger than 1.0 kcal mol1) to the binding free energy were derived predominantly from the residues of Leu704, Leu707, Met742, Met745, and Phe764 in common in both systems. As the above residues were nonpolar, it was reasonable to speculate that these residues can form greater van der Waals interactions with more hydrophobic ligand, and exhibit more favorable nonpolar interaction contribution to the binding free energy. Moreover, special attention was paid to several residues with relatively large different energy contribution. It is observed that the van der Waals energy of Asn705, Val746, Met749, Met780, and Ala877 were mainly responsible for the difference of binding free energy between both systems. Meanwhile, the electrostatic energy of residue Met745 took part in yielding energy difference. The above results showed that hydrophobic contacts and electrostatic interactions played an important role on the binding affinity.

801

The above analysis suggested that the average structure extracted from MD simulations was in good agreement with that of the docked model of the complex. But it should be pointed out that there were some differences between the results of 3D-QSAR model and those of MD simulations, such as the major contribution to the activity. In 3D-QSAR model, electrostatic field was the largest contribution to the activity, but in MD simulations, van der Waals interactions provided significant contribution. This difference may be resulted from the different calculation methods of the parameters and modeling principle. Yang et al. (2009) performed molecular docking and CoMSIA to explore the possible anti-androgenicity of 18 PBDEs. The statistical results of their model (q2 = 0.642 R2 = 0.973) were similar to our study. For molecular docking did not take the impact of solvents into consideration, MD simulation, which consider the influence of solvents and much closer to actual conditions was performed in the present study. MD simulation could not only verify the reliability of molecular docking but also give more detailed ligand–receptor interactions (Liu et al., 2011). Furthermore, the subsequent binding free energy analysis could identify key residues related with the binding process and reveal the contribution of various energy components.

4. Conclusion In this work, a molecular modeling study combining 3D-QSAR, molecular docking, and MD simulation was performed on a series of PCBs to understand the structural characteristic affecting the anti-androgen activity and explore the binding mode between ligands (PCBs) and receptor (AR). The developed receptor-based CoMSIA model yielded high q2 and r2 values. The CoMSIA contour maps visualized the regions of structural features that could explain the variance in the anti-androgen activities. Molecular docking studies revealed that all the PCBs were located at a hydrophobic pocket. MD simulation and binding free energy calculation confirmed the binding models of compound PCB-52-AR and compound PCB-168-AR complexes, and depicted the key interactions. The order of the binding free energies determined by the MM–PBSA method was in good agreement with that of experimental activities. The analysis of the components in binding free energy suggested that the van der Waals interactions provide the substantial driving force for the binding of PCBs and AR. The decomposition of binding free energy to each residue revealed that strong van der Waals interactions with residues Asn705, Val746, Met749, Met780, and Ala877 in the hydrophobic pocket affected the difference of binding and activity between the two compounds. It is anticipated that the results obtained from the computational approach may predict the anti-androgen activity of PCBs accurately and provide a better structural understanding of PCBs-AR interaction mechanism.

Acknowledgements This work was supported by the National Natural Science Foundation of PR China (Grant No. 20737001) and Program for Environment Protection in Jiangsu Province (201140). All authors thank anonymous reviewers and editors for their valuable suggestions on revising and improving the work.

Appendix A. Supplementary material Supplementary data associated with this article can be found, in the online version, at http://dx.doi.org/10.1016/j.chemosphere. 2013.04.022.

802

X. Li et al. / Chemosphere 92 (2013) 795–802

References Agency for Toxic Substances and Disease Registry (ATSDR). 2000. Toxicological Profile for Selected Polychlorinated Biphenyls-update. US Department of Health and Social Services, Public Health Service, Agency for Toxic Substances and Disease Registry, Atlanta, GA. Besselink, H.T., Denison, M.S., Hahn, M.E., Karchner, S.I., Vethaak, A.D., Koeman, J.H., Brouwer, A., 1998. Low inducibility of CYP1A activity by polychlorinated biphenyls (PCBs) in flounder (Platichthys flesus): characterization of the Ah receptor and the role of CYP1A inhibition. Toxicol. Sci. 43, 161–171. Bremle, G., Okla, L., Larsson, P., 1995. Uptake of PCBs in fish in a contaminated river system: bioconcentration factors measured in the field. Environ. Sci. Technol. 29, 2010–2015. Case, D.A., Cheatham, T.E., Darden, T., Gohlke, H., Luo, R., Merz, K.M., Onufriev, A., Simmerling, C., Wang, B., Woods, R.J., 2005. The amber biomolecular simulation programs. J. Comput. Chem. 26, 1668–1688. Cheng, F.X., Shen, J., Yu, Y., Li, W.H., Liu, G.X., Lee, P.W., Tang, Y., 2011. In silico prediction of Tetrahymena pyriformis toxicity for diverse industrial chemicals with substructure pattern recognition and machine learning methods. Chemosphere 82, 1636–1643. Clark, M., Cramer, R.D., Vanopdenbosch, N., 1989. Validation of the general purpose tripos 5.2 force field. J. Comput. Chem. 10, 982–1012. Decastro, B.R., Korrick, S.A., Spengler, J.D., Soto, A.M., 2006. Estrogenic activity of polychlorinated biphenyls present in human tissue and the environment. Environ. Sci. Technol. 40, 2819–2825. Diao, J., Li, Y., Shi, S., Sun, Y., Sun, Y., 2010. QSAR models for predicting toxicity of polychlorinated dibenzo-p-dioxins and dibenzofurans using quantum chemical descriptors. Bull. Environ. Contam. Toxicol. 85, 109–115. Duan, Y., Wu, C., Chowdhury, S., Lee, M.C., Xiong, G.M., Zhang, W., Yang, R., Cieplak, P., Luo, R., Lee, T., Caldwell, J., Wang, J.M., Kollman, P., 2003. A point-charge force field for molecular mechanics simulations of proteins based on condensedphase quantum mechanical calculations. J. Comput. Chem. 24, 1999–2012. Essmann, U., Perera, L., Berkowitz, M.L., Darden, T., Lee, H., Pedersen, L.G., 1995. A smooth particle mesh ewald method. J. Chem. Phys. 103, 8577–8593. Gasteiger, J., Marsili, M., 1980. Iterative partial equalization of orbital electronegativity – a rapid access to atomic charges. Tetrahedron 36, 3219– 3228. Gilbert, K.M., Boos, T.L., Dersch, C.M., Greiner, E., Jacobson, A.E., Lewis, D., Matecka, D., Prisinzano, T.E., Zhang, Y., Rothman, R.B., Rice, K.C., Venanzi, C.A., 2007. DAT/ SERT selectivity of flexible GBR 12909 analogs modeled using 3D-QSAR methods. Bioorgan. Med. Chem. 15, 1146–1159. Hamers, T., Kamstra, J.H., Cenijn, P.H., Pencikova, K., Palkova, L., Simeckova, P., Vondracek, J., Andersson, P.L., Stenberg, M., Machala, M., 2011. In Vitro toxicity profiling of ultrapure non-dioxin-like polychlorinated biphenyl congeners and their relative toxic contribution to PCB mixtures in humans. Toxicol. Sci. 121, 88–100. Hao, G.F., Tan, Y., Yu, N.X., Yang, G.F., 2011a. Structure-activity relationships of diphenyl-ether as protoporphyrinogen oxidase inhibitors: insights from computational simulations. J. Comput. Aided Mol. Des. 25, 213–222. Hao, M., Li, Y., Wang, Y., Yan, Y., Zhang, S., 2011b. Combined 3D-QSAR, molecular docking, and molecular dynamics study on piperazinyl-glutamate-pyridines/ pyrimidines as potent P2Y(12) antagonists for inhibition of platelet aggregation. J. Chem. Inf. Model. 51, 2560–2572. Hevener, K.E., Ball, D.M., Buolamwini, J.K., Lee, R.E., 2008. Quantitative structureactivity relationship studies on nitrofuranyl anti-tubercular agents. Bioorg. Med. Chem. 16, 8042–8053. Jain, A.N., 2007. Surflex-Dock 2.1: robust performance from ligand energetic modeling, ring flexibility, and knowledge-based search. J. Comput. Aided Mol. Des. 21, 281–306. Jakalian, A., Bush, B.L., Jack, D.B., Bayly, C.I., 2000. Fast, efficient generation of highquality atomic charges. AM1-BCC model: I. Methods. J. Comput. Chem. 21, 132– 146. Jorgensen, W.L., Chandrasekhar, J., Madura, J.D., Impey, R.W., Klein, M.L., 1983. Comparison of simple potential functions for simulating liquid water. J. Chem. Phys. 79, 926–935. Kollman, P.A., Massova, I., Reyes, C., Kuhn, B., Huo, S.H., Chong, L., Lee, M., Lee, T., Duan, Y., Wang, W., Donini, O., Cieplak, P., Srinivasan, J., Case, D.A., Cheatham, T.E., 2000. Calculating structures and free energies of complex molecules: combining molecular mechanics and continuum models. Acc. Chem. Res. 33, 889–897. Layton, A.C., Sanseverino, J., Gregory, B.W., Easter, J.P., Sayler, G.S., Schultz, T.W., 2002. In vitro estrogen receptor binding of PCBs: measured activity and detection of hydroxylated metabolites in a recombinant yeast assay. Toxicol. Appl. Pharm. 180, 157–163. Li, F., Li, X.H., Liu, X.L., Zhang, L.B., You, L.P., Zhao, J.M., Wu, H.F., 2011. Docking and 3D-QSAR studies on the Ah receptor binding affinities of polychlorinated

biphenyls (PCBs), dibenzo-p-dioxins (PCDDs) and dibenzofurans (PCDFs). Environ. Toxicol. Pharmacol. 32, 478–485. Li, F., Xie, Q., Li, X.H., Li, N., Chi, P., Chen, J.W., Wang, Z.J., Hao, C., 2010. Hormone activity of hydroxylated polybrominated diphenyl ethers on human thyroid receptor-beta: in vitro and in silico investigations. Environ. Health Perspect. 118, 602–606. Liu, G.Y., Ju, X.L., Cheng, J., Liu, Z.Q., 2010. 3D-QSAR studies of insecticidal anthranilic diamides as ryanodine receptor activators using CoMFA, CoMSIA and DISCOtech. Chemosphere 78, 300–306. Liu, J., Li, Y., Zhang, H.-X., Zhang, S.-W., Yang, L., 2012. Studies of H4R antagonists using 3D-QSAR, molecular docking and molecular dynamics. J. Mol. Model. 18, 991–1001. Liu, J.L., Wang, F.F., Ma, Z., Wang, X., Wang, Y.H., 2011. Structural determination of three different series of compounds as Hsp90 inhibitors using 3D-QSAR modeling, molecular docking and molecular dynamics methods. Int. J. Mol. Sci. 12, 946–970. Mouchlis, V.D., Melagraki, G., Mavromoustakos, T., Kollias, G., Afantitis, A., 2012. Molecular modeling on pyrimidine-urea inhibitors of TNF-alpha production: an integrated approach using a combination of molecular docking, classification techniques, and 3D-QSAR CoMSIA. J. Chem. Inf. Model. 52, 711–723. Nimmanpipug, P., Khampa, C., Lee, V.S., Nangola, S., Tayapiwatana, C., 2011. Identification of amino acid residues of a designed ankyrin repeat protein potentially involved in intermolecular interactions with CD4: analysis by molecular dynamics simulations. J. Mol. Graphics Modell 31, 65–75. Oh, S.M., Ryu, B.T., Lee, S.K., Chung, K.H., 2007. Antiestrogenic potentials of orthoPCB congeners by single or complex exposure. Arch. Pharm. Res. 30, 199–209. Portigal, C.L., Cowell, S.P., Fedoruk, M.N., Butler, C.M., Rennie, P.S., Nelson, C.C., 2002. Polychlorinated biphenyls interfere with androgen-induced transcriptional activation and hormone binding. Toxicol. Appl. Pharm. 179, 185–194. Shen, X.L., Takimoto-Kamimura, M., Wei, J., Gao, Q.Z., 2012. Computer-aided de novo ligand design and docking/molecular dynamics study of Vitamin D receptor agonists. J. Mol. Model. 18, 203–212. Sun, X.Q., Li, Y.Z., Li, W.H., Xu, Z.J., Tang, Y., 2011. Computational investigation of interactions between human H(2) receptor and its agonists. J. Mol. Graphics Modell 29, 693–701. Takigami, H., Suzuki, G., Sakai, S., 2010. Screening of dioxin-like compounds in biocomposts and their materials: chemical analysis and fractionation-directed evaluation of AhR ligand activities using an in vitro bioassay. J. Environ. Monit. 12, 2080–2087. Trossini, G.H.G., Guido, R.V.C., Oliva, G., Ferreira, E.I., Andricopulo, A.D., 2009. Quantitative structure-activity relationships for a series of inhibitors of cruzain from Trypanosoma cruzi: molecular modeling, CoMFA and CoMSIA studies. J. Mol. Graphics Modell 28, 3–11. Wang, J.M., Wolf, R.M., Caldwell, J.W., Kollman, P.A., Case, D.A., 2004. Development and testing of a general amber force field. J. Comput. Chem. 25, 1157–1174. Wang, R.X., Lu, Y.P., Wang, S.M., 2003. Comparative evaluation of 11 scoring functions for molecular docking. J. Med. Chem. 46, 2287–2303. Wold, S., Ruhe, A., Wold, H., Dunn, W.J., 1984. The collinearity problem in linear regression. the partial least squares (PLS) approach to generalized inverses. SIAM J. Sci. Stat. Comput. 5, 735–743. Yang, W.H., Liu, X.H., Liu, H.L., Wu, Y., Giesy, J.P., Yu, H.X., 2010. Molecular docking and comparative molecular similarity indices analysis of estrogenicity of polybrominated diphenyl ethers and their analogues. Environ. Toxicol. Chem. 29, 660–668. Yang, W.H., Mu, Y.S., Giesy, J.P., Zhang, A.Q., Yu, H.X., 2009. Anti-androgen activity of polybrominated diphenyl ethers determined by comparative molecular similarity indices and molecular docking. Chemosphere 75, 1159–1164. Yang, W.H., Shen, S., Mu, L., Yu, H.X., 2011a. Structure–activity relationship study on the binding of PBDEs with thyroxine transport proteins. Environ. Toxicol. Chem. 30, 2431–2439. Yang, Y., Qin, J., Liu, H., Yao, X.J., 2011b. Molecular dynamics simulation, free energy calculation and structure-based 3D-QSAR studies of B-RAF kinase inhibitors. J. Chem. Inf. Model. 51, 680–692. Yang, Y., Shen, Y., Liu, H., Yao, X.J., 2011c. Molecular dynamics simulation and free energy calculation studies of the binding mechanism of allosteric inhibitors with p38 alpha MAP kinase. J. Chem. Inf. Model. 51, 3235–3246. Zhao, Y., Li, W., Zeng, J., Liu, G., Tang, Y., 2008. Insights into the interactions between HIV-1 integrase and human LEDGF/p75 by molecular dynamics simulation and free energy calculation. Proteins: Struct., Funct., Bioinf. 72, 635–645. Zhou, Z.G., Madrid, M., Evanseck, J.D., Madura, J.D., 2005. Effect of a bound nonnucleoside RT inhibitor on the dynamics of wild-type and mutant HIV-1 reverse transcriptase. J. Am. Chem. Soc. 127, 17253–17260. Zhu, H., Ye, L., Richard, A., Golbraikh, A., Wright, F.A., Rusyn, I., Tropsha, A., 2009. A novel two-step hierarchical quantitative structure–activity relationship modeling work flow for predicting acute toxicity of chemicals in rodents. Environ. Health Perspect. 117, 1257–1264.