A study of possible substitutes for the endocrine disruptor DEHP in two hormone receptors

Abstract Bis(2-ethylhexyl) phthalate (DEHP) has been widely used for the production of plastics, and the compound has also been found to act as endocrine disruptor. Exposure to DEHP has been found to cause several hormonal problems, including decreased fertility. Due to the environmental and health risks posed by the use of DEHP, the present study employed molecular docking, molecular dynamics, and free energy analyses (MM-GBSA, MM-PBSA, and SIE) aiming at evaluating the action of DEHP and that of two other compounds (ATEC and DL9TH), tested as potential DEHP substitutes, on two hormone receptors (sex hormone-binding globulin – SHBG – and progesterone receptor – PR). The results obtained showed that ATEC may be a good substitute for DEHP in the production of plastics, such as PVC, considering that the compound recorded the greatest free energy values with respect to binding with SHBG (−31.36 kcal/mol obtained from MM-GBSA; −20.28 kcal/mol for MM-PBSA, and −7.40 for SIE) and PR (−36.40 kcal/mol for MM-GBSA; −27.00 kcal/mol for MM-PBSA, and −8.51 kcal/mol for SIE) – this shows that ATEC presented the least activity in the two hormone receptors. The findings of this study provide relevant insights on potential substitutes for DEHP and help shed light on the action of these new efficient substances, which have similar properties to DEHP (ATEC and DL9TH) yet do not act as endocrine disruptors. Communicated by Ramaswamy H. Sarma


Introduction
The normal functioning of the internal system of the human body is regulated by the endocrine system, which is composed of organs that secrete messenger cells, i.e. hormones. These organs are divided into three groups: (i) endocrine cells; (ii) organs that contain a large part of the endocrine cells (which also play specific roles in other systems); and (iii) organs with some endocrine cells. Hormones are constituted by amino acids and steroid hormones (Marieb & Hoehn, 2007).
The messenger cells act in the human organism by binding hormones to specific receptors (Marieb & Hoehn, 2007), which can be linked to other types of moleculesfor example, endocrine disruptors (ED). These molecules are able to bind to hormone receptors and can cause adverse reactions in the human organism. This class of molecules can be found in drugs, pesticides, and in materials used for the production of plastics, and have become commonly referred to as persistent organic pollutants (POP) (Bila & Dezotti, 2007;Deblonde et al., 2011;Ghiselli & Jardim, 2007;Kabir et al., 2015).
Studies published in the literature have shown that frequent exposure of animals to POPs causes several problems in their organisms; these include decreased egg hatching and reproductive problems (Bila & Dezotti, 2007;Hoffman et al., 2003). In humans, exposure to these molecules, even in low concentrations, can give rise to a wide range of health problems including testicular, prostate, and breast cancers (Bila & Dezotti, 2007;Ghiselli & Jardim, 2007;Johnson & Sumpter, 2001). Bis(2-ethylhexyl)phthalate (DEHP) is an organic compound known to be widely used for the production of plastics, such as PVC (polyvinyl chloride). DEHP has been classified as an endocrine disruptor (ED) and has been linked to a wide range of health problems in living beings. This endocrine disruptor can impair the functioning of steroid organs, and has also been found to inhibit the proper functioning of the reproductive system, leading to decreased fertility, termination of pregnancy, abnormal onset of puberty, endometriosis, and abnormal breast development (Grishkovskaya et al., 2000;Hee-Su et al., 2018;. DEHP can act in hormone receptors, such as progesterone receptor (PR) and sex hormone binding globulin (SHBG), causing adverse effects . In view of the risks posed by the exposure of human organism to DEHP, there has been considerable interest among researchers regarding the discovery and development of molecules that can act as substitutes for DHEP, particularly in the area of plastics manufacturing. Some viable substitutes for DEHP which have been reported in the literature include ATEC (triethyl 2 acetylcitrate) and DL9TH (4-cyclohexene-1,2-dicarboxylic acid dinonyl ester) (Hee-Su et al., 2018;Morishita et al., 2018). Several approaches can be adopted for the study and development of possible substitutes for DEHP; one of such approaches involves the use of computer simulation techniques. The present study employed molecular modeling methods with the aim of studying the action mechanisms of DEHP, ATEC, and DL9TH, and evaluating which of the three molecules has the greatest affinity for hormone receptors. For the insertion of DEHP, ATEC and DL9TH into the hormone receptor binding sites (progesterone receptor and sex hormone-binding globulin), molecular docking was used to help obtain the possible bioactive conformations of these compounds.
To gain a better understanding of the interactions between the three molecules under study and the hormone receptor binding sites, molecular dynamics (MD) simulations were performed. (Adcock & McCammon, 2006;Almeida et al., 2018;Namba et al., 2008) Finally, to help us predict which substance exhibits the greatest activity in the selected hormone receptors and determine which of the two molecules -ATEC and DL9TH can be used as the best possible substitute for DEHP, simulations were conducted aiming at estimating the values of free energy (DG) using Molecular Mechanics -Poisson-Boltzmann Surface Area (MM-PBSA), Molecular Mechanics -Generalized Born Surface Area (MM-GBSA), and Solvated Interaction Energy (SIE) (Genheden & Ryde, 2015;Hou et al., 2011;Lill & Thompson, 2011;Massova & Kollman, 2000;Naïm et al., 2007). The findings of this study will provide relevant insights for the conduct of further experiments targeted at developing new substances with optimized properties that can act as substitutes for DEHP and which do not present any harmful side effects on the endocrine system.

Dataset
For the analysis of the dataset, DEHP and the two possible substitutes for DEHP (namely, ATEC and DL9TH) were selected. The structures (2 D and 3 D) of these molecules are shown in Table 1.
To obtain the 3 D structures of the hormone receptors selected for this study (progesterone receptor -PRand sex hormone-binding globulin -SHBG), a thorough search was conducted in the Protein Data Bank (PDB) (Berman et al., 2000), and the following structures were chosen: 1A28 (PR) and 1D2S (SHBG). The 3 D crystallographic structures of these biological targets are presented in Figure S1 (Grishkovskaya et al., 2000;Williams & Sigler, 1998). The molecular docking method was used to obtain the possible bioactive conformations of DEPH, ATEC, and DL9TH in the receptor binding sites.

Molecular docking
First, the three-dimensional structures of the two biological targets (1A28 and 1D2S) were chosen due to their low resolution values (1.80 Å and 1.55 Å, respectively). Afterwards, the structures were constructed using the Sybyl 8.0 program; in other words, hydrogen atoms were inserted and their atomic charges were calculated by the MMF94 method (Halgren, 1996;Sybyl, 2008). The next step involved the construction of the compound structures from the Avogadro1.2.0 software and calculating the atomic charges using the PM3 semiempirical method from Sybyl 8.0 (Hanwell et al., 2012;Stewart, 2007). It is worth noting that, in the case of the progesterone receptor (PR), the literature shows that there is an important structural water (Verdonk et al., 2005) in the region of the connection site; as such, this was maintained in the simulations.
After preparing the structures, a redocking analysis was performed in order to validate the protocols that were to be used in other molecular docking simulations. Thus, the crystallographic ligands (progesterone and 5-alpha-dihydrotestosterone) were docked in the target binding sites using GOLD 5.2; the protocols employed can be found in the Supporting Information (Tables S1 and S2) (Verdonk et al., 2003). Based on the results obtained, the crystallographic conformations were compared with the new conformations generated from the redocking analysis and the RMSD (Root-Mean-Square Deviation) values calculated from the UCSF Chimera 1.4 were subjected to analysis (Pettersen et al., 2004).
After choosing the redocking parameters with the lowest RMSD values (presented in the Section 3), all other molecular docking simulations were performed using the GOLD 5.2 software. Subsequently, an analysis was conducted aiming at investigating the possible intermolecular interactions between DEHP, ATEC, and DL9TH with the binding sites of the targets. The best conformations (taking into account the higher number of hydrogen bonds (H-bonds) with the main residues mentioned in the literature) of the substance-receptor complex were chosen for the conduct of the molecular dynamics simulations.

Molecular dynamics simulations
The structures of the systems (hormone receptors þ DEHP, ATEC, and DL9TH) evaluated in the molecular docking simulations were used for the conduct of molecular dynamics calculations. To this end, the charges of the three molecules -DEHP, ATEC, and DL9TH, were obtained from the RESP (Restrained Electrostatic Potential) model, and their structures were maintained in the conformation derived from the molecular docking simulations. Single-point energy calculations were only performed using HF (Hartree-Fock) and 6-31 G Ã implemented in the Gaussian 09 software (Ditchfield et al., 1971;Echenique & Alonso, 2007;Frisch et al., 2009). After parameterizing the three substances, the Tleap module (implemented in Ambertools18) was used for the preparation of the systems (D.A. Case, 2018). For the complexes with progesterone hormone (PR), the water model TIP3P was used in a 10 Å octahedral box, while the octahedral box of 12 Å was used for the sex-hormone binding globulin (SHBG).
The box was defined based on inspection in the UCSF Chimera 1.4 (Pettersen et al., 2004). Thereafter, Amber18 was used for performing the molecular dynamics calculations, with the parameters PMEMD (Particle Mesh Ewald Molecular Dynamics) and force field 99SB (ff99SB) (Case et al., 2005;Hornak et al., 2006;Lindorff-Larsen et al., 2010). For all the complexes, five minimizations were carried out with 10,000 steps aiming at removing bad contacts from the structures obtained previously. In the first minimization, the structural movements of the TIP3P water molecules (excluding the structural water in the PR) and the ions were simulated. In the subsequent three minimizations, the structural movements of the receptors were simulated and the DEHP, ATEC, and DL9TH molecules were kept fixed. In the last minimization step, the three substances were simulated using all the degrees of freedom.
After these initial steps, all the complexes were subjected to a thermal bath, with 50 ps variation in temperature (0 to 300 K). The equilibrium was simulated using the following: 10 ns (10,000 ps), 1 atm, Isothermal-Isobaric (NPT) as ensemble, Langevin thermostat (ntt ¼ 3), and the Monte Carlo barostat. After obtaining the equilibrium of the systems, the steps involving the molecular dynamics analysis were carried out using 150 ns and based on the same parameters applied in the equilibrium step. Molecular dynamics simulations were also carried out for the apo form (structures without ligands in the active site) of the two receptors using the same parameters mentioned previously.
All the results obtained were initially analyzed from the CPPTRAJ module (implemented in Ambertools 18); and the RMSD and RMSF (Root Mean Square Fluctuations) plots were obtained. The RMSF graphs indicate the regions of stability. It is worth pointing out that the conformations of the subtrajectories were used to obtain the DG values of the systems under investigation.

Free energy calculations
Based on the most stable conformations of the MD subtrajectories, the MM-GBSA, MM-PBSA and SIE methods were used to calculate the values of DG between DEHP/ATEC/DL9TH and the biological receptors; the results helped assess the degree of interaction of the substances (helped determine which substance presented low interaction) and their binding energy with the receptors. This analysis is very important in the sense that it helps one to determine which molecule is likely to cause the least damage to the endocrine system. For both methodologies, 250 snapshots of the most stable conformations were selected for the calculations of DG.
For the MM-GBSA and MM-PBSA analyses, the simulations were performed using the MMPBSA.py module. Equation (1) shows the parameters used for these calculations: where E MM is the molecular mechanical energy; E El is the energy of electrostatic interactions; E vdW is the energy of van der Waals' interactions; G Pol and G Np are polar and non-polar contributions, respectively; T is the absolute temperature; and S is the entropy (Genheden & Ryde, 2015;Massova & Kollman, 2000). The difference between the methods lies in the calculation of E Elwhile MM-PBSA calculates this energy by the Poisson Boltzman method, the MM-GBSA uses the generalized method of Born to calculate the energy (Genheden & Ryde, 2015).
The SIE method is executed using the SIETRAJ program which calculates the Poisson Boltzman equation by the BEM (Boundary Element Method) method; the application of this technique enhances the efficiency of the method (see Equation (2)).
Based on Equation (2), one can verify the terms used for conducting the SIE methodological analysis: E Coul inter is the Coulomb's intermolecular energy; DG R desolv is the free energy of desolvation; E vdW inter is the van der Waals's energy; DG np desolv is the non-polar free energy of desolvation; q is the derivation factor of Born's atomic radius; Din is the inner dielectric constant of the solute; c is the stress coefficient of the molecular surface; and DMSA is the molecular surface area of the solute on the bond. Based on the three methods mentioned above, free energy calculations were performed with respect to binding between DEHP/ATEC/DL9TH and the biological receptors in order to analyze whether ATEC and DL9TH can be employed as substitutes for DEHP in plastics manufacturing.

Molecular docking
The first step involving the molecular docking simulations is the validation of the parameters from the redocking analyses. Figure S2 shows the two crystallographic ligands that were removed from the PR (PDB 1A28) and SHBG (PDB 1D2S) binding sites. In the redocking analyses, 18 tests were performed with variations in protocol, taking into account rigid and flexible simulations of Ser42 (in the case of SHBG) and Gln726 at the binding site of PR. It is worth mentioning that, for the analyses involving the PR, the presence of a structural water molecule at the binding site was also considered.
All the protocols employed and the RMSD values obtained in the redocking analyses are presented in Tables S1 and S2 (Supporting Information). Based on these results, one can observe that the best protocol for SHBG was the one that employed the CHEMPLP scoring function with the rigid receptor structure and 10 Å cavity as the binding site (RMSD ¼ 0.36 Å). For PR, the best protocol was the one that used the CHEMPLP scoring function with the flexible Gln726 residue, 10 Å of cavity, and a structural water molecule (REMS ¼ 0.47 Å). Thus, these protocols were used for the insertion of DEHP, ATEC, and DL9TH into the hormone receptor binding sites investigated in this study.
All docking simulations were carried out using the GOLD 5.2 software, and several conformations of the three molecules were generated. The best conformation for each molecule was chosen by analyzing the number of hydrogen bonds formed with important residues, as described in the literature (Grishkovskaya et al., 2000;Williams & Sigler, 1998). Figure 1 shows the possible bioactive conformations of the three substances inserted in the binding site of SHBG. Looking at Figure 1, one can observe that the endocrine interferent DEHP forms two hydrogen bonds with Ser128, and ATEC and DL9TH form a hydrogen bond with Ser128. These results show that DEHP possesses a greater ability to interact with SHBG, and this leads to greater binding capacity and greater adverse effects on the organism. With regard to PR, Figure 2 shows the possible bioactive conformations of DEHP, ATEC and DL9TH at the binding site of the receptor.
The intermolecular interactions shown in Figure 2 indicate that the structural water and DEHP at the binding site form two H-bonds with Arg766 and one H-bond with Gln725. This structural water molecule is essentially important when it comes to promoting interactions with the receptor (Verdonk et al., 2005). Along with ATEC and DL9TH, the water molecule at the active site forms two H-bonds with Arg766. Based on this result, one can conclude that DEHP has a greater bond with this receptor; in other words, it possibly has a greater endocrine interfering character.
In short, the results obtained from the docking simulations enabled us to evaluate the possible action mechanism of the three substances at the binding site of the two hormone receptors investigated in this study. Molecular dynamics (MD) simulations were used to assess the dynamic behavior of the substance-receptor complexes.

Molecular dynamics simulations
The possible bioactive conformations of DEHP, ATEC and DL9TH at the active site of the hormone receptors derived from the molecular docking analysis were used for the conduct of MD simulations. After obtaining the trajectories (150 ns) for the six complexes, the CPPTRAJ module was used to obtain the RMSD plots and the interactions between the three molecules and SHBG. The results obtained are shown in Figure 3.
Based on the RMSD plots, Figure 3 shows the conformational changes that ATEC, DEHP and DL9TH caused in the SHBG receptor. One can observe that ATEC presented greater instability in the trajectories, once it exhibited greater variations in the RMSD plot. In the case of DEHP (compound that should be replaced), one can observe the structural changes caused by this substance. The data obtained showed that the presence of DEHP at the SHBG binding site possibly contributed to greater stability of the target; this implies a greater action of DEHP on SHBG, leading to an increase in the binding capacity of this molecule to this receptor, thus causing several adverse effects on the organisms (with DEHP acting as an endocrine disruptor -ED). With regard to DL9TH, the RMSD plot shows that this compound did not cause many conformational changes in the apo form; furthermore, the compound presented a certain degree of instability in the trajectories. A comparison of the action mechanism of the three molecules shows that ATEC presented the greatest instability; this means that the compound does not present a significant action as ED compared to DEHP. Thus, of the three substances under study, ATEC may be the best substitute for DEHP. The results obtained for the interaction of these substances with PR are presented in Figure 4. Figure 4 shows the most stable conformations obtained from the DM trajectories and the RMSD plots. Analyzing the complexes formed with ATEC and DL9TH in relation to the apo form, one can observe that both molecules did not increase the stability of PR. ATEC did not cause many significant conformational changes; this shows that the structure of PR did not have to undergo a significant adaptation in order to accommodate ATEC. With regard to DEHP, one will notice from the RMSD plots that this ED increased the stability of PR; essentially, this points to the fact that the molecule displayed a greater action on PR. In addition, the RMSF plots were also generated for the biological targets (SBGH and PR) in the presence of ATEC, DEHP, and DL9TH. The results obtained for the complexes with SBGH are presented in Figure S3.
In the case of ATEC ( Figure S3), one will notice the occurrence of greater fluctuations in the amino acid residues from the RMSF plots and the presence of superimposed structures. These results show that there was instability in the SBHG structure in the presence of ATEC; this implies that the molecule was unable to fit into the binding site of SBHG; in other words, ATEC can be a good substitute for DEHP, since its action as ED was relatively lower. With regard to DL9TH, in general, the fluctuations were also found to be relatively greater when the receptor was bound to this molecule than in the apo form; however, the fluctuations were still less than those found to occur at the binding site of SBHG upon the complexation of ATEC.
In the case of DEHP, the fluctuations observed in the RMSF plots were found to be generally smaller when SBHG was complexed with this molecule (except for the fluctuations of the last residues which were larger in both cases and were associated with loop residues). This shows that, in the presence of DEHP, SBGH exhibited relatively greater stabilization and less fluctuations in the residue positions; this observation was also corroborated by the results obtained from the RMSD analysis. For PR, the results obtained for the RMSF analysis are displayed in Figure S4.
From Figure S4, one can compare the fluctuations of the complexes formed by PR and ATEC and DL9TH with its apo form. One will notice that, in general, the fluctuations observed for PR were greater in the presence of these molecules. In the case of DEHP, one will observe that the fluctuations were generally smaller or very similar; this shows that some changes occurred in the structure of PR to enable the DEHP to fit better in it. To obtain additional details regarding the action mechanism of ATEC, DEHP and DL9TH on the hormone receptors investigated, free energy calculations were performed using the MM-GBSA, MM-PBSA, and SIE techniques.

Calculations of free energy of binding
After the analysis of the DM simulations from the RMSD and RMSF plots and the generation of the most stable conformations for the substance-receptor complexes, we calculated the DG values for all the systems. To perform these calculations, the following techniques were employed: MM-GBSA, MM-PBSA, and SIE. Table 2 shows the results obtained for the SHBG systems.
From the results shown in Table 2, the DG values obtained for ATEC (À31.36 kcal/molobtained from MM-GBSA; À20.28 kcal/mol for MM-PBSA, and À7.40 kcal/mol for SIE) suggest that this molecule can be a good substitute for DEHP, since its capacity to interact with the hormone receptor is relatively lower. The DG values obtained for DL9TH AND DEHP were found to be close to each other, though the latter exhibited relatively smaller DG values compared to the former (À47.53 kcal/mol for MM-GBSA; À42.16 kcal/mol for MM-PBSA, and À9.37 kcal/mol for SIE); this shows that, of the three compounds investigated, DEHP presented the best affinity with SHBG; this molecule exhibited the greatest action as ED. Table 3 presents the DG values obtained for PR.
From Table 3, one can see that the highest DG values were obtained for PR with ATEC at its binding site. The DG values obtained for the PR þ ATEC complex were: À36.40 kcal/mol (MM-GBSA), À27.00 kcal/mol (MM-PBSA), and À8.51 kcal/mol (SIE). As noted for SHBG, ATEC presented the least capacity of interaction with PR; this shows that this molecule may be a good substitute for DEHP. With regard to DL9TH, the DG values obtained for this molecule were greater than those of DEHP. DEHP exhibited the lowest DG values (À60.91 kcal/mol for MM-GBSA; À44.03 kcal/mol for MM-PBSA, and À10.67 kcal/mol for SIE). The results obtained in this analysis show that DEHP presented the greatest capacity of interaction with PR; as such, exposure to this compound can cause adverse effects in humans and in the environment as a whole.

Conclusions
The present study investigated the possible interactions between DEHP (an endocrine disruptor used in the production of plastics) and two hormone receptorssex hormonebinding globulin (SHBG) and progesterone receptor (PR). The study also evaluated the mechanism of action of two other molecules (ATEC and DL9TH) and their suitability to be used as possible substitutes for DEHP considering their effects on the endocrine system. With the aid of the molecular docking technique, the three substances -DEHP, ATEC, and DL9TH, evaluated in this study were inserted into the binding sites of SHBG and PR; the technique also enabled us to assess the possible interactions that occurred between these three compounds and the two hormone receptors. The application of docking simulations helped obtain possible bioactive conformations for the systems investigated; the results obtained were used for the conduct of molecular dynamics (MD) simulations. These simulations helped evaluate the dynamics of these systems in relation to time. Based on the most stable trajectories, the DG values related to the binding between the two hormonal receptors and the three compounds -ATEC, DEHP and DL9TH, were estimated using Molecular Mechanics -Poisson-Boltzmann Surface Area (MM-PBSA), Molecular Mechanics -Generalized Born Surface Area (MM-GBSA), and Solvated Interaction Energy (SIE) techniques. The results obtained showed that ATEC presented the lowest DG values, indicating that this molecule may be a good substitute for DEHP in plastics manufacturing. The present study enabled us to conduct a physicochemical analysis regarding the action mechanism of the three substances on the hormonal receptors investigated -SHBG and PR. The findings of the study can contribute meaningfully to the search for molecules that can play the same role as DEHP in the production of plastics while not posing serious risks to human health and to the environment.