The Molecular Mechanism of Bisphenol A (BPA) as an Endocrine Disruptor by Interacting with Nuclear Receptors: Insights from Molecular Dynamics (MD) Simulations

Bisphenol A (BPA) can interact with nuclear receptors and affect the normal function of nuclear receptors in very low doses, which causes BPA to be one of the most controversial endocrine disruptors. However, the detailed molecular mechanism about how BPA interferes the normal function of nuclear receptors is still undiscovered. Herein, molecular dynamics simulations were performed to explore the detailed interaction mechanism between BPA with three typical nuclear receptors, including hERα, hERRγ and hPPARγ. The simulation results and calculated binding free energies indicate that BPA can bind to these three nuclear receptors. The binding affinities of BPA were slightly lower than that of E2 to these three receptors. The simulation results proved that the binding process was mainly driven by direct hydrogen bond and hydrophobic interactions. In addition, structural analysis suggested that BPA could interact with these nuclear receptors by mimicking the action of natural hormone and keeping the nuclear receptors in active conformations. The present work provided the structural evidence to recognize BPA as an endocrine disruptor and would be important guidance for seeking safer substitutions of BPA.


Introduction
The item endocrine disruptor chemicals (EDCs) were firstly coined in 1991 at the Wingspread Conference Center in Wisconsin. The paper by Theo Colborn et al. in 1993 was one of the earliest papers about this phenomenon [1]. EDCs are referred to those exogenous substances that can interfere with the endocrine system and then lead to a range of developmental, reproductive, immune, neurological, or metabolic diseases in human and animals [2,3]. Many EDCs are man-made chemicals produced and are released into the environment by industry production such as plasticizers, organotins, pesticides, or alkylphenols [4]. In recent years, many efforts function of its natural agonist of nuclear receptors, the interaction between hERα with estradiol (E2) was also simulated.
Before MD simulations, geometry optimization was performed on BPA and 17β-estradiol (E2) (Fig. 1). The electrostatic potential was calculated using Gaussian 09 program [42] at the Hartree-Fock level with 6-31G Ã basis set. Restrained electrostatic potential (RESP) [43][44][45] was generated using the antechamber module of AMBER10 to describe the partial atomic charges. The LEaP module of AMBER 10 software package was used to add all missing hydrogen atoms of the proteins. Parameters of the ligands and proteins were described by the general AMBER force field (GAFF) [46] and the standard AMBER force field (ff03) [47], respectively. All systems were added appropriate number of sodium counter ions to keep electroneutrality. Each system was solvated using TIP3P water [48] molecules in a octahedron box with at least 10 Å distance around the complex. As a result, the hERα-estradiol system contains 29436 atoms with 8747 waters; the hERα-BPA system contains 32548 atoms with 9515 waters, the hERRγ-BPA contains 30862 atoms with 9049 waters and the hPPARγ-BPA system contains 45521 atoms with 13649 waters. These four complexes were then used as the initial structures for the following molecular dynamics simulations.

Molecular dynamics simulation methods
Molecular dynamics simulations were performed using the AMBER 10.0 software package [49]. For each system, three steps of energy minimization were adopted. Firstly, by restraining the protein with a force constant of 0.5 kcal/mol Å 2 , the solvent molecules and counter ions were minimized by using 2500 steps of steepest descent algorithm first and then switched to 2500 steps of conjugated gradient algorithm. Secondly, by adopting a force constant of 0.5 kcal/mol Å 2 on the heavy atoms of the protein, the minimization was performed using the same method as the first step. Followed this step, the whole system was minimized by using 5000 steps of steepest descent method followed by 5000 steps of conjugated gradient method without any restraint. After the minimization, each system was gradually heated from 0 K to 310 K over a period of 100 ps and maintained at 310 K with a coupling coefficient of 1.0/ps. Finally, a production simulation of 100 ns was performed in the NPT ensemble at a temperature of 310 K and a pressure of 1 atm. During the simulation, periodic boundary conditions were applied and the particle-mesh Ewald (PME) method was used to handle all electrostatic interactions [50] with a cutoff of 10.0 Å using a time step of 2 fs. The SHAKE [51] algorithm was employed on all atoms covalently bonded to hydrogen atoms.

Binding free energy calculations
Binding free energies between the studied ligands and three nuclear receptors were calculated using MM-GBSA method implemented in AMBER 10.0 software package. This method has been used successfully to binding free energy calculations [52][53][54]. 500 snapshots extracted from the last 15 ns MD trajectory of each system were used for MM-GBSA calculations. For each snapshot, the binding free energy is estimated as follows: where G complex , G receptor , and G ligand are the free energy of the complex, receptor and ligand molecules, respectively. Free energy (G) was calculated based on an average over the extracted snapshots from the single MD trajectories. The free energy can be obtained from the molecular mechanics energy E gas , the solvation free energy G sol , and the solute entropy S, respectively [55].
G ¼ E gas þG sol À TS ð2Þ Where E gas is the gas-phase energy; E int is the internal energy; E ele and E vdw are the Coulomb and van der Waals energies, respectively. E gas was calculated using the AMBER03 molecular mechanics force field. G sol is the solvation free energy and can be decomposed into polar and nonpolar contributions. G GB is the polar solvation contribution calculated by solving the GB equation. Dielectric constants for solute and solvent were set to 1 and 80, respectively. G nonpolar is the nonpolar solvation contribution and was estimated by the solvent accessible surface area (SAS) determined using a water probe radius of 1.4 Å. The surface tension constant γ was set to 0.0072 kcal/mol/Å 2 [56].
In order to obtain the detailed interaction of two ligands with nuclear receptors, MM-GBSA calculated binding free energy were decomposed to the contribution of each residue by only considering molecular mechanics energies and solvation energies without considering entropy.

Stability of the simulation systems
To monitor the equilibration of each system, the root-mean-square-deviation (RMSD) of protein backbone atoms, protein Cα atoms, Cα atoms of the residues within 5 Å around the ligand and the heavy atoms of BPA and E2 with respect to the initial structures were monitored and shown in Fig  2.2Å and ultimately attained equilibrium at about 80ns. However, for of hPPARγ-BPA and hERα-E2 systems, the RMSD values of protein backbone increased sharply at the first 15ns to about 3.0Å and then fluctuated around 2.5Å~3.5Å at the rest time, indicating the conformation change of the proteins. As shown in Fig. 2c, the RMSD values of the active sites for hERα-BPA and hERα-E2 were less than 1 Å and were stable through the simulations, indicating no significant structural fluctuations of the active sites. However, for hERRγ-BPA and hPPARγ-BPA, the evolution of RMSD values was not so stable. For the complex of hPPARγ-BPA, it is easy to understand since hPPARγ has a much larger pocket and the complex of hPPARγ-BPA was obtained by molecular docking. In the complex of hERRγ-BPA, RMSD values of the active sites were stable around 0.5~0.8 Å until the significant increase to above 1.0 Å at 60 ns, which was mainly caused by the conformational change of H6, β1, β2 and the link loops adjacent to them (residues 317-341) ( Fig. 3b and 3d). By visual inspection of the conformation change, overturn of one ring of BPA to the opposite direct induced the H6 conformation change. This was in accordance with the RMSD evolution of BPA in hERRγ (Fig. 2d), which revealed the multi-conformations of BPA in the binding pocket (Fig. 3c). For the ligands (Fig. 2d), the RMSD evolution in the four complexes had quite different trends. E2 was always stable in the pocket throughout the simulation, suggesting that E2 could steadily bind with hERα. However, BPA in the three binding pockets had much larger fluctuations with quite different trends. In spite of this, BPA could achieve stable binding pose in three different pockets.
In addition, the root-mean-square-fluctuations (RMSF) of Cα atoms of the proteins were also calculated and drawn in Fig. 4. As indicated in Fig. 4a, the RMSF of hERα in complex with BPA and E2 had similar values, indicating BPA and E2 bound to hERα in the same manner and interacted with hERα in a similar mechanism. When interacting with BPA, residues in hERRγ and hPPARγ ( Fig. 4b and 4c) had small RMSF values. As a whole, most of the residues in four complexes were quite stable with minor mobility, except for some flexible regions that locate in the loops or two ends of proteins.

Active conformations of the nuclear receptors
When binding with the agonists or activators, nuclear receptors were activated by switching from inactive state (Fig. 5b) to active state with the twelfth helix (H12) sealed the hydrophobic binding pocket constituting with H2, H5, H7, H11, β1 and β2 (Fig. 5a). As an activator, E2 had the ability to stabilize hERα in the active state with H12 covering the binding pocket through the simulations and then form a shallow hydrophobic groove (helices H3-H5 and H12) for peptide coactivators binding. The receptors were then activated (Fig. 5a) and triggered downstream signaling pathways. In contrast to the active conformation, H12 was shifted away from the ligand binding pocket to occupy the shallow groove of coactivators in the inactive conformation of the receptors (Fig. 5b). When taken up by H12, coactivators could no longer bind to the groove and the downstream signaling pathway was obstructed.
The binding of BPA with hERα, hERRγ and hPPARγ could also hold H12 (As for hERRγ and hPPARγ, it refers to the helix with the same position as in hERα.) in the similar position as that in hERα-E2 complex so as to maintain the receptors activated and facilitate the binding of coactivators.
In order to further inspect the stability of active conformations of four complexes, the distances between H12 and H11 were monitored and depicted in Fig. 6. In spite of large fluctuations at 30-40 ns of hERα-E2, there were no structural clashes and all the four complexes were able to keep the active state of proteins through the simulations with fluctuations less than 4 Å. It indicated that BPA could stably bind and interact with active nuclear receptors and keep the receptors activated. In order to gain insight into the contribution spectrum of binding energy for E2 to hERα, BPA to hERα, hERRγ and hPPARγ, the enthalpy contributions during the last 15ns of simulation were calculated for each system using the MM-GBSA method and were shown in Fig. 7. Due to long time-consuming, the entropy contributions were not calculated. The obtained binding free energies were shown in Table 1. It can be seen from Fig. 8 and Table 1, as a natural estrogen, E2 bound to the receptor with the strongest affinity with a binding free energy of -34.88 kcal/mol. Though with lower binding affinities, BPA binding to hERα, hERRγ and hPPARγ still had moderate potencies, with the binding free energies of -23.77 kcal/mol, -26.69 kcal/mol and -23.39 kcal/mol, respectively. According to the order of binding free energies, the binding potencies of BPA towards hERα, hERRγ, hPPARγ can be ranked as: hERRγ > hERα > hPPARγ, which was almost consistent with the experiment results [11,17,36]. From the contribution of various energy components of binding free energies (Table 1), the nonpolar interaction (ΔG nonpolar = ΔE vdw + ΔG sol-np ) contribution is the main driving force for ligand binding. Contrary to the hydrophobic contribution, the polar interaction (ΔG polar = ΔE ele + ΔG sol-ele )   contributed unfavorably to the binding of ligands. In fact, the direct intermolecular electrostatic interactions (ΔE ele ) made favorable contribution to the binding process, while they were compensated by the large desolvation penalties (ΔG sol-ele ) ( Table 1).

Identification of the key residues responsible for the binding of BPA
In order to obtain a detailed profile about the binding process, the total binding free energy was further decomposed to each residue by using MM-GBSA method. The corresponding results were depicted in Fig. 8. For the binding of E2 and BPA, several key residues were  identified in three receptors. For example, in hERα, the hotspot residues include M343, L346, A350, L384, L387, L391 F404, M421, H524 and L525. In hERRγ, residues L268, L271, A272, L309, V313, Y326, L345 and F435 contribute obviously to the binding of BPA. As for hPPARγ, the main hotspot residues are I281, F282, C285, I326, Y327, L330 V339, L353 and M364. The identified key residues were in well consistent with previous studies [4,11,19]. Residues with energy contribution larger than 0.5 kcal/mol were indicated for each complex in Fig. 8. The polar and nonpoalr contributions of each identified key residues were given in Fig. 9. From Fig. 7a and Fig. 7b, it can be seen that for the binding of E2 and BPA to hERα, the hotspot residues are almost the same, suggesting that BPA can interact with hERα in a similar way to E2. From Fig. 7c and Fig. 7d, the identified residues of hERRγ and hPPARγ are mostly hydrophobic. Furthermore, as shown in Fig. 8 and Fig. 9, the nonpolar interactions of a few hydrophobic residues make main contribution to the binding, indicating the hydrophobic interactions are important to the binding process of BPA to several targets. In addition, residues R394, H524 from hERα, L271, E275, V313, R316 from hERRγ and F282, I326, M329, V339 from hPPARγ also have obvious polar contributions to E2 and BPA binding.

The comparison of binding mode of E2 and BPA to hERα
In order to understand the detailed mechanism about how BPA affects the function of nuclear receptors, a 100 ns MD simulations of E2-hERα complex was performed as a comparison because E2 is a strong endogenous estrogen receptor agonist. In Fig. 10a, the core ring section of E2 formed strong interaction with several hydrophobic residues of hERα, such as Leu346, Leu384, Leu387, Met388, Phe404, Met421 and Leu525. The T-shape stacking between the phenyl group of E2 and the phenyl group of Phe404 could still be favorable to the ligand binding. Moreover, the hydrogen bonds formed between the ring A and Glu353 and Arg394 as well as the hydrogen bonds between the ring D and His524 further strengthened the ligand-receptor interaction throughout the simulation. Compared to E2, BPA binding to hERα adopted a fashion similar to E2 with two phenol fragments pointing to two ends of the hydrophobic pocket, one ring to residues Glu353 and the other to His524 (Fig. 10b). Even so, the hydrogen bond interactions formed between BPA and Glu353, Arg394 and His524 greatly reduced compared with the coresponding hydrogen bond interaction for E2 ( Table 2). The nonpoalr component of the binding free energy contributing to BPA binding mainly came from the van der Waals contibution of residues Met343, Leu346, Leu387, Leu391, Phe404, Met421 and Leu525. However, BPA has a smaller structure than that of E2 and occupies less space in the binding pocket. Thus, the interaction between BPA and several residues was impaired to a certain degree, such as the van der Waals contribution of M343, Leu384, Met388, I424, L428 and Leu525. The electrostatic contirbution comes from Met343, Glu353, Arg394 and His524 ( Fig. 8 and Fig. 9), leading that the direct interaction between BPA and hERα weakened. Structural basis of BPA binding to hERRγ and hPPARγ Fig. 11 shows the direct interaction between BPA and hERRγ. As can be seen, BPA binds to hERRγ with two nonplanar phenol fragments pointing toward two ends of the binding pocket, wchi is similar to that of E2 and BPA in the pocket of hERα. In the initial structure, BPA   Table 2, the occupancy is about 70%). As the ring B of BPA rotated away from Asn346 and Tyr326 to a more stable position, the hydrogen bond network between these two residues were interrupted ( Fig. 3c and Table 2). Herein, the hydrogen bond between BPA and Glu275 played key role in the ligand-receptor interaction. In addition, most residues of the binding pocket surrounded BPA were hydrophobic. These hydrophobic residues such as Leu268, Leu271, Ala272, Leu309, Ile310, Val313, and Phe435, form strong hydrophobic interactions with BPA. Besides, the face-to-face interaction formed between the B ring of BPA and Phe435 was also strong (Fig. 11). The interaction between BPA and hPPARγ was given in Fig. 12. As the initial complex structure was obtained using docking method, BPA in the binding pocket moved a lot to have a more appropriate location (Fig. 12). Although BPA is moved in the pocket, the two rings of BPA extended toward two directions, similar to that in hERα and hERRγ. BPA has van der Waals and electrostatic interactions with residues Ile281, Phe282, Cys285, Ile326, Leu330, Val339, Met364 and Tyr327 of hPPARγ (Fig. 12c). Because the buried binding pocket in hPPARγ has a volume larger than that in hERα and hERRγ [36], the direct hydrogen bond interactions were crippled. In spite of this, the van der Waals and electrostatic contributions could still favor the binding and trap BPA in the binding pocket.
Overall, the obtained results show that BPA can bind and interact stably with three nuclear receptors and keep the twelfth alpha helix sealing the hydrophobic ligand binding pocket. With the twelfth alpha helix covering the pocket, the receptors were trapped in an active conformation which allowed the binding of coactivators and would transduct the downstream gene transcription signaling. Once activated abnormally by exogenous chemicals, the normal regulation of endocrine system was disorganized and normal functions were disordered.

Conclusions
In this study, the molecular mechanism of BPA binding to hERα, hERRγ and hPPARγ was investigated by using molecular dynamics simulations and MM-GBSA calculations. The simulation results demonstrate that BPA can bind to the hydrophobic pocket of the three studied nuclear receptors. The twelfth helix (H12) of the binding pocket was in activated state when BPA binding to receptor. BPA also has stable interaction with the nuclear receptors by mimicking the behavior of natural hormone E2. The calculated binding free energies are favorable to BPA binding and the binding process is mainly driven by van der Waals and hydrogen bond interactions. The calculated binding free energy of BPA to hERRγ is the lowest. The binding free energy to hPPARγ is the highest. Binding mode analysis suggests that BPA can stay in the pocket with two rings extended to interact with the residues in the pocket. The obtained results can provide structural evidence of BPA as an endocrine disruptor and will be great importance in the guidance of searching for safer BPA substitute.