Molecular docking for screening chemicals of environmental health concern: insight from a case study on bisphenols

Abstract

To explore the use of molecular docking as a high throughput in silico screening tool for identifying chemicals of environmental health concern, we conducted a case study to assess endocrine disruption effects due to targeting of nuclear receptors (NRs) by chemicals with backbone structures like bisphenols, but with varied functional groups. The molecular docking analysis elucidates how functional groups of the chemicals, such as NH₂, Cl, and OCH₃, influence their interaction with the human estrogen receptor alpha (hERα), a key player in endocrine regulation. Through comparative docking analysis, we examined how bisphenol analogs interact with three distinct conformations of hERα: the apo structure and two structures with bound agonist and antagonist ligands. Water molecules within the protein and surrounding the ligand binding domain (LBD) were found to have little impact on the affinity of compounds binding to the receptor across various conformations. This can be attributed to the hydrophobic nature of the ligand-binding pocket, which consists mainly of hydrophobic amino acid residues and binding sites. In the assessment of bisphenol analogs compared to well established endocrine disrupting chemicals (EDCs), it was observed that these analogs exhibit characteristics commonly associated with endocrine disruptors. While compounds like BPA and BPF exhibited partial agonist activity, stimulating hERα activity to varying degrees, other compounds displayed non-agonist behavior, suggesting a different mode of interaction with the receptor. Further analysis revealed the significance of specific functional groups, such as hydroxyl or amine groups, on the aromatic ring of these compounds in modulating their binding affinity to hERα. Within the ligand binding site of hERα, amino acid residues Glu353, Arg394, and His524 have the capacity to form hydrogen bonds with hydroxyl or amine groups. Protonation or deprotonation of these groups can further alter their binding affinity, thereby influencing their interaction with estrogen receptors and subsequent estrogenic effects. Via this case study, we demonstrate the potential and provide best practices of using molecular docking as a new approach methodology (NAM) for chemical assessments and regulations.

---

Environmental significance

Environmental processes and impact assessment of a large number of synthetic chemicals needs high throughput screening and prioritization. Various in silico models have been developed and used as cost-effective approaches to prioritize chemicals for environmental health risk assessment but there is a need to extend external exposures to chemical binding to specific receptors causing toxic effects within the body. This study explores the application of molecular docking as a New Approach Methodology for chemical assessment. Binding of chemicals with similar backbone structures and different functional groups to an estrogen receptor linked to endocrine-disrupting effects was explored. The influences of receptor conformations, ligand conformers, protonation/deprotonation states, and inclusion/exclusion solvent molecules within the binding pocket of a receptor on the ligand–receptor binding are explored and necessary considerations for the applications of molecular docking are highlighted. Molecular docking as a mechanistic based model that captures physicochemical interactions between ligands and receptors has its advantage in screening chemicals in structures because it is independent of data used for training other in silico models.

Introduction

An increasing number and larger quantities of chemicals have been produced and used in industry, agriculture, and consumer products around the globe.¹ Comprehensive environmental and health risk assessment is not feasible for all these chemicals due to the high cost associated with obtaining experimental data, essential for risk characterization of each chemical.² Although regulatory frameworks such as REACH have started to address the safety evaluation of chemicals produced before 1981, which account for over 97% of major chemicals in use, there is still a substantial data gap, with toxicity information missing for about 86% of these substances.³ Traditional experimental methods, such as in vivo and in vitro testing, are often too slow, costly, and ethically constrained to keep up with the large number of new chemicals being introduced into the environment. Furthermore, these tests do not always predict how chemicals will affect humans very well.^3–5 Because of the large number of chemicals that need to be tested, it is essential to develop faster and robust methods that can screen many chemicals for their safety.³ These considerations together have resulted in a shift from in vivo methods to in vitro or in silico approaches for chemicals' risk assessment.⁶

In silico approaches based on quantitative structure–activity relationships (QSARs),⁷ environmental fate and human exposure modeling,⁸ and machine learning⁹ have proved effective in screening large numbers of chemicals and prioritizing them for more comprehensive risk assessment.¹⁰ While process-based mechanistic models have been effectively used for screening the environmental persistence, long-range transport and bioaccumulation¹¹ of many chemicals, in silico toxicity screening for a large number of chemicals remains challenging because available QSARs developed based on experimental data are not available to capture different modes of action.¹² Available computational approaches such as molecular dynamic simulations allow modelling of the role of xenobiotics in biochemical processes,¹³ but such approaches are often computationally too intensive to be applied to large numbers of chemicals over a reasonable amount of time. For identification purposes, a high throughput method is preferred.¹⁴

Molecular docking is a powerful high-throughput in silico approach widely used in drug discovery and environmental toxicology to predict how small molecules, including environmental contaminants, fit and interact with biological targets such as proteins and receptors based on molecular mechanics and force fields.^15,16 Within environmental toxicology, molecular docking plays a key role in clarifying the molecular mechanisms of toxicity by demonstrating how contaminants interact and disrupt the structure and function of vital proteins that regulate physiological functions.^17–19 By calculating different binding poses that show the orientation of a small molecule within the active site of a protein, molecular docking provides information about binding affinity, specificity, and the potential for disrupting normal protein function. These insights are essential for understanding how environmental chemicals may trigger adverse biological outcomes, such as endocrine disruption or enzyme inhibition, even before in vitro or in vivo testing.²⁰ For example, endocrine disruptors can bind to estrogen receptors by occupying hormone-binding sites and prevent normal hormone–receptor interactions and hormonal signaling pathways.^21,22 Additionally, the rapid screening capability of molecular docking enables the evaluation of large libraries of environmental chemicals, supporting early hazard identification and prioritization for further investigation of contaminants against relevant biological targets.^17,18,23 Among these, endocrine-related proteins involved in hormone regulation and susceptible to disruption by environmental chemicals are frequently studied, with particular attention to the estrogen receptor (ER).^21,23–26

ERs are functional proteins belonging to the nuclear receptor (NR) superfamily. These proteins serve as ligand-activated transcription factors and are essential to numerous biological processes, such as lipid metabolism, embryonic development, homeostasis, and cell death.^27–31 They typically consist of three functional domains; the N-terminal transactivation domain, the central DNA binding domain, and the C-terminal ligand binding domain (LBD), which contains the activation function-2 (AF-2).^32–34 The activity of the human estrogen-related receptor (hERR) is significantly influenced by the conformation of the LBD, which is dependent on ligand binding. X-ray crystallography of the nuclear receptors has revealed three primary structural conformations of the LBD – apo, agonist and antagonist bound conformations (Fig. S1) – with distinct positions of the 12 helix determined by the biological properties of the ligand, whether it functions as an agonist or antagonist.^25,35–38 Because the structural conformation of the LBD controls receptor activity, any compound capable of binding to this domain has the potential to alter normal receptor signaling.^35,37 This property makes estrogen receptors a key target for endocrine-disrupting chemicals (EDCs), which can mimic or inhibit endogenous hormones and thereby modulate physiological pathways.^24,29 Among endocrine-disrupting chemicals, bisphenol A (BPA) has emerged as a prominent compound of concern, with significant impact on human health, causing e.g. diabetes, obesity, cardiovascular diseases, reproductive disorders, and cancer.^39–43 For decades, bisphenol compounds exemplified by BPA have been pervasive in plastics widely utilized for food storage purposes. Bisphenol polymers can release monomers into food and drinks, potentially disrupting endocrine pathways by acting like estrogen.^44,45 As a result, numerous countries have imposed restrictions or completely banned the use of BPA in food-related applications.^46,47 In response to these concerns, a new generation of materials marketed as “BPA-free” has emerged. However, many of these alternative compounds, such as bisphenol F (BPF) and bisphenol S (BPS), may still possess endocrine-disrupting capabilities by interacting with ERs.^46,48–50

Besides BPA alternatives, there are many other industrial chemicals and their transformation products that share similar backbone chemical structures but with different functional groups. These chemicals are likely to act as endocrine-disrupting compounds akin to BPA.^48,50 For example, chlorinated BPA has been reported as an endocrine disruptor.^51,52 Evidence also suggests that hydroxyl polychlorinated biphenyls⁵³ and hydroxyl polycyclic aromatic compounds^54,55 as environmental transformation products are also likely to interact with ERs and cause endocrine disrupting effects. As a result, hERRs are considered as one potential target of chemicals that possibly act as EDCs.^25,27,28,56 Given the large number of chemicals produced and used, a high-throughput in silico screening method is needed to identify chemicals with potential risk of endocrine disruption by targeting NRs.²⁶ In this study, bisphenol A, bisphenol F, and bisphenol S were selected due to their widespread use, well-documented endocrine-disrupting effects, and extensive existing literature, which provides a solid foundation for molecular docking studies.^{40–45,48,49,57–62}

This article highlights the potential and versatility of molecular docking as a high-throughput in silico screening tool for fast assessment of chemicals with potential NR-related endocrine disruption effects. We employed three distinct conformations of the estrogen receptor alpha (ERα) to model NRs. As a case study, we selected thirty chemicals with structures resembling bisphenol A (BPA), bisphenol F (BPF), and bisphenol S (BPS) and 40 well-known EDCs⁶³ as benchmarks for the screening. Through this case study, we examine the various factors that cause variability in the predicted docking scores and how the results can be interpreted to identify chemicals with high binding affinities to proteins, and thus implications for human population and environmental health.

Method

Study design

Fig. 1 shows the workflow of using molecular docking to screen chemicals with binding affinities to health impact related protein. The Open-Eye software was used to conduct molecular docking simulations for assessment of bisphenol analogs. Receptor preparation involved two versions: one with water molecules removed and another with water molecules preserved to assess the impact of solvent presence on docking results. Ligands were geometry-optimized, protein-docked, and the resulting binding affinity was evaluated. Post-docking analysis involved the influence of water molecules on binding interactions, the effects of K_OW values, and protonation/deprotonation states of ligands on binding predictions. Additionally, we applied the same molecular docking workflow to known EDCs and use the results as benchmarks to screen our targeted chemicals (BPA analogs). This benchmark based approach helps eliminate systematic error possibly originated from assumptions and simplifications of molecular docking.

Fig. 1 Study design of using the molecular docking method to characterize interactions between ligands (chemical contaminants) and receptors (protein related to health impact).

Ligand selection and preparation

The thirty selected compounds with structural frameworks similar to bisphenol A (BPA), bisphenol F (BPF), and bisphenol S (BPS) are listed in Table S1. These compounds share the backbone structure of bisphenols, with different functional groups such as NH₂, Cl, and OCH₃ as OH substituents. Additionally, 40 well-known EDCs⁶³ with the same functional groups on the aromatic ring and without long carbon chains (Table S2) were selected as benchmarks for comparison and evaluation of the estrogenic binding affinity of the investigated compounds. The .sdf files representing chemical structures of the ligands were retrieved from PubChem. In instances where structures of the molecules were not accessible via PubChem, they were drawn using ACD ChemSketch Freeware Version.⁶⁴ Both neutral and protonated/deprotonated forms of the chemicals were included to account for the different ionization states and corresponding interactions with the receptors. The 3D geometry of all ligands was optimized using SZYBKI 2.7.0.3 with default approaches and parameters – the MMFF94 force field and the Sheffield solvation model with AM1BCC solution charges. This approach adds an additional electrostatic term to accurately model the solution environment.⁶⁵ Conformers for each ligand were generated using OMEGA 5.0.0.3,^66,67 a tool for generating a diverse set of 3D conformers of molecules. Estradiol (E2) was employed as a positive control for both the apo and agonist conformers due to its well-established role as a potent agonist of estrogen receptors. In fact, E2 is chosen as the positive control because it is the native ligand that naturally binds to estrogen receptors, providing a biologically relevant and standardized reference point for evaluating the effects of experimental compounds on estrogen receptor activity.⁶⁸ On the other hand, 4-hydroxytamoxifen (4-OHT) was selected as the positive control specifically for the antagonist conformer. 4-OHT is a well-known selective estrogen receptor modulator (SERM) that acts as an antagonist by blocking the binding of estrogen and inhibiting receptor activation.^38,69

Protein preparation

Precise evaluation of receptor–ligand binding affinities through molecular docking calculations requires the selection of suitable template structures.⁷⁰ It is often impractical to select a single ideal receptor structure for in silico investigations, and focusing on molecular docking studies with a single receptor structure may not be suitable for assessing the binding affinities of diverse chemicals. In such scenarios, employing statistical analyses of binding affinity values obtained from molecular docking calculations with multiple receptor structures could effectively minimize prediction errors.²⁵ Hence, the protein receptor structures selected in this work includes the apo conformation of hERR (complexed to E2), the agonist conformation with bound E2, and the antagonist conformation with bound 4-OHT. These structures represent different functional forms of the receptor, and allow a comprehensive investigation of how different ligands interact with the hERR LBD. The binding of an agonist induces conformational change in the receptor, transitioning it from its inactive form to an active conformation in which the receptor is capable of initiating a signal transduction pathway or cellular response. On the other hand, the receptor antagonist conformation is stabilized upon binding of an antagonist, causing inactivation or blocking its transition to an activated form.⁷¹ The protein structures of the apo conformation (1A52),³⁶ the agonist conformation (1GWR),³⁷ and the antagonist conformation (3ERT)³⁸ were retrieved from the Protein Data Bank (PDB). The receptors were prepared using the Spruce 1.5.3.3 (ref. 72) module of the Open Eye Scientific Software package, which corrects for structural gaps, such as incomplete side chains and loops, refinement of hydrogen atom positions and assessment of various tautomer states of the ligands and cofactors to ensure compatibility within the biomolecule structure. The binding pocket size was determined based on a default distance of 5.0 Å from the ligand, and maximum limit for processing atoms in the entire system was set to 50 000. The ligand binding pockets of all receptors were generated around the structure of active ligands and all receptors were visualized by the Make Receptor software. The ligand binding pockets generated for the 1A52, 1GWR, and 3ERT receptor conformations have volumes of 2602 ų, 2444 ų, and 8967 ų, respectively. These volumes represent the confined spaces within which interactions between the receptors and other molecules can be studied or simulated. To assess the impact of solvent molecules in the LBD on molecular docking results, the receptors were prepared in two different ways: one with water molecules removed and another with water molecules preserved. This approach allows for a detailed investigation of how the inclusion or exclusion of water affects the receptor–ligand binding interactions and overall molecular docking performance.

Molecular docking

The docking simulations were conducted using FRED 4.3.0.3.^73–75 FRED carries out docking of different conformations of a molecule into a single receptor through an exhaustive search method. It methodically examines all possible rotations and translations of each ligand conformation within the receptor's active site. Chemgauss3 (Ch3) served as the default exhaustive scoring function in FRED. The Chemgauss3 scoring function assesses how well ligand poses fit into the active site using Gaussian smoothed potentials. It takes into account several types of interactions, including: shape complementarity, hydrogen bonding, and hydrogen bonds with implicit solvent and metal–chelator interactions. In the docking calculations, binding energies are reported such that negative values indicate stabilizing interactions, with more negative energies corresponding to stronger predicted binding. This convention clarifies that compounds with the lowest (most negative) energies are expected to have the strongest interactions with the target receptor. The standard approach, using a 1.0 Å translational step size and a 1.5 Å rotational step size, was employed to manage the resolution of docking throughout both the exhaustive search and optimization phases.⁷⁴ VIDA 5.0.5.3 (ref. 76) was employed for the visualization of structures, analysis of molecular surfaces, and examination of protein–ligand interactions. To investigate the influence of solvent on the docking results, we carried out docking calculations using the same LBD, with and without water molecules. Water molecules play a crucial role as spacers or mediators of hydrogen bonds between ligands and residues within the ligand binding pocket (LBP).⁷⁷ To evaluate the accuracy of docking protocol, docked ligands, including the positive controls estradiol (E2) and 4-hydroxytamoxifen (4-OHT), were aligned to their corresponding crystal structures using PyMOL, and RMSD values were calculated with rms_cur command. Only heavy atoms were included, excluding hydrogens, to avoid bias from differences in protonation.

Conformational preference of ligands

Ligands were docked to all three hERR protein conformations. The ligand binding preference for one conformation or the other could significantly influence their biological activity, i.e. whether they behave as agonist or antagonist. To establish a quantitative measure for discriminating between agonist and antagonist behavior, we adopted the conformational preferences (Cpf) factor, defined as:²⁴

Cpf = ΔG (agonist conformation) − ΔG (antagonist conformation) (1)

where ΔG (agonist conformation) is the binding affinity obtained from molecular docking calculations using the hERR agonist conformation (1GWR), while ΔG (antagonist conformation) is that derived using the hERR antagonist conformation (3ERT). This parameter can be used to determine the specificity of ligands towards hERR, whether they tend to behave as agonists or antagonists of the receptor.

Results and discussion

Docking accuracy and validation

To assess the reliability of our docking protocol, we first validated it with well-known reference ligands. Estradiol (E2) and 4-hydroxytamoxifen (4-OHT) were docked into the receptor's ligand-binding domain, and their RMSD values were calculated. The RMSD analysis of these positive controls indicated that the docked ligands matched their crystal structures very closely. Estradiol (E2) showed an RMSD of 0.7 Å, and 4-hydroxytamoxifen (4-OHT) had an RMSD of 1.1 Å (Fig. S2). These low values show that the docking method can reproduce the experimentally observed binding poses. The close agreement also indicates that the receptor/ligand preparation methods and docking settings used in this study were appropriate, which gives confidence in the docking and structure–activity analysis of bisphenol analogs.⁷⁸

Effects of water molecules within the LBP on ligand docking calculations

Comparison of molecular docking results with and without water molecules reveal that the presence of solvent within the protein and LBP has negligible influence on the binding affinity of the compounds (Tables S3–S5). The maximum difference in binding affinity between the results obtained with water molecules and those without water molecules was found to be 1.06 kcal mol⁻¹, which is obtained for 4-[2-(4-hydroxyphenyl)propan-2-yl]phenol docked in the apo conformer (Fig. S3). The average difference and the root mean square deviation (RMSD) between the binding affinities with and without water were 0.19 and 0.24 kcal mol⁻¹, respectively, which are relatively small compared to the large deviation observed between bisphenol analogs and the positive control E2, which showed an RMSD of 5.34 kcal mol⁻¹ and absolute differences ranging from 4.38 to 7.2 kcal mol⁻¹. This suggests that water molecules do not dramatically affect binding interactions. Similar findings have been reported in another study.²⁵ This observation can be related to the predominantly hydrophobic nature of the LBP, which is mainly made up of multiple hydrophobic amino acid residues.⁷⁹

Also, the molecular docking analysis indicates that the binding affinity of compounds is influenced by several factors associated with the binding site's characteristics. These factors include shape energy, hydrogen bond energy, protein desolvation, and ligand desolvation energy. Shape refers to the volume of a molecule, and in Chemgauss, shape energy is related to the ligand's spatial arrangement and interactions of its heavy atoms within a protein's active site.^75,80 As detailed in Fig. S4, shape energy of all ligands accounts for more than 80% of the total binding energy, while hydrogen bond energy contributes around 20%. This major contribution of shape energy indicates that the main determinant of binding strength is the spatial alignment between the ligand and the receptor, which is influenced by the properties and arrangement of the residues within the ligand-binding pocket. These findings show that the inclusion of water molecules does not impact the binding affinity of the ligands, which highlights the robustness of the spatial complementarity and the hydrophobic interactions within the ligand-binding pocket (Fig. 2). Thus, for the following analysis, water molecules in the LBP crystal structures were retained.

Fig. 2 Comparison of ligand binding affinities with and without water molecules towards three ER conformations: (A) agonist conformation (PDB ID: 1GWR), (B) apo conformation (PDB ID: 1A52), (C) antagonist conformation (PDB ID: 3ERT).

Influence of ligand conformers on binding affinity

We examined the impact of conformational analysis of the compounds on binding affinity toward the receptors. The maximum difference in binding affinity observed between various conformers of a ligand was −2.44 kcal mol⁻¹, as illustrated by 4-[(4-aminophenyl) methyl]aniline ligand. This variance in binding affinity is a result of multiple factors that influence how the ligand interacts with the receptor. As previously discussed, these factors include shape energy, hydrogen bond energy, protein desolvation energy, and ligand desolvation energy. For compounds capable of forming hydrogen bonds, the difference in binding affinity (ΔG) between the conformers with the highest and lowest affinities is predominantly attributed to variations in hydrogen bond energy (Fig. 3). In fact, different conformations of a ligand can significantly alter the positioning of functional groups within the active site of a receptor, affecting the formation and strength of hydrogen bonds with amino acid residues. This conformational variability can lead to substantial differences in binding affinity. As illustrated in Fig. S5, which presents the binding “fingerprints” of 4-[(4-aminophenyl)methyl]aniline conformers with the highest and lowest affinities, these differences are evident. The conformer with the highest binding affinity is positioned within the active site in such a way that its functional groups are optimally aligned to form hydrogen bonds with several key amino acid residues, including Glu353, His524, Gly521, and Leu346. This favorable interaction pattern enhances its overall binding strength. In contrast, the conformer with the lowest binding affinity is only able to form hydrogen bonds with Glu353, limiting its interaction with the receptor and resulting in a weaker binding affinity. For ligands containing chlorine or methoxide groups, shape energy becomes the predominant factor influencing the ΔG difference between conformers. This indicates that, for these types of ligands, the geometric fit and shape complementarity between the ligand and the receptor are more critical in determining binding affinity than hydrogen bonding. The alignment and spatial arrangement of the ligand's functional groups in relation to the receptor's binding site are crucial for achieving optimal interaction. Consequently, the differences in binding affinity observed across various conformers of these ligands are closely tied to how well each conformer matches the receptor's shape and structural requirements.

Fig. 3 Differences in interaction energies between BPA analogs ligand conformers with the highest and lowest binding affinities against the apo form of the human estrogen receptor.

Comparison of the estrogenic activity of bisphenol analogs with well-known EDCs

The estrogenic activity of bisphenol derivatives were assessed in comparison to those of established the Endocrine Disrupting Chemicals.⁶³ Our evaluation focused on assessing the binding affinity of selected EDCs, such as bisphenol B, 4-tert-octylphenol, phenol, etc., among a total of 40 compounds, with three distinct conformers of the human estrogen receptor α. The benchmark EDCs included in this analysis have well-documented estrogenic activity. For example, bisphenol B (BPB) has been shown to activate both ERα and ERβ, disrupt steroidogenesis, and induce reproductive toxicity in vivo via both genomic and non-genomic pathways.^81–83 Similarly, 4-tert-octylphenol has been widely used as a model alkylphenol EDC. It binds ERα with high affinity and stimulates uterine growth in uterotrophic assays.⁸⁴ An integrated approach to testing and assessment (IATA) case study demonstrated that substituted phenols, including alkyl- and aryl-substituted analogues, can be reliably predicted to activate the estrogen receptor using combined in vitro reporter gene assays and structure–activity relationship models.⁸⁵ Among the investigated EDCs when tested against the apo, agonist, and antagonist conformers of the estrogen receptor we observed diverse binding affinities. Particularly, binding affinities ranged between −13.65 kcal mol⁻¹ to −9.44 kcal mol⁻¹ for the apo conformer, −14.11 kcal mol⁻¹ to −9.23 kcal mol⁻¹ for the agonist conformer, and −14.18 kcal mol⁻¹ to −9.87 kcal mol⁻¹ for the antagonist conformer. Contrarily, the binding affinity ranges for the bisphenol analogs exhibited almost same range values, spanning from −13.83 kcal mol⁻¹ to −9.57 kcal mol⁻¹ for the apo conformer, −13.72 kcal mol⁻¹ to −8.91 kcal mol⁻¹ for the agonist conformer, and −14.16 kcal mol⁻¹ to −9.95 kcal mol⁻¹ for the antagonist conformer (Fig. 4). To visualize binding trends between benchmark EDCs and bisphenol analogs, we mapped shape energy against H-bond energy (Fig. S6). As shown in Fig. S6, most compounds clustered within a shape energy range of −13 to −10 kcal mol⁻¹and an H-bond energy range of −3 to 0 kcal mol⁻¹. Some bisphenol analogs displayed stronger steric complementarity with the ERα binding pocket, as reflected by more negative shape energy values around −14 kcal mol⁻¹. These lower shape energy values indicate a better geometric fit between the ligand and the receptor's hydrophobic cavity, and suggests that these compounds may adopt conformations that maximize van der Waals interactions and minimize steric clashes. In parallel, other analogs demonstrated more favorable hydrogen bonding interactions, with H-bond energy values below −3 kcal mol⁻¹. This increase in polar interaction strength indicates the formation of stable hydrogen bonds with key polar amino acid residues in the ERα ligand-binding domain, such as Glu353, Arg394, or His524. The overlap in docking energies, particularly in the region of shape energy around −12 to −11 kcal mol⁻¹ and H-bond energy around −2 kcal mol⁻¹, suggests that some bisphenol analogs mimic the binding behavior of known EDCs, and analogs that fall outside the main EDC cluster may possess distinct binding profiles. These findings together suggest that the most bisphenol analogs analyzed in this study demonstrate characteristics commonly associated with endocrine disruptor chemicals, which can interfere with hormonal signaling pathways. Such interference can disrupt normal endocrine function and have adverse effects on physiological processes regulated by hormones.⁴⁴

Fig. 4 Binding affinity (ΔG (kcal mol⁻¹)) of known EDCs, bisphenols, E2, 4-OHT, and 4-tert-butylphenol against three conformations of human estrogen-related receptor ligand binding domain (LBD) (A) binding affinity of known EDCs, binding affinity of bisphenols against (B) apo conformation (PDB ID: 1A52), (C) agonist conformation (PDB ID: 1gwr), (D) Antagonist conformation (PDB ID: 3ert). E2 and 4-OHT are designated as positive controls, while 4-tert-butylphenol represents the lowest-affinity compound in the benchmark set.

Is octanol–water partitioning coefficient a good predictor for the binding affinities?

After examining the docking results, it became evident that BPA analogs showed an overall strong binding affinity toward three different forms of estrogen receptors (ER) when compared to their counterpart BPF and BPS analogs (Fig. 4B–D). One plausible explanation for the observed differences lies in the nature of the binding site on the ER. Among the amino acids within the ligand binding site, Glu353, Arg394, and His524 possess side chains capable of forming hydrogen bonds with the ligands. Apart from these three residues, the binding site mainly consists of hydrophobic amino acids, making it well-suited for accommodating the endogenous steroid ligand, estradiol (E2).⁷⁹ The calculated log K_OW values for BPA, BPF, and BPS were 3.25, 2.63, and 0.65, respectively (Table S6), and among these three compounds, BPA with higher log K_OW showed higher binding affinity (Fig. S7). The K_OW values indicate the level of hydrophobicity of the bisphenols, which could influence their interactions with ERα. Since the binding site on the ER tends to prefer hydrophobic interactions, it is anticipated that analogs with higher hydrophobicity, like BPA analogs, are more likely to bind strongly. However, this rule does not apply to each category of bisphenol analogues with different functional groups (Tables S3–S5). For example, in the BPA analog, 4-[2-(4-aminophenyl)propan-2-yl]aniline has the lowest log K_OW value (2.70) but high binding affinity (−13.64, −13.72, and −13.99 kcal mol⁻¹ against apo, agonist, and antagonist conformers of ERα, respectively) in the category (Fig. S8). This indicates that increased hydrophobicity alone does not significantly impact the binding affinities without considering the functional groups that can form hydrogen bonds and the steric effects, which appear to be additional factors in the reduction of estrogenic activity (EA).

Screening for agonist/antagonist differentials

Docking analysis of chemicals behavior towards ERα conformations indicated that the activation conformation offers greater affinity for agonist E2 (−16.59 kcal mol⁻¹ toward agonist conformation and −14.16 kcal mol⁻¹ toward antagonist conformation), whereas antagonist 4-OHT exhibited higher affinity in the inactivation conformation (−16.5 kcal mol⁻¹ toward antagonist conformation and −11.18 kcal mol⁻¹ toward agonist conformation). To classify the agonist or antagonist activity of chemicals, we used the agonist/antagonist factor as a parameter to evaluate conformational preferences based on binding affinity. Agonists display an affinity for binding to templates that already bound by agonists, whereas antagonists exhibit a preference for templates bound by antagonists. This differential binding preference significantly influences the biological activity of various ligands.⁷¹ By plotting the binding affinity (ΔG) and the agonist/antagonist factor (Cpf) for E2, 4-OHT, known EDCs, and bisphenol analogs, distinct patterns were observed (Fig. 5). E2 appeared as the most effective agonist, with the highest and the strongest affinity toward agonist activity. Conversely, 4-OHT was the most potent antagonist, displaying the strongest preference. When comparing the absolute values of the agonist/antagonist factor, it was observed that most BPA and BPF analogs fall within the range of −1 to 1. This range is used as a tolerance window that reflects the limitations of molecular docking accuracy. Docking scoring functions typically have an error margin of about ±1 kcal mol⁻¹, and many studies have noted that differences smaller than 1–2 kcal mol⁻¹ fall within the noise of the method and should not be over interpreted.^86–88 Energy differences within this range are generally not considered reliable for distinguishing between binding modes or affinities. This suggests that BPA and BPF derivatives have the potential to act as both agonists and antagonists of the ER based on their structural and environmental context. This ability emphasizes the dual nature of BPA and BPF derivatives, which can lead to varied biological effects and potential endocrine-disrupting properties. On the other hand, BPS analogs displayed a higher tendency against antagonism. This observation suggests that BPS derivatives are more likely to inhibit ER activity than promote it. This is supported by experimental studies showing BPA acts as an ER agonist in many human cell lines, promoting transcription and proliferation at micromolar levels, but can also display partial antagonist activity depending on concentration and context.^59,61 BPF similarly shows agonistic effects in breast cancer cells and zebrafish assays, though with less antagonist activity than BPA.⁵⁹ Additionally, bisphenol S (BPS) demonstrates both agonist and antagonist behaviors, with some evidence suggesting it acts more antagonistically toward ERα by reducing receptor protein levels in certain biological systems.^59,60 These findings validate the predicted binding preferences between ER conformations observed in computational studies.

Fig. 5 The analysis of conformational preferences method involved plotting the agonist/antagonist factor (ΔΔG) versus the calculated binding affinity values. ΔΔG was determined by subtracting the binding affinity of the agonist-bound receptor from that of the antagonist-bound receptor, based on the docking calculations. The small dots represent known endocrine-disrupting compounds (EDCs).

To strengthen the dual activity of BPA and BPF, a comparative analysis of docking fingerprints was performed in the agonist and antagonist conformations of ERα (Fig. S9 and S10). In both cases, the conformational preference (CPF) values were close to zero, indicating no strong preference for either conformation and supporting their dual binding capability. The results showed that BPA and BPF engage similar key residues in both conformations, displaying nearly identical interaction profiles. For BPA, hydrogen bonding occurs with Arg394, Glu353, and Thr347 in both conformations (Fig. S9), suggesting a stable binding mode that does not favor either the agonist or antagonist state, whereas BPF in the agonist conformation forms hydrogen bonds with Glu353, Gly521, and His524, but lacks the optimal hydrogen donor interaction with His524 in antagonist conformation (Fig. S10), which may account for its less antagonistic activity compared to BPA.

As shown in Fig. 5 bisphenol analogs functionalized with NH₂–Cl and NH₂–OCH₃ groups demonstrate enhanced binding affinity toward the antagonist conformation of the estrogen receptor (ER). The amino (NH₂) group serves as a strong electron-donating moiety, influencing the electron density distribution on the aromatic ring and promoting stronger polar interactions with receptor residues. Also, the presence of chlorine and methoxy substituents (steric bulk) increases the overall molecular volume, which may hinder the proper positioning of helix 12 in the AF-2 domain, and stabilize the receptor's antagonist conformation.⁸⁹ Bisphenol analogs with NH₂–NH₂, OH–OH, and NH₂–OH functional groups exhibit comparable binding affinities toward both the agonist and antagonist conformations of the estrogen receptor (ER). This dual affinity arises from the flexibility of these groups to act as both hydrogen bond donors and acceptors, which allows them to form favorable interactions within the distinct binding pocket of each receptor. Their relatively small size and polarity minimize steric hindrance and electronic bias, enabling accommodation without inducing significant conformational preference. Unlike bulkier or strong electron-donating/withdrawing substituents, these groups cannot selectively stabilize one receptor conformation over the other, resulting in similar binding energies for both agonist and antagonist forms.

The structural differences between BPA, BPF, and BPS are likely related to their varying degrees of antagonistic activity. Understanding these differences is crucial for assessing the safety and biological impact of these chemicals, especially given their widespread use in consumer products and potential health implications.

Protonation and deprotonation states of ligands on binding affinity

Protonation and deprotonation alter the electrostatic charges on the ligand can affect ionic interactions and hydrogen bonding with the receptor.⁹⁰ Protonated/deprotonated forms of the chemicals showed different binding affinity compared to their neutral form. In most cases, both the deprotonated and protonated forms of the ligands showed a decrease in binding affinity (Fig. S11, Tables S7–S9). This decrease can be attributed to hydrogen binding and electrostatic interactions as well as ligand desolvation. A significant decrease in the binding affinity among BPA analogs is related to BPA when it undergoes deprotonation, resulting in a −2 charge. In its neutral form, BPA exhibited binding energies of −13.35 kcal mol⁻¹, −13.46 kcal mol⁻¹, and −13.37 kcal mol⁻¹ against the apo, agonist, and antagonist forms of the estrogen receptor (ER), respectively. However, when ionized, the binding energies shifted to −10.26 kcal mol⁻¹, −9.78 kcal mol⁻¹, and −11.36 kcal mol⁻¹, respectively. This reduction is mainly due to the loss of hydrogen bonding. In its neutral form, BPA forms hydrogen bonds with Glu353, His524, and Gly521, with a hydrogen bonding energy of −2.86 kcal mol⁻¹. Upon deprotonation, however, BPA loses the capacity to form these hydrogen bonds, leading to a drastic reduction in hydrogen bonding energy to −0.28 kcal mol⁻¹. Additionally, the ligand desolvation energy in ionized form is higher compared to its neutral form (2.35 and 1.19 kcal mol⁻¹). This means the increased negative charge enhances solvation energy and makes the ligand more likely to remain in the solvent than binding to the protein, further diminishing the binding affinity. Also, in residue fingerprint showing which residues (amino acids) in the receptor site the ligand is interacting with (Fig. S12), the number of interacting amino acids decreased in overall binding strength, or it may bind differently, leading to a change in the biological activity of the ligand. These observations are consistent almost across BPF and BPS ionized ligands. The significant difference in binding affinity between the neutral/ionized BPF and BPS analogs is related to 4-[(4-aminophenyl)methyl]aniline and 4-(4-aminophenyl)sulfonyl aniline (with differences of 3.1 and 2.86 kcal mol⁻¹, respectively) when they undergo protonation, resulting in a +2 charge. This difference is mainly due to the ability to form hydrogen bonds, as the hydrogen bond energy decreases from −2.91 to −0.93 kcal mol⁻¹ in ionized 4-[(4-aminophenyl)methyl]aniline and from −2.60 to −0.07 kcal mol⁻¹ in ionized 4-(4-aminophenyl)sulfonyl aniline. This indicates that ionization can significantly alter binding affinity. This alteration, in turn, can influence the chemical's endocrine-disruptive properties. Ionization may affect how a ligand interacts with its target receptor, potentially enhancing or diminishing its ability to disrupt endocrine function.

Structure–activity analysis of bisphenol analogs

The common ERα agonist estradiol (E2) and the antagonist (4-OHT) were utilized as positive controls. The binding of E2 with binding affinity −16.59 kcal mol⁻¹ (against ERα agonist conformer) involves the formation of hydrogen bonds with Arg394, His524, Glu353, and Gly521 within the binding pocket. In contrast, 4-OHT with binding affinity −16.50 kcal mol⁻¹ (against ERα antagonist conformer) establishes hydrogen bonds with Glu353 in the binding site with the alkyl amine side chain of 4-OHT fitting into the antagonist structure (Fig. S13).

Building on these findings, BPA analogs, characterized by the presence of hydroxyl and amine groups on the rings, including compounds like 4-[2-(4-aminophenyl)propan-2-yl]phenol, 4-[2-(4-aminophenyl)propan-2-yl]aniline, and bisphenol A exhibited remarkable binding affinities across apo and agonist conformers of ERα by the formation of hydrogen bonds with Gly521, Arg394, His524, Leu346, and Glu353 in the binding site. They exhibited binding affinities to the apo conformer at −13.83 kcal mol⁻¹, −13.64 kcal mol⁻¹, and −13.35 kcal mol⁻¹, and to agonist conformers at −13.05 kcal mol⁻¹, −13.72 kcal mol⁻¹, and −13.46 kcal mol⁻¹, respectively. Notably, for antagonist conformer of ERα, compounds 4-[2-(4-aminophenyl)propan-2-yl]phenol, 4-[2-(4-aminophenyl)propan-2-yl]aniline, and 4-[2-(4-chlorophenyl)propan-2-yl]aniline showed highest affinity with binding affinities of −14.16 kcal mol⁻¹, −13.99 kcal mol⁻¹, and −13.72 kcal mol⁻¹, respectively (Fig. 6).

Fig. 6 3D visualization of bisphenols interactions within the active site of three conformers of ERα, highlighting residues involved in binding. Hydrogen bonds between bisphenols and residues are indicated by green dashed lines. (A–C) Display the binding pose of 4-[2-(4-aminophenyl)propan-2-yl]phenol, 4-[(4-aminophenyl) methyl]aniline, and 4-(4-aminophenyl)sulfonyl phenol within the active site of apo conformer. (D–F) Show the binding pose of 4-[2-(4-aminophenyl)propan-2-yl]aniline, 4-[(4-aminophenyl)methyl]phenol, and 4-(4-aminophenyl)sulfonyl phenol within the active site of agonist conformer. (G–I) Present the binding pose of 4-[2-(4-aminophenyl)propan-2-yl]phenol, 4-[(4-chlorophenyl) methyl]aniline, and 4-(4-aminophenyl)sulfonyl aniline within the active site of agonist conformer.

Interestingly, among the BPF analogs, compounds such as 4-[(4-aminophenyl)methyl]aniline, 4-[(4-aminophenyl)methyl]phenol, and bisphenol F exhibited significant binding affinities towards both the apo and agonist conformers of ERα. This affinity is primarily due to the formation of hydrogen bonds with key residues His524, Gly521, Leu346, and Glu353 within the binding site. For the apo conformer, the binding affinities were recorded at −13.41 kcal mol⁻¹, −13.28 kcal mol⁻¹, and −13.07 kcal mol⁻¹, respectively. Similarly, for the agonist conformer, these compounds showed binding affinities of −12.97 kcal mol⁻¹, −13.32 kcal mol⁻¹, and −12.61 kcal mol⁻¹, respectively. In addition to their interaction with the apo and agonist conformers, these compounds also demonstrated notable binding affinities (establishing hydrogen bonds with Glu353 and Leu346 in the binding site) towards the antagonist conformer of ERα. Specifically, 4-[(4-chlorophenyl) methyl]aniline, 4-[(4-aminophenyl)methyl]aniline, and 4-[(4-methoxyphenyl)methyl]aniline exhibited binding affinities of −12.94 kcal mol⁻¹, −12.77 kcal mol⁻¹, and −12.76 kcal mol⁻¹, respectively (Fig. 6).

For BPS analogs, 4-(4-aminophenyl)sulfonyl phenol, bisphenol S, and 4-(4-chlorophenyl)sulfonyl phenol showed highest binding affinity towards apo and agonist conformers. They exhibited binding affinities to the apo conformer at −11.06 kcal mol⁻¹, −10.65 kcal mol⁻¹, and −10.49 kcal mol⁻¹, and to agonist conformers at −10.78 kcal mol⁻¹, −10.52 kcal mol⁻¹, and −10.27 kcal mol⁻¹, respectively. While compounds 4-(4-aminophenyl)sulfonyl aniline, bisphenol S, and 4-(4-aminophenyl)sulfonyl phenol with binding affinity of −12.69 kcal mol⁻¹, −12.29 kcal mol⁻¹, and −12.24 kcal mol⁻¹, respectively, showed more affinity to antagonist conformer (Fig. 6). Fig. S14 shows the superimposition of BPA, BPF, and BPS analogs with the highest binding affinity toward the three conformers of ER, comparing how each ligand fits into the binding site.

These analysis revealed that bisphenol analogs with high agonistic activity exhibited interactions with many of the same key amino acid residues observed for the natural agonist 17β-estradiol (E2), including Glu353, Arg394, His524, Leu387, Leu391, and Gly521, and compounds with higher affinity for the antagonist conformation display interaction fingerprints resembling those of 4-hydroxytamoxifen (4-OHT), with amino acid residues within or near the AF-2 region of ERα (Fig. S13 and S14). This similarity suggests that these compounds may be associated with H12 positioning. Since H12 acts as a conformational switch controlling coactivator recruitment, ligand interactions with these residues likely modulate the receptor's functional state.^27,79 Like E2 and 4-OHT, which are known to stabilize agonist and antagonist conformations by altering H12 orientation respectively, our compounds may similarly influence the dynamic positioning of H12, thereby affecting agonist versus antagonist activity. From the data, it is evident that bisphenol analogs with hydroxyl (OH) and amine (NH₂) groups at the para-position exhibit higher binding affinities compared to those with methoxide (OCH₃) and Cl groups at the same position. This difference in binding affinity can be attributed to the hydrogen bonding capabilities of the hydroxyl and amine groups. These functional groups can form strong hydrogen bonds with key amino acid residues within the binding site of the target receptor, such as Glu353 and Arg394. Hydroxyl groups are known to be excellent hydrogen bond donors due to the presence of a polar oxygen–hydrogen bond, which can interact favorably with electron-rich acceptor atoms in the binding pocket. Similarly, amine groups can act as both hydrogen bond donors and acceptors, further enhancing their ability to form stable interactions with the receptor.

The number of hydrogen bond acceptors and donors in a ligand can correlate with its binding affinity to a target protein. Generally, ligands with higher counts of both hydrogen bond acceptors and donors are expected to exhibit higher binding affinity. 4-[2-(4-Aminophenyl)propan-2-yl]phenol that has 2 hydrogen bond acceptors and 2 donors, showed the highest binding affinity among the tested compounds. The lowest binding affinity, however, was not observed in compounds with chlorine at the R1 and R2 positions, which have 0 acceptors and 0 donors. Instead, it was found in compounds featuring methoxide groups at the R1 and R2 positions, which have 2 acceptors and 0 donors. Methoxide groups, while capable of participating in hydrogen bonding to a lesser extent, primarily act as hydrogen bond acceptors and lack the versatility of hydroxyl and amine groups. Chlorine atoms, being highly electronegative, do not participate in hydrogen bonding and instead rely on van der Waals interactions, which are generally weaker than hydrogen bonds. Consequently, derivatives with methoxide and chlorine groups at the para-position demonstrate lower binding affinities due to their reduced capacity to form strong and stable interactions with the receptor. These observations highlight the critical role of functional group chemistry in determining the binding affinity of bisphenol derivatives and underscore the importance of hydrogen bonding in receptor–ligand interactions.

Concluding remarks and outlook on the uses of molecular docking for contaminant assessment

This study provides a comprehensive ranking of bisphenol analogs and well-known endocrine-disrupting chemicals (EDCs) based on docking scores. It facilitates a better understanding of their agonist and antagonist potential relative to known reference compounds. These findings show how minor structural differences among investigated compounds can affect their interaction with the receptor's active site. Among the investigated compounds, several bisphenol analogs demonstrated binding affinity almost similar to traditional EDCs. Within this benchmark, 4-[2-(4-aminophenyl)propan-2-yl]phenol exhibited the strongest binding affinity comparable to that of reference EDCs such as bisphenol B against three conformers of ER. Conversely, 1-methoxy-4-(4-methoxyphenyl)sulfonyl benzene showed the weakest binding affinity, indicating lower estrogenic potential relative to the benchmark set. Our analysis further revealed that adding certain functional groups, such as amines and hydroxyl groups, can enhance EDC activity by increasing binding affinity toward the estrogen receptor through hydrogen bonds. These stronger interactions often lead to increased agonistic activity, which means these compounds can activate the receptor more effectively and disrupt normal hormonal functions. In contrast, introducing groups like methoxy (–OCH₃) or bulky substituents can change the compounds' interaction within the receptor, sometimes promoting antagonistic behavior that blocks receptor activation. The presence of amine groups in these antagonists can also help form alternative hydrogen bonds that stabilize the receptor in its inactive state.

From this study, we demonstrated that it is essential to consider factors that affect binding affinity to ensure proper setting up of the parameters when applying molecular docking. Factors to consider when applying molecular docking include solvent water molecules within the binding pocket within the receptor, ligand conformers, ligands' ionization states, and structural features of ligand and receptors. Water molecules within the receptor's ligand-binding pocket can play important roles due to hydrogen bonds between ligands and residues.⁷⁷ In cases when hydrogen bond has lower contributions to the binding than the hydrophobic interactions, whether to include water molecules has little effect on the predicted binding affinity. Therefore, to decide whether water molecules within the binding pocket should be included, it is essential to understand the interactions between the ligands and binding sites at the receptor. The ionization state of a ligand also has an impact on its binding affinity at the receptor. Protonation and deprotonation states can change charge distributions and alter electrostatic interactions as well as hydrogen bonds within the binding site.⁹⁰ Similar to the treatment of different conformers of a ligand by adopting the maximum binding affinity among those derived for all conformers as the binding score, the common practice for virtual screening using molecular docking only considers the neutral form of ionizable chemicals,^21,24,25 because the neutral species would result in higher binding affinity due to lower desolvation energy from the water. Theoretically, we can derive the binding affinity of ioniazable chemicals by considering the speciation and binding affinities of natural and ionizable forms with their pK_a and pH. However, the pH at the microenvironment within the binding pocket of the receptor is often unknown and can be highly dynamic and influenced by numerous factors. For screening purpose, we argue that it is unnecessary and for most cases impossible to derive a bulk binding affinity considering the speciation of ionization chemicals. However, instead of deriving the binding affinities for neutral species only, deriving the binding affinity and pK_a for ionization species would help determine chemicals whose bulk binding affinity can hardly vary due to speciation at pH range of organisms. Additional analysis from this study also suggests that octanol–water partition coefficient, a useful value for predicting hydrophobicity, is not always a good predictor for binding affinity, as other factors also play important roles in ligand binding affinity.

As a high throughput in silico tool, molecular docking is invaluable for screening anthropogenic chemicals on their potential toxicity determined by well-known mechanism of actions and toxicity pathways.^20,22 This approach would be useful to assess newly synthesized chemicals as replacements of existing commercial chemicals or environmental transformation products following similar backbone structures. Despite of alternative approaches such as QSARs and machine learning models to characterize interactions between ligands and receptors, molecular docking quantifies the interactions of atom pairs between ligands and their target proteins as a function of distance.¹⁵ Such mechanistical simulations of intermolecular interactions make molecular docking not restricted to a given protein or group of ligands used to develop QSARs or machine learning models. As such, molecular docking is applicable to a larger chemical space. Additionally, molecular docking can tackle different proteins relevant to different biochemical processes or mechanisms of actions for toxic effects. From this study, we demonstrate the importance of the 3D shape of molecules and their conformations would affect their binding to proteins such interactions can be described by QSARs based on 2D molecular descriptors.

Despite capturing 3D shape related energy, it is worth noting that the force field used for molecular docking is a simplification. Molecular docking commonly assumes rigid structures for both ligands and receptors and the dynamic processes of the chemical system are not considered. This treatment makes molecular docking unable to capture the full complexity of molecular interactions in vivo. However, the dynamic nature of the receptor can be partly accounted for in molecular docking by considering receptors of different conformations as snapshots of the highly dynamic structures of the protein.^15,70 Due to the simplifications of molecular docking, the binding energy output should not be used directly and docking results from different docking software with different force fields and treatment of solvation processes should not be compared. Instead, the binding energy from molecular docking should be viewed as a scoring system for ranking and screening chemicals. To apply docking to screen new sets of chemicals, it is essential to include control and benchmark chemicals that are well known for their interactions with the protein and mechanism of action to intervene biochemical processes causing effects such as endocrine disruption.⁹¹

Overall, molecular docking has emerged as a valuable approach for screening chemicals in environmental and exposure assessments.^18,92 Docking allows rapid screening of large chemical libraries and prioritizes compounds for further experimental validation based on predicted binding affinities and interactions. Combining docking results with experimental data (e.g., toxicity assays, bioavailability studies) enhances confidence in predictions and improves decision-making.^92–94 There is great potential for molecular docking to be used as part of NAMs for assessment and screening of new and existing industrial chemicals.⁹⁵

Conflicts of interest

There are no conflicts to declare.

Data availability

Additional text, figures and tables supporting this article have been included as part of the SI. Supplementary information is available. See DOI: https://doi.org/10.1039/d5em00084j.

Acknowledgements

We acknowledge funding support from Natural Science and Engineering Research Council (NSERC) of Canada (Discovery Grant RGPIN-2021-04132) and Environment and Climate Change Canada (GCXE22S068) to X. Z. We thank OpenEye Scientific for the Academic Licence of OEDocking.

References

1. R. Naidu, B. Biswas, I. R. Willett, J. Cribb, S. B. Kumar and N. C. Paul, et al., Chemical pollution: A growing peril and potential catastrophic risk to humanity, Environ. Int., 2021, 156, 106616.
2. D. Nikolopoulou, E. Ntzani, K. Kyriakopoulou, C. Anagnostopoulos and K. Machera, Priorities and Challenges in Methodology for Human Health Risk Assessment from Combined Exposure to Multiple Chemicals, Toxics, 2023, 11(5), 401.
3. T. Hartung, Toxicology for the twenty-first century, Nature, 2009, 460(7252), 208–212.
4. F. A. P. C. Gobas, P. Mayer, T. F. Parkerton, R. M. Burgess, D. v. d. Meent and T. Gouin, A chemical activity approach to exposure and risk assessment of chemicals, Environ. Toxicol. Chem., 2018, 37(5), 1235–1251.
5. S. Auerbach, W. Casey, D. Chang, J. Cowden, S. Davidson-Fritz, and J. Dean, et al., Scientific Studies Supporting Development of Transcriptomic Points of Departure for EPA Transcriptomic Assessment Products (ETAPs), EPA/600/X-23/084, U.S. Environmental Protection Agency, Washington, DC, 2024,  DOI:10.23645/epacomptox.25365550.
6. R. R. Diniz, T. F. S. Domingos, G. R. Pinto, L. M. Cabral, M. d. Pádula and A. M. T. d. Souza, Use of in silico and in vitro methods as a potential new approach methodologies (NAMs) for (photo)mutagenicity and phototoxicity risk assessment of agrochemicals, Sci. Total Environ., 2023, 904, 167320.
7. B. J. Neves, R. C. Braga, C. C. Melo-Filho, J. T. Moreira-Filho, E. N. Muratov and C. H. Andrade, QSAR-Based Virtual Screening: Advances and Applications in Drug Discovery, Front. Pharmacol, 2018, 9, 1275.
8. A. D. Guardo, T. Gouin, M. MacLeod and M. Scheringer, Environmental fate and exposure models: advances and challenges in 21 st century chemical risk assessment, Environ. Sci.: Processes Impacts, 2018, 20(1), 58–71.
9. X. Sun, X. Zhang, D. C. G. Muir and E. Y. Zeng, Identification of Potential PBT/POP-Like Chemicals by a Deep Learning Approach Based on 2D Structural Features, Environ. Sci. Technol., 2020, 54(13), 8221–8231.
10. M. T. D. Cronin, J. C. Madden, C. Yang and A. P. Worth, Unlocking the potential of in silico chemical safety assessment – A report on a cross-sector symposium on current opportunities and future challenges, Comput. Toxicol., 2019, 10, 38–43.
11. X. Zhang, X. Sun, R. Jiang, E. Y. Zeng, E. M. Sunderland and D. C. G. Muir, Screening New Persistent and Bioaccumulative Organics in China's Inventory of Industrial Chemicals, Environ. Sci. Technol., 2020, 54(12), 7398–7408.
12. M. P. Gleeson, S. Modi, A. Bender, L. M. R. Robinson, J. Kirchmair and M. Promkatkaew, et al., The Challenges Involved in Modeling Toxicity Data In Silico: A Review, Curr. Pharm. Des., 2012, 18(9), 1266–1291.
13. J. Kirchmair, M. J. Williamson, J. D. Tyzack, L. Tan, P. J. Bond and A. Bender, et al., Computational Prediction of Metabolism: Sites, Products, SAR, P450 Enzyme Dynamics, and Mechanisms, J. Chem. Inf. Model., 2012, 52(3), 617–648.
14. F. Gentile, J. C. Yaacoub, J. Gleave, M. Fernandez, A.-T. Ton and F. Ban, et al., Artificial intelligence–enabled virtual screening of ultra-large chemical libraries with deep docking, Nat. Protoc., 2022, 17(3), 672–697.
15. J. M. Paggi, A. Pandit and R. O. Dror, The Art and Science of Molecular Docking, Annu. Rev. Biochem., 2024, 93(1), 389–410.
16. M. T. Muhammed and E. Aki-Yalcin, Molecular Docking: Principles, Advances, and Its Applications in Drug Discovery, Lett. Drug Des. Discovery, 2024, 21(3), 480–495.
17. C. A. Ng and K. Hungerbuehler, Exploring the Use of Molecular Docking to Identify Bioaccumulative Perfluorinated Alkyl Acids (PFAAs), Environ. Sci. Technol., 2015, 49(20), 12306–12314.
18. S. D. Walker and S. McEldowney, Molecular docking: a potential tool to aid ecotoxicity testing in environmental risk assessment of pharmaceuticals, Chemosphere, 2013, 93(10), 2568–2577.
19. Y. Hong, Y. Wang, D. Wang, Q. Yuan, Z. Yang and C. Deng, Assessing male reproductive toxicity of environmental pollutant di-ethylhexyl phthalate with network toxicology and molecular docking strategy, Reprod. Toxicol., 2024, 130, 108749.
20. Z. Liu, Y. Liu, G. Zeng, B. Shao, M. Chen and Z. Li, et al., Application of molecular docking for the degradation of organic pollutants in the environmental remediation: A review, Chemosphere, 2018, 203, 139–150.
21. L. Celik, J. D. D. Lund and B. Schiøtt, Exploring Interactions of Endocrine-Disrupting Compounds with Different Conformations of the Human Estrogen Receptor α Ligand Binding Domain: A Molecular Docking Study, Chem. Res. Toxicol., 2008, 21(11), 2195–2206.
22. H. Tan, Q. Chen, H. Hong, E. Benfenati, G. C. Gini and X. Zhang, et al., Structures of Endocrine-Disrupting Chemicals Correlate with the Activation of 12 Classic Nuclear Receptors, Environ. Sci. Technol., 2021, 55(24), 16552–16562.
23. D. Trisciuzzi, D. Alberga, K. Mansouri, R. Judson, S. Cellamare and M. Catto, et al., Docking-based classification models for exploratory toxicology studies on high-quality estrogenic experimental data, Future Med. Chem., 2015, 7(14), 1921–1936.
24. T. Nose, T. Tokunaga and Y. Shimohigashi, Exploration of endocrine-disrupting chemicals on estrogen receptor alpha by the agonist/antagonist differential-docking screening (AADS) method: 4-(1-adamantyl)phenol as a potent endocrine disruptor candidate, Toxicol. Lett., 2009, 191(1), 33–39.
25. H. Kesamaru, K. Suyama and T. Nose, Importance of Receptor Conformations in Docking Calculation-Based Risk Assessment for Endocrine Disruptors against Estrogen Receptor α, ACS Omega, 2019, 4(4), 6620–6629.
26. A. Vuorinen, A. Odermatt and D. Schuster, In silico methods in the discovery of endocrine disrupting chemicals, J. Steroid Biochem. Mol. Biol., 2013, 137, 18–26.
27. P. Germain, L. Altucci, W. Bourguet, C. Rochette-Egly and H. Gronemeyer, Nuclear receptor superfamily: Principles of signaling, Pure Appl. Chem., 2003, 75(11–12), 1619–1664.
28. H. Gronemeyer, J.-A. Gustafsson and V. Laudet, Principles for modulation of the nuclear receptor superfamily, Nat. Rev. Drug Discovery, 2004, 3(11), 950–964.
29. J. T. Moore, J. L. Collins and K. H. Pearce, The nuclear receptor superfamily and drug discovery, ChemMedChem, 2006, 1(5), 504–523.
30. A. Berkenstam and J.-A. Gustafsson, Nuclear receptors and their relevance to diseases related to lipid metabolism, Curr. Opin. Pharmacol., 2005, 5(2), 171–176.
31. W. Bourguet, P. Germain and H. Gronemeyer, Nuclear receptor ligand-binding domains: three-dimensional structures, molecular interactions and pharmacological implications, Trends Pharmacol. Sci., 2000, 21(10), 381–388.
32. R. Kumar and E. B. Thompson, Transactivation functions of the N-terminal domains of nuclear hormone receptors: protein folding and coactivator interactions, Mol. Endocrinol., 2003, 17(1), 1–10.
33. D. E. Frigo, M. Bondesson and C. Williams, Nuclear receptors: from molecular mechanisms to therapeutics, Essays Biochem., 2021, 65(6), 847–856.
34. M. Pawlak, P. Lefebvre and B. Staels, General Molecular Biology and Architecture of Nuclear Receptors, Curr. Top. Med. Chem., 2012, 12(6), 486–504.
35. J. M. Wurtz, W. Bourguet, J. P. Renaud, V. Vivat, P. Chambon and D. Moras, et al., A canonical structure for the ligand-binding domain of nuclear receptors, Nat. Struct. Biol., 1996, 3(1), 87–94.
36. D. M. Tanenbaum, Y. Wang, S. P. Williams and P. B. Sigler, Crystallographic comparison of the estrogen and progesterone receptor's ligand binding domains, Proc. Natl. Acad. Sci. U. S. A., 1998, 95(11), 5998.
37. A. Wärnmark, E. Treuter, J.-Å. Gustafsson, R. E. Hubbard, A. M. Brzozowski and A. C. W. Pike, Interaction of Transcriptional Intermediary Factor 2 Nuclear Receptor Box Peptides with the Coactivator Binding Site of Estrogen Receptor α, J. Biol. Chem., 2002, 277(24), 21862–21868.
38. A. K. Shiau, D. Barstad, P. M. Loria, L. Cheng, P. J. Kushner and D. A. Agard, et al., The structural basis of estrogen receptor/coactivator recognition and the antagonism of this interaction by tamoxifen, Cell, 1998, 95(7), 927–937.
39. A. Shankar and S. Teppala, Relationship between urinary bisphenol A levels and diabetes mellitus, J. Clin. Endocrinol. Metab., 2011, 96(12), 3822–3826.
40. F. S. Vom Saal, S. C. Nagel, B. L. Coe, B. M. Angle and J. A. Taylor, The estrogenic endocrine disrupting chemical bisphenol A (BPA) and obesity, Mol. Cell. Endocrinol., 2012, 354(1–2), 74–84.
41. D. Melzer, N. J. Osborne, W. E. Henley, R. Cipelli, A. Young and C. Money, et al., Urinary bisphenol A concentration and risk of future coronary artery disease in apparently healthy men and women, Circulation, 2012, 125(12), 1482–1490.
42. E. Jašarević, P. T. Sieli, E. E. Twellman, T. H. Welsh, T. R. Schachtman and R. M. Roberts, et al., Disruption of adult expression of sexually selected traits by developmental exposure to bisphenol A, Proc. Natl. Acad. Sci. U. S. A., 2011, 108(28), 11715–11720.
43. S. Jenkins, J. Wang, I. Eltoum, R. Desmond and C. A. Lamartiniere, Chronic oral exposure to bisphenol A results in a nonmonotonic dose response in mammary carcinogenesis and metastasis in MMTV-erbB2 mice, Environ. Health Perspect., 2011, 119(11), 1604–1609.
44. J. Dueñas-Moreno, A. Mora, M. Kumar, X.-Z. Meng and J. Mahlknecht, Worldwide risk assessment of phthalates and bisphenol A in humans: The need for updating guidelines, Environ. Int., 2023, 181, 108294.
45. A. Beronius, C. Rudén, H. Håkansson and A. Hanberg, Risk to all or none? A comparative analysis of controversies in the health risk assessment of Bisphenol A, Reprod. Toxicol., 2010, 29(2), 132–146.
46. V. Mustieles, J.-P. Arrebola and M. Porta, From old pollutants to the regulation of bisphenol A: Lessons learned for health promotion and disease prevention, Prev. Med., 2023, 169, 107460.
47. L. N. Vandenberg, Exposure to bisphenol A in Canada: invoking the precautionary principle, Can. Med. Assoc. J., 2011, 183(11), 1265–1270.
48. G. D. Bittner, M. S. Denison, C. Z. Yang, M. A. Stoner and G. He, Chemicals having estrogenic activity can be released from some bisphenol A-free, hard and clear, thermoplastic resins, Environ. Health: Global Access Sci. Source, 2014, 13, 103.
49. S. Eladak, T. Grisin, D. Moison, M.-J. Guerquin, T. N'Tumba-Byn and S. Pozzi-Gaudin, et al., A new chapter in the bisphenol A story: bisphenol S and bisphenol F are not safe alternatives to this compound, Fertil. Steril., 2015, 103(1), 11–21.
50. Y. Niu, B. Wang, R. Yang, Y. Wu, Y. Zhao and C. Li, et al., Bisphenol Analogues and Their Chlorinated Derivatives in Breast Milk in China: Occurrence and Exposure Assessment, J. Agric. Food Chem., 2021, 69(4), 1391–1397.
51. R. Viñas, R. M. Goldblum and C. S. Watson, Rapid estrogenic signaling activities of the modified (chlorinated, sulfonated, and glucuronidated) endocrine disruptor bisphenol A, Endocr. Disruptors, 2013, 1(1), e25411.
52. N. Plattard, A. Dupuis, V. Migeot, S. Haddad and N. Venisse, An overview of the literature on emerging pollutants: Chlorinated derivatives of Bisphenol A (ClxBPA), Environ. Int., 2021, 153, 106547.
53. L. Gustavson, T. M. Ciesielski, J. Bytingsvik, B. Styrishave, M. Hansen and E. Lie, et al., Hydroxylated polychlorinated biphenyls decrease circulating steroids in female polar bears (Ursus maritimus), Environ. Res., 2015, 138, 191–201.
54. K. Hayakawa, Y. Onoda, C. Tachikawa, S. Hosoi, M. Yoshita and S. W. Chung, et al., Estrogenic/Antiestrogenic Activities of Polycyclic Aromatic Hydrocarbons and Their Monohydroxylated Derivatives by Yeast Two-Hybrid Assay, J. Health Sci., 2007, 53(5), 562–570.
55. X. Ning, N. Du, X. Zhang, S. Wang, Y. Zhi and Z. Li, et al., Metastatic effects of hydroxy-polycyclic aromatic hydrocarbons on liver cancer cells mediated by estrogen receptor α, Sci. Total Environ., 2024, 952, 175878.
56. S. Takayanagi, T. Tokunaga, X. Liu, H. Okada, A. Matsushima and Y. Shimohigashi, Endocrine disruptor bisphenol A strongly binds to human estrogen-related receptor γ (ERRγ) with high constitutive activity, Toxicol. Lett., 2006, 167(2), 95–105.
57. L. Barraza, A New Approach for Regulating Bisphenol A for the Protection of the Public's Health, J. Law Med. Ethics, 2013, 41(S1), 9–12.
58. M. Cao, J. Wei, Y. Pan, L. Wang, Z. Li and Y. Hu, et al., Antagonistic mechanisms of bisphenol analogues on the estrogen receptor α in zebrafish embryos: Experimental and computational studies, Sci. Total Environ., 2023, 857(Pt 1), 159259.
59. C. Park, H. Song, J. Choi, S. Sim, H. Kojima and J. Park, et al., The mixture effects of bisphenol derivatives on estrogen receptor and androgen receptor, Environ. Pollut., 2020, 260, 114036.
60. H. Lee, J. Park and K. Park, Mixture Effects of Bisphenol A and Its Structural Analogs on Estrogen Receptor Transcriptional Activation, Toxics, 2023, 11(12), 986.
61. Y. Li, K. A. Burns, Y. Arao, C. J. Luh and K. S. Korach, Differential estrogenic actions of endocrine-disrupting chemicals bisphenol A, bisphenol AF, and zearalenone through estrogen receptor α and β in vitro, Environ. Health Perspect., 2012, 120(7), 1029–1035.
62. J.-M. Molina-Molina, E. Amaya, M. Grimaldi, J.-M. Sáenz, M. Real and M. F. Fernández, et al., In vitro study on the agonistic and antagonistic activities of bisphenol-S and other bisphenol-A congeners and derivatives via nuclear receptors, Toxicol. Appl. Pharmacol., 2013, 272(1), 127–136.
63. Substances identified as endocrine disruptors at EU level | Endocrine Disruptor List.
64. Free Chemical Drawing Software for Students | ChemSketch.
65. SZYBKI 2.7.0.3, OpenEye, Cadence Molecular Sciences, Santa Fe, NM, https://www.eyesopen.com.
66. OMEGA 5.0.0.3, OpenEye, Cadence Molecular Sciences, Santa Fe, NM, https://www.eyesopen.com.
67. P. C. D. Hawkins, A. G. Skillman, G. L. Warren, B. A. Ellingson and M. T. Stahl, Conformer generation with OMEGA: algorithm and validation using high quality structures from the Protein Databank and Cambridge Structural Database, J. Chem. Inf. Model., 2010, 50(4), 572–584.
68. T. Long, J. K. Yao, J. Li, Z. Z. Kirshner, D. Nelson and G. G. Dougherty, et al., Estradiol and selective estrogen receptor agonists differentially affect brain monoamines and amino acids levels in transitional and surgical menopausal rat models, Mol. Cell. Endocrinol., 2019, 496, 110533.
69. J. S. Lewis and V. C. Jordan, Selective estrogen receptor modulators (SERMs): Mechanisms of anticarcinogenesis and drug resistance, Mutat. Res. Fundam. Mol. Mech. Mutagenesis, 2005, 591(1), 247–263.
70. A. M. Khan, J. Shawon and M. A. Halim, Multiple receptor conformers based molecular docking study of fluorine enhanced ethionamide with mycobacterium enoyl ACP reductase (InhA), J. Mol. Graphics Modell., 2017, 77, 386–398.
71. D. P. McDonnell, The molecular determinants of estrogen receptor pharmacology, Maturitas, 2004, 48, 7–12.
72. Spruce 1.5.3.3, OpenEye, Cadence Molecular Sciences, Santa Fe, NM, https://www.eyesopen.com.
73. M. McGann, FRED pose prediction and virtual screening accuracy, J. Chem. Inf. Model., 2011, 51(3), 578–596.
74. M. McGann, FRED and HYBRID docking performance on standardized datasets, J. Comput.-Aided Mol. Des., 2012, 26(8), 897–906.
75. OEDOCKING 4.3.1.0, OpenEye, Cadence Molecular Sciences, Santa Fe, NM, https://www.eyesopen.com.
76. VIDA 5.0.5.3, OpenEye, Cadence Molecular Sciences, Santa Fe, NM, https://www.eyesopen.com.
77. B. C. Roberts and R. L. Mancera, Ligand-protein docking with water molecules, J. Chem. Inf. Model., 2008, 48(2), 397–408.
78. S. Mukherjee, T. E. Balius and R. C. Rizzo, Docking Validation Resources: Protein Family and Ligand Flexibility Experiments, J. Chem. Inf. Model., 2010, 50(11), 1986–2000.
79. L. Celik, J. D. D. Lund and B. Schiøtt, Conformational Dynamics of the Estrogen Receptor α: Molecular Dynamics Simulations of the Influence of Binding Site Structure on Protein Dynamics, Biochemistry, 2007, 46(7), 1743–1758.
80. J. A. Grant and B. T. Pickup, A Gaussian Description of Molecular Shape, J. Phys. Chem., 1995, 99(11), 3503–3510.
81. N. Chen, X. Li, S. Zhou, X. Peng, S. Xue and Y. Liu, et al., Prenatal Bisphenol B Exposure Induces Adult Male Offspring Reproductive Dysfunction via ERα Inhibition-Triggered MHC I-Mediated Testicular Immunological Responses, Toxics, 2025, 13(6), 423.
82. H. Serra, C. Beausoleil, R. Habert, C. Minier, N. Picard-Hagen and C. Michel, Evidence for Bisphenol B Endocrine Properties: Scientific and Regulatory Perspectives, Environ. Health Perspect., 2019, 127(10), 106001.
83. H.-M. Li, Y.-Y. Li, Y.-C. Zhang, J.-B. Li, H.-M. Xu and Y.-M. Xiong, et al., Bisphenol B disrupts testis differentiation partly via the estrogen receptor-mediated pathway and subsequently causes testicular dysgenesis in Xenopus laevis, Ecotoxicol. Environ. Saf., 2022, 236, 113453.
84. S. C. Laws, S. A. Carey, J. M. Ferrell, G. J. Bodman and R. L. Cooper, Estrogenic activity of octylphenol, nonylphenol, bisphenol A and methoxychlor in rats, Toxicol Sci., 2000, 54(1), 154–167.
85. F. Webster, M. Gagné, G. Patlewicz, P. Pradeep, N. Trefiak and R. Judson, et al., Predicting estrogen receptor activation by a group of substituted phenols: An integrated approach to testing and assessment case study, Regul. Toxicol. Pharmacol., 2019, 106, 278–291.
86. D. L. Mobley and M. K. Gilson, Predicting Binding Free Energies: Frontiers and Benchmarks, Annu. Rev. Biophys., 2017, 46(46), 531–558.
87. I. Muegge and Y. Hu, Recent Advances in Alchemical Binding Free Energy Calculations for Drug Discovery, ACS Med. Chem. Lett., 2023, 14(3), 244–250.
88. Y. Fukunishi, Y. Yamashita, T. Mashimo and H. Nakamura, Prediction of Protein−compound Binding Energies from Known Activity Data: Docking-score-based Method and its Applications, Mol. Inf., 2018, 37(6–7), 1700120.
89. D. J. Kojetin and T. P. Burris, Small Molecule Modulation of Nuclear Receptor Conformational Dynamics: Implications for Function and Drug Discovery, Mol. Pharmacol., 2013, 83(1), 1–8.
90. M.-S. Park, C. Gao and H. A. Stern, Estimating binding affinities by docking/scoring methods using variable protonation states, Proteins:Struct., Funct., Bioinf., 2011, 79(1), 304–314.
91. P. H. M. Torres, A. C. R. Sodero, P. Jofily and F. P. Silva-Jr, Key Topics in Molecular Docking for Drug Design, Int. J. Mol. Sci., 2019, 20(18), 4574.
92. J. Wang, Q. Deng and L. Qi, Exploring the toxicological effects of DOTP exposure on periodontitis by combining molecular docking and molecular dynamics simulations, Sci. Rep., 2025, 15(1), 20915.
93. R. Yang, Y. Ye, Y. Chen, Y. Yang, L. Yang and Y. Yao, et al., First Insight into the Formation of In Vivo Transformation Products of 2-Ethylhexyl diphenyl phosphate in Zebrafish and Prediction of Their Potential Toxicities, Environ. Sci. Technol., 2023, 57(1), 451–462.
94. R. Cárdenas, J. Martínez-Seoane and C. Amero, Combining Experimental Data and Computational Methods for the Non-Computer Specialist, Molecules, 2020, 25(20), 4783.
95. OECD, Overview of Concepts and Available Guidance related to Integrated Approaches to Testing and Assessment (IATA), 2020.

---