Discovery of New Protein Targets of BPA Analogs and Derivatives Associated with Noncommunicable Diseases: A Virtual High-Throughput Screening.

BACKGROUND: Bisphenol A analogs and derivatives (BPs) have emerged as new contaminants with little or no information about their toxicity. These have been found in numerous everyday products, from thermal paper receipts to plastic containers, and measured in human samples.
OBJECTIVES: The objectives of this research were to identify in silico new protein targets of BPs associated with seven noncommunicable diseases (NCDs), and to study their protein-ligand interactions using computer-aided tools.
METHODS: Fifty BPs were identified by a literature search and submitted to a virtual high-throughput screening (vHTS) with 328 proteins associated with NCDs. Protein-protein interactions between predicted targets were examined using STRING, and the protocol was validated in terms of binding site recognition and correlation between in silico affinities and in vitro data.
RESULTS: According to the vHTS, several BPs may target proteins associated with NCDs, some of them with stronger affinities than bisphenol A (BPA). The best affinity score (the highest in silico affinity absolute value) was obtained after docking 4,4'-bis(N-carbamoyl-4-methylbenzensulfonamide)diphenylmethane (BTUM) on estradiol 17-beta-dehydrogenase 1 (-13.7 kcal/mol). However, other molecules, such as bisphenol A bis(diphenyl phosphate) (BDP), bisphenol PH (BPPH), and Pergafast 201 also exhibited great affinities (top 10 affinity scores for each disease) with proteins related to NCDs. DISCUSSION: Molecules such as BTUM, BDP, BPPH, and Pergafast 201 could be targeting key signaling pathways related to NCDs. These BPs should be prioritized for in vitro and in vivo toxicity testing and to further assess their possible role in the development of these diseases. https://doi.org/10.1289/EHP7466.

Introduction

Bisphenol A (BPA) is a high production-volume chemical used in the fabrication of polycarbonate plastics and epoxy resins and is a ubiquitous contaminant with endocrine-disrupting activity (Liao et al. 2012b; Takayanagi et al. 2006; Trasande et al. 2012). This compound is present in many everyday products (Santangeli et al. 2017; dos Santos Rosa 2018) and has been recently banned in certain goods, such as baby bottles, in some countries (EC 2011; Government of Canada 2010). This prohibition and the concern regarding its possible negative effects to human health have resulted in the increasing production of BPA analogs and derivatives (BPs) with similar structural or functional features (Lee et al. 2015; Liao et al. 2012d; Wang et al. 2017a).

Several BPs have been detected in different environmental matrices (Lee et al. 2015; Liao et al. 2012d; Yamazaki et al. 2015), everyday products (Catenza et al. 2021), thermal paper receipts (Liao et al. 2012c), food (Liao and Kannan 2013, 2014), indoor dust (Liao et al. 2012b), human breast milk (Niu et al. 2017), human blood serum (Owczarek et al. 2018), and human urinary samples (Hines et al. 2017; Liao et al. 2012a; Wang et al. 2019c). In addition, some BPs have been shown to be more bioaccumulated and biomagnified in the trophic chain than BPA due to their octanol–water coefficient (Wang et al. 2017a).

Some of these compounds can be present in products labeled as BPA-free (Inadera 2015; Rochester and Bolden 2015) because their use is not currently regulated or restricted. However, it is not clear if these are safer alternatives given that some of these molecules have induced similar effects as BPA in animal models (Rosenfeld 2017) and have presented endocrine-disrupting activity in vitro and in vivo, with even higher potency than BPA (Moreman et al. 2017; Rochester and Bolden 2015). Several human observational studies have suggested a possible association of BPs with diabetes (Duan et al. 2018), obesity (Andújar et al. 2019), oxidative stress (Kataria et al. 2017; Wang et al. 2019c; Zhang et al. 2016), neurodevelopmental effects (Jiang et al. 2020), and genotoxicity (Pelch et al. 2019). Furthermore, some of these compounds, such as bisphenol E (BPE) have been associated with negative effects after prenatal exposure in mice, such as inhibition of germ cell nest breakdown and a reduction of primary and secondary follicles (Shi et al. 2019). However, the available toxicological information about most of these emerging pollutants is very limited.

The study of the interactions of these chemicals with key proteins is crucial given that their molecular mechanisms are poorly understood (Sharma et al. 2018; Zhuang et al. 2014). Furthermore, it is a subject of emerging concern because there is epidemiologic evidence that the exposure to certain environmental pollutants, such as endocrine-disrupting chemicals (EDCs), might result in the development of noncommunicable diseases (NCDs) (Norman et al. 2013; Zarean and Poursafa 2019). NCDs are responsible for of all death worldwide, with an important contribution in health conditions such as cardiovascular disorders, cancers, respiratory diseases, and diabetes (WHO 2018b). These four groups of diseases account for of all premature death associated with NCDs (WHO 2018b). Therefore, the selection of the seven NCDs included in this study was made taking into account the health conditions highly related to these groups.

Computational chemistry approaches (such as molecular docking and molecular dynamics simulations) have been used to elucidate the estrogenicity of some BPs and to assess their ability to bind several nuclear receptors, such as the retinoid-related orphan nuclear receptors and the glucocorticoid receptor, by showing a correlation with in vitro data (Ng et al. 2015; Nishigori et al. 2012; Zhang et al. 2018). Therefore, the aim of this research was to identify in silico potential protein targets of BPs involved in the development of seven NCDs of high prevalence worldwide. These are cardiovascular diseases, lung cancer, breast cancer, cervical cancer, prostate cancer, diabetes, and thyroid disorders.

Methods

All calculations and data analysis described herein were performed on a Linux Ubuntu 18.04.5 LTS operating system. The system was run on a Dell Precision 3630 Tower workstation equipped with Intel Core i7-9700K CPU at 3.60 GHz (8 cores), 64 GB RAM, and GPU (NVIDIA Quadro P620 with 2 GB memory).

Data Extraction

Several BPs are new compounds (Xiong et al. 2020). Therefore, their names and structural features are only available in scientific articles. In order to create a list of them, a literature search was carried out using the following databases and search engines: Science Direct (http://www.sciencedirect.com/), PubMed (https://www.ncbi.nlm.nih.gov/pubmed/), and Google Scholar (https://scholar.google.com/) (Tober 2011), as well as the website of the U.S. Environmental Protection Agency (EPA; https://www.epa.gov/) and PubMed PubReMiner (https://hgserver2.amc.nl/cgi-bin/miner/miner2.cgi). The results were obtained using the following query: “BPA derivative” OR “BPA analog” OR “BPA analogue” OR “BPA substituent” OR “bisphenol A derivative” OR “bisphenol A analog” OR “bisphenol A analogue” OR “bisphenol A substituent”. Only articles in English were considered, and their abstracts were screened to recognize those reporting BPs. Each article was manually reviewed to identify the names and structures of the BPs.

On the other hand, human proteins associated with cardiovascular diseases, lung cancer, breast cancer, cervical cancer, prostate cancer, diabetes, and thyroid disorders were identified in the same databases and search engines by using generic queries with the following structure: the word “protein” AND “the name of the disease” (e.g., “protein” AND “breast cancer”). The selection criteria was the association of each protein with these diseases in at least one scientific article, which was manually verified.

In addition, we examined the human proteins related to each of the studied diseases deposited in the “gene–disease associations” section of the Comparative Toxicogenomics Database (CTD) (Davis et al. 2021) and selected those that also presented more than 500 hits in PubMed when searching for the name of the protein and the disease (e.g., “caspase-3” and “breast cancer”). To do that, a National Center for Biotechnology Information (NCBI) PubMed search was performed to retrieve the number of results (hits) for each protein and the seven NCDs included in this study (cardiovascular diseases, lung cancer, breast cancer, cervical cancer, prostate cancer, diabetes, and thyroid disorders), separately. Text files containing the lists of genes/proteins associated with each of the studied diseases obtained from the CTD were used as input to perform the searches on PubMed. Due to the large amount of data, the Entrez module and Bio.Entrez package of Biopython (Cock et al. 2009) were used with the following generic query structure: “protein name” AND “disease” to search across the lists of proteins associated with each of the NCDs included in this study. The number of results in PubMed generated by each protein–disease pair were recorded, and only those proteins with more than 500 hits were selected.

Preparation of the Ligand and Protein Structures

The three-dimensional (3D) structures of BPA and BPs with coordinates available in PubChem were downloaded from this database. The rest of the BPs were drawn using GaussView 3.0 program (Frisch et al. 2004). Subsequently, all the structures were optimized using the density functional theory (DFT) method at the B3LYP/6-31G level in Gaussian 09 (Frisch et al. 2009). The resultant geometries were translated to pdbqt files by Open Babel (version 2.3.1) (O’Boyle et al. 2011).

Selected human protein structures, with amino acids and resolution of , with 3D structures available in the Protein Data Bank (PDB), were downloaded from this database in pdb format (Berman et al. 2000). The crystallographic structures were prepared for molecular docking by removing all water molecules, ions, and other substructures, and by using the biopolymer structure preparation tool of Sybyl X-2.0 (Tripos) with default parameters. The optimization of the proteins was carried out in the same software by using the Powell method, as well as the Kollman United and Kollman all-atoms force fields with AMBER charges, , , , and (Montes-Grajales and Olivero-Verbel 2013). The resultant optimized structures were submitted to AutoDock Tools (MGL Tools) (Morris et al. 2009; Sanner 1999). This software was used to add Kollman charges and hydrogen polar atoms, as well as to convert them to pdbqt files. Furthermore, the size and coordinates of the center of the grid containing the whole protein surface were determined in the same program, with spacing of . These parameters were used in the molecular docking process.

Virtual High-Throughput Screening

The molecular docking calculations were carried out in AutoDock Vina 1.1.1 (Trott and Olson 2010) in triplicate with the following parameters: , , and . The means of the affinity scores resulting from the three replicates were calculated and used to rank the compounds in order to identify novel protein targets of BPs with strong theoretical binding affinity (high in silico affinity absolute value). A cutoff affinity of was used to distinguish the complexes with good affinity scores. This binding affinity estimate, which splits strong-binders and nonbinders or active and inactive compounds with a specificity of , was calculated in two of our previous articles (Montes-Grajales et al. 2018; Montes-Grajales and Olivero-Verbel 2020). The cutoff affinity of was obtained through a receiver operating characteristic (ROC) analysis by using docking affinity scores calculated with the same protocol used in this study and in vitro experimental activity information available online (Stojić et al. 2010; Wang et al. 2017c) of a similar data set of ligands and proteins. This data set comprised endocrine disruptors and some of the proteins included in this article, such as estrogen receptor alpha (ESR1), progesterone receptor (PRGR), androgen receptor (ANDR), retinoic acid receptor RXR-alpha (RXRA), and thyroid hormone receptor beta (THB) (Montes-Grajales et al. 2018; Montes-Grajales and Olivero-Verbel 2020).

Furthermore, heatmaps with clustering trees were generated by using the heatmap.2 function of R (version 3.6.3; R Development Core Team; Warnes et al. 2016), as reported in our previous article (Montes-Grajales and Olivero-Verbel 2020). The color key for the predicted affinity scores was established according to the affinity cutoff value of , calculated by ROC analysis in our previous publications (Montes-Grajales et al. 2018; Montes-Grajales and Olivero-Verbel 2020). Therefore, good affinity scores () were colored in red, weak affinities ( to ) in white, and the predicted nonbinders () in blue. The R code to generate the heatmaps is available in the Supplemental Materials (“Heatmaps_Script.R”).

Protein–Ligand Interactions

The interactions between BPs and proteins related to NCDs that exhibited the best affinity scores (top 10 affinity scores for each disease) were determined in silico using the LigandScout 3.0 program (Wolber and Langer 2005) with default parameters. The best ligand pose resultant from the molecular docking with AutoDock Vina 1.1.1. (Trott and Olson 2010) was extracted by using AutoDock Tools (Morris et al. 2009; Sanner 1999), translated to pdb format, and merged with the optimized protein in Sybyl X-2.0 (Tripos). Each protein–ligand complex was opened in LigandScout 3.0. The ligand was selected from the graphic interface of the software and analyzed using the structure-based tool to determine the contact residues of the protein interacting with the ligand, as well as the nature of these interactions (hydrophobic interactions, aromatic ring interactions, or hydrogen bonds). The retrieved information was used to create images of the 3D view of the protein–ligand complexes with the highest in silico affinity (absolute value) for each disease resultant from the virtual high-throughput screening (vHTS), as well as their interactions.

Protein–Protein Interaction Network

The functional association between the proteins that presented the best affinity scores (top 10 affinity scores for each disease) with BPs were carried out using STRING 11.0. (Szklarczyk et al. 2019). The short names of the proteins belonging to the protein–ligand complexes with the top 10 affinity scores with BPs were used as query in the search box for “Multiple proteins” of STRING 11.0. (Szklarczyk et al. 2019). The parameter used were as follows: Species: Homo sapiens; Network type: full STRING network; required score: high confidence: 0.700; maximum number of interactors to show: none/query proteins only; meaning of network edges: by evidence; and active interaction sources: text mining, neighborhood, experiments, databases, and co-occurrence. The software establishes protein–protein interactions from multiple sources based on the relationship between the input macromolecules and their reported biological processes. The interactions are determined using pathway knowledge, experimental data, and text mining in biological databases, among others (Szklarczyk et al. 2019).

Validation and Molecular Dynamics Simulation

In vitro experimental values of half-maximal activity concentration () of 55 protein–ligand complexes containing BPs with proteins included in the data set of the vHTS [ESR1, ESR2, peroxisome proliferator-activated receptor gamma (PPARG), peroxisome proliferator-activated receptor delta (PPARD), vitamin D3 receptor (VDR), and thyroid-stimulating hormone receptor (TSHR)] were obtained from PubChem Bioassay (Wang et al. 2017c) (see Table 9), and used to perform a correlation with the calculated molecular docking affinity scores for the same protein–ligand complexes. The Pearson’s correlation coefficient and -value were determined to evaluate the association between the molecular docking results and in vitro data. Furthermore, the validation of our protocol to assess the power of prediction to correctly determine the binding site was carried out. The crystallographic structure of the complexes bisphenol B (BPB)/estrogen-related receptor gamma [ERR3; PDB identification (ID): 1I61] and BPE/ERR3 (PDB ID: 6I64) were compared with the resultant structures from the molecular docking. Molecular docking simulations were carried out following the protocols used in this study for protein and ligand preparation and for vHTS. The root mean square deviation (RMSD) was calculated by using DockRMSD (Bell and Zhang 2019).

Table 9

Data set of compounds and proteins used for validation showing the calculated molecular docking affinity scores and experimental values of half-maximal activity concentration () obtained from PubChem BioAssay (Wang et al. 2017c).

BPA analogs and derivatives	Protein	Affinity [(kcal/mol) mean±SD]	AC_50 (μM )	PubChem bioassay	Reference
BPB	ESR1	−9.0	0.138	743077	NCBI 2014b
BPAF	ESR1	−9.3	0.373	743077	NCBI 2014b
BPA	ESR1	−9.8	1.127	743077	NCBI 2014b
BPE	ESR1	−9.5	4.249	743077	NCBI 2014b
BPZ	ESR1	−9.3	6.813	743077	NCBI 2014b
4-CP	ESR1	−9.9	7.128	743077	NCBI 2014b
BPS	ESR1	−8.5	9.106	743077	NCBI 2014b
BPF	ESR1	−8.2	22.038	743077	NCBI 2014b
BzP	ESR1	−8.8	23.891	743077	NCBI 2014b
2,2-BPF	ESR1	−8.3	27.539	743075	NCBI 2014a
TBBPA	ESR1	−6.8	43.492	743077	NCBI 2014b
BFDGE	ESR1	−7.1	43.828	743075	NCBI 2014a
DBP	ESR1	−6.3	49.002	743077	NCBI 2014b
TCBPA	ESR1	−7.1	51.880	743077	NCBI 2014b
TGSA	ESR1	−7.1	54.941	743075	NCBI 2014a
BPFL	ESR1	−7.6	54.954	743077	NCBI 2014b
BADGE	ESR1	−7.2	62.460	743075	NCBI 2014a
BPAF	ESR1	−9.3	0.171	743079	NCBI 2014c
BPE	ESR1	−9.5	1.812	743079	NCBI 2014c
BPC	ESR1	−8.4	0.243	743079	NCBI 2014c
BADGE	ESR2	−6.9	61.940	1259380	NCBI 2018a
TBBPA	ESR2	−6.2	68.311	1259380	NCBI 2018a
DMBPA	ESR2	−7.6	29.179	1259380	NCBI 2018a
BPAF	ESR2	−7.0	60.820	1259382	NCBI 2018b
BPC	ESR2	−7.0	68.450	1259380	NCBI 2018a
BFDGE	ESR2	−6.8	48.558	1259380	NCBI 2018a
BPZ	ESR2	−6.6	57.912	1259380	NCBI 2018a
BPE	ESR2	−8.2	10.871	1259396	NCBI 2018c
BPAP	ESR2	−6.7	61.131	1259382	NCBI 2018b
BPPH	ESR2	−8.4	27.306	1259382	NCBI 2018b
TBBPA	PPARG	−7.4	43.545	743194	NCBI 2014e
TCBPA	PPARG	−7.7	24.541	743194	NCBI 2014e
BPB	PPARG	−7.5	30.869	743191	NCBI 2014d
BPAF	PPARG	−7.6	43.647	743194	NCBI 2014e
2,2-BPF	PPARG	−7.1	30.899	743191	NCBI 2014d
BPFL	PPARG	−8.0	17.374	743194	NCBI 2014e
TMBPA	PPARG	−8.2	29.882	743191	NCBI 2014d
BPZ	PPARG	−7.8	41.345	743194	NCBI 2014e
BPE	PPARG	−7.3	36.462	743191	NCBI 2014d
BADGE	PPARD	−8.8	34.670	743211	NCBI 2014f
BPB	PPARD	−7.4	54.894	743211	NCBI 2014f
BPAF	PPARD	−7.9	54.948	743211	NCBI 2014f
BPFL	PPARD	−9.3	10.962	743215	NCBI 2014g
TMBPA	PPARD	−8.1	47.359	743211	NCBI 2014f
BPZ	PPARD	−8.4	41.345	743211	NCBI 2014f
BPE	PPARD	−7.6	40.911	743215	NCBI 2014g
TGSA	PPARD	−8.0	43.641	743211	NCBI 2014f
BPC	PPARD	−7.8	19.456	743215	NCBI 2014g
BPZ	VDR	−8.7	36.849	743225	NCBI 2014h
BPAF	VDR	−8.8	38.900	743225	NCBI 2014h
DMBPA	TSHR	−6.2	46.246	1224843	NCBI 2016a
BPAF	TSHR	−6.1	38.375	1224895	NCBI 2016b
TMBPA	TSHR	−6.7	33.226	1224843	NCBI 2016a
BPFL	TSHR	−7.1	54.483	1224843	NCBI 2016a
BPC	TSHR	−5.6	54.372	1224895	NCBI 2016b

Note: Experimental values of half-maximal activity concentration () were obtained from PubChem BioAssay. BADGE, bisphenol A diglycidyl ether; BFDGE, bisphenol F diglycidyl ether; BP, bisphenol A derivative; BPA, bisphenol A; BPAF, bisphenol AF; BPAP, bisphenol AP; BPB, bisphenol B; BPC, bisphenol C; BPE, bisphenol E; BPF, bisphenol F; BPFL, bisphenol FL; BPPH, bisphenol PH; BPS, bisphenol S; BzP, benzylparaben; DBP, dibutyl phthalate; DMBPA, 3,3’-dimethylbisphenol A; ESR1, estrogen receptor alpha; ESR2, estrogen receptor beta; PPARD, peroxisome proliferator activated receptor delta; PPARG, peroxisome proliferator activated receptor gamma; SD, standard deviation; TBBPA, tetrabromobisphenol A; TCBPA, tetrachlorobisphenol A; TGSA, 2,2’-diallyl-4,4’-sulfonyldiphenol; TMBPA, tetramethylbisphenol A; TSHR, thyroid-stimulating hormone receptor; VDR, vitamin D3 receptor; 2,2-BPF, 2,2’-bisphenol F; 4-CP, 4-cumylphenol.

For comparison, the affinity scores of 10 proteins with co-crystallized ligands available in PDB were calculated. The proteins were randomly selected from the top 10 protein–ligand complexes that presented the best affinity scores for each of the studied NCDs (top 10 affinity scores for each disease). The co-crystallized ligands were extracted from the pdb file and employed for molecular docking studies using the same protocol used for the vHTS.

In addition, a molecular dynamics (MD) simulation of the protein–ligand complex with the highest in silico affinity (absolute value) was performed using Gromacs (version 2020.2) (Abraham et al. 2015) to confirm the stability of the system. The Chemistry at Harvard Macromolecular Mechanics (CHARMM) General Force Field (CGenFF) (Vanommeslaeghe et al. 2010) and the CHARMM force field (MacKerell et al. 1998) were used for the ligand and protein, respectively. The protein–ligand complex was solvated by inserting it into the center of a cubic box filled with water, from the edges of the complex. Subsequently, ions were added to neutralize the system. Constant pressure (NVT) equilibrium was performed for 1 ns with a time step of 2 fs and reference temperature of 300 K, a second equilibrium step was carried out by using a constant particle number, pressure, and temperature (NPT) ensemble for 1 ns. The production step of the MD simulation was performed during 10 ns under isothermal–isobaric conditions, with a time step of , reference temperature of , pressure of , van der Waals cutoff of , and grid spacing of using the leap-frog integrator and Verlet cutoff scheme. The atomic coordinates and velocities of the systems were recorded every 10 ps.

Results

A total of 328 human proteins associated with cardiovascular diseases, lung cancer, breast cancer, cervical cancer, prostate cancer, diabetes, and thyroid disorders were identified by a literature search using several databases and search engines, such as Science Direct (http://www.sciencedirect.com/), PubMed (https://www.ncbi.nlm.nih.gov/pubmed/), and Google Scholar (https://scholar.google.com/) (Table 1). These proteins have been associated with these NCDs in scientific reports or the CTD (Excel Table S1), and have 3D structures available in PDB (Berman et al. 2000). In addition, 50 BPs were found in scientific articles (Excel Table S2). The chemical structures of these molecules are available in Figure S1.

Table 1

Number of proteins associated with noncommunicable diseases (NCDs) selected by literature search.

NCDs	Number of proteins
Cardiovascular diseases	84
Lung cancer	38
Breast cancer	75
Cervical cancer	29
Prostate cancer	34
Diabetes	43
Thyroid disorders	25

Note: Human proteins associated with seven NCDs were identified through a literature search using Science Direct (http://www.sciencedirect.com/), PubMed (https://www.ncbi.nlm.nih.gov/pubmed/), Google Scholar (https://scholar.google.com/), the website of the U.S. Environmental Protection Agency (EPA; https://www.epa.gov/), and PubMed PubReMiner (https://hgserver2.amc.nl/cgi-bin/miner/miner2.cgi). The queries included the word “protein” and the name of the disease (e.g., “protein” and “breast cancer”).

vHTS and Protein–Ligand Interactions

In order to identify new protein targets of BPs associated with NCDs, a total of 16,728 different protein–ligand pairs were assessed through a vHTS. The in silico affinity scores (mean of triplicates) obtained from the molecular docking simulations between BPA and 50 of its analogs and derivatives with 328 proteins associated with seven prevalent NCDs are shown in Excel Table S3. This study presents BPs that exhibited high in silico affinity absolute values for the proteins associated with cardiovascular diseases, lung cancer, breast cancer, cervical cancer, prostate cancer, diabetes, and thyroid disorders as compounds to be prioritized for in vitro and in vivo testing.

Cardiovascular Diseases

The in silico docking affinity scores of proteins related to cardiovascular diseases with BPs ranged from to (Excel Table S3). The best affinity score (the highest in silico affinity absolute value) in this group was obtained for 4,4′-bis(N-carbamoyl-4-methylbenzensulfonamide)diphenylmethane (BTUM) with pro protein convertase subtilisin/kexin type 9 (PCSK9) (Figure 1). The BPs that exhibited the best affinity scores (top 10 affinity scores) for proteins related to cardiovascular diseases were BTUM, bisphenol A bis(diphenyl phosphate) (BDP), bisphenol PH (BPPH), and Pergafast 201 (Table 2).

Figure 1.

(A) Three-dimensional view of the 4,4′-bis(N-carbamoyl-4-methylbenzensulfonamide)diphenylmethane/proprotein convertase subtilisin/kexin type 9 (BTUM/PCSK9) complex, showing (B) the binding site and interactions predicted by LigandScout 3.1. Contact residues: VAL359B, CYS358B, ILE416B, THR459B, VAL460B, ALA475B, VAL655B, THR653B, VAL650B, and THR623B. The black arrows with circles represent hydrogen-bond donor features.

Table 2

Top 10 docking affinity scores of bisphenol A (BPA) analogs and derivatives (BPs) targeting proteins related to cardiovascular diseases.

BPA analogs and derivatives	Proteins	Short names	PDB ID	Affinity [(kcal/mol) mean±SD]
BTUM	Proprotein convertase subtilisin/kexin type 9a	PCSK9	6U26	−12.6±0.3
BTUM	Fibroblast growth factor receptor 3b	FGFR3	4K33	−12.5±0.1
BTUM	Nitric oxide synthase, endothelial	NOS3	4D1P	−12.5±0.2
BDP	Nitric oxide synthase, endothelial	NOS3	4D1P	−12.1±0.2
BDP	Cholesteryl ester transfer proteina	CETP	2OBD	−12.1±0.0
BTUM	Cholesteryl ester transfer proteina	CETP	2OBD	−12.0±0.4
BTUM	Angiotensin-converting enzymec	ACE	6H5W	−11.9±0.1
BTUM	15-Hydroxyprostaglandin dehydrogenase [NAD^+]	PGDH	2GDZ	−11.9±0.2
BPPH	Carbamoyl-phosphate synthase [ammonia], mitochondrial	CPSM	5DOT	−11.7±0.0
Pergafast 201	Nitric oxide synthase, endothelial	NOS3	4D1P	−11.6±0.4

Note: Molecular docking calculations were performed with Autodock Vina 1.1.1. BDP, bisphenol A bis(diphenyl phosphate); BPPH, bisphenol PH; BTUM, 4,4′-bis(N-carbamoyl-4-methylbenzensulfonamide)diphenylmethane; , nicotinamide adenine dinucleotide; PDB ID, Protein Data Bank identification; SD, standard deviation.

Proteins also associated with diabetes.

Protein also associated with diabetes and lung cancer.

Protein also associated with breast cancer, diabetes, lung cancer, and prostate cancer.

In addition, complexes formed by other BPs, such as bisphenol FL (BPFL), also presented strong affinity scores for several targets. Among them are the PPARG, E3 ubiquitin–protein ligase MIB1 (MIB1), lysine-specific demethylase 6A (KDM6A), desmoplakin (DESP), 15-hydroxyprostaglandin dehydrogenase [] (PGDH), fibroblast growth factor receptor 3 (FGFR3), lysosomal acid glucosylceramidase (GLCM), carbamoyl-phosphate synthase [ammonia] mitochondrial (CPSM), folate hydrolase 1 (FOLH1), cholesteryl ester transfer protein (CETP), cyclin-dependent kinase 13 (CDK13), lipoprotein lipase (LIPL), neurogenic locus notch homolog protein 3 (NOTC3), and calcium/calmodulin-dependent protein kinase type II delta chain (CAMK2D). A heatmap with dendrograms showing the docking affinity scores predicted between BPs with proteins of this category is shown in Figure S2.

Lung Cancer

BPs had the potential to interact in silico with a broad range of proteins associated with lung cancer (Excel Table S3). The complex that presented the best affinity score (the highest in silico affinity absolute value) in this group was formed by Pergafast 201 with l-lactate dehydrogenase A chain (LDHA), which exhibited an affinity score of (Figure 2).

Figure 2.

(A) Three-dimensional view of the 3-(3-tosylureido)phenyl -toluenesulfonate/l-lactate dehydrogenase A chain (Pergafast 201/LDHA) complex, showing (B) the binding site and interactions predicted by LigandScout 3.1. Contact residues: TRP187B, ARG170B, ALA167B, LEU69A, PHE70A, LEU182B, ARG268D, LEU182D, TRP187D, LEU69C, ALA167D, and ARG170D. The red arrows represent hydrogen-bond acceptor features, and the blue double-sided arrows symbolize aromatic ring interactions.

However, other complexes also showed high in silico affinity absolute values (Table 3), with a predominant occurrence of Pergafast 201, BTUM, and BDP. The complete set of results of BPs that interacted in silico with proteins involved in lung cancer is better visualized in Figure S3, through a heatmap with clustering trees.

Table 3

Top 10 docking affinity scores of bisphenol A (BPA) analogs and derivatives (BPs) targeting proteins related to lung cancer.

BPA analogs and derivatives	Proteins	Short names	PDB ID	Affinity [(kcal/mol) mean±SD]
Pergafast 201	l-Lactate dehydrogenase A chain	LDHA	5W8J	−11.9±0.1
Pergafast 201	Fructose-bisphosphate aldolase A	ALDOA	5KY6	−11.7±0.3
BTUM	Receptor of activated protein C kinase 1	RACK1	4AOW	−11.6±0.1
BDP	Fructose-bisphosphate aldolase A	ALDOA	5KY6	−11.6±0.3
BTUM	l-Lactate dehydrogenase A chain	LDHA	5W8J	−11.4±0.1
BDP	Arf-GAP with SH3 domain, ANK repeat, and PH domain-containing protein 3	ASAP3	2B0O	−11.3±0.3
BTUM	Epidermal growth factor receptora	EGFR	3POZ	−11.2±0.1
BDP	Epidermal growth factor receptora	EGFR	3POZ	−10.9±0.0
Pergafast 201	Arf-GAP with SH3 domain, ANK repeat, and PH domain-containing protein 3	ASAP3	2B0O	−10.9±0.3
BTUM	Arf-GAP with SH3 domain, ANK repeat, and PH domain-containing protein 3	ASAP3	2B0O	−10.7±0.1
Pergafast 201	GTPase Kras	KRAS	6P0Z	−10.7±0.1

Note: Molecular docking calculations were performed with Autodock Vina 1.1.1. BDP, bisphenol A bis(diphenyl phosphate); BTUM, 4,4′-bis(N-carbamoyl-4-methylbenzensulfonamide)diphenylmethane; PDB ID, Protein Data Bank identification; SD, standard deviation.

Proteins also associated with cardiovascular disease, breast cancer, cervical cancer, prostate cancer, and diabetes.

Breast Cancer

BPs had a strong potential to interact in silico with proteins associated with breast cancer (Excel Table S3). The best affinity score (the highest in silico affinity absolute value) in this group () was obtained for BTUM with estradiol 17-beta-dehydrogenase 1 (DHB1; Figure 3).

Figure 3.

(A) Three-dimensional view of the 4,4′-bis(N-carbamoyl-4-methylbenzensulfonamide)diphenylmethane/estradiol 17-beta-dehydrogenase 1 (BTUM/DHB1) complex, showing (B) the binding site and interactions predicted by LigandScout 3.1. Contact residues: VAL225X, PHE259X, LEU149X, VAL143X, MET147X, GLY144X, SER142X, PHE192X, VAL188X, ILE14X, THR190X, ARG37X, THR140X, ALA91X, VAL113X, and VAL66X.

The protein–ligand complexes with the top 10 docking affinity scores between BPs with proteins related to breast cancer are presented in Table 4. However, many other complexes also showed good affinity scores () (Figure S4).

Table 4

Top 10 docking affinity scores of bisphenol A (BPA) analogs and derivatives (BPs) targeting proteins related to breast cancer.

BPA analogs and derivatives	Proteins	Short names	PDB ID	Affinity [(kcal/mol) mean±SD]
BTUM	Estradiol 17-beta-dehydrogenase 1	DHB1	3HB5	−13.7±0.0
BTUM	Cytochrome P450 2D6	CP2D6	3TBG	−13.5±0.4
BDP	Nuclear receptor ROR-gamma	RORG	3L0L	−12.1±0.3
BTUM	Serine/threonine-protein kinase Chk2	CHK2	2W0J	−12.0±0.0
BTUM	Nuclear receptor ROR-gamma	RORG	3L0L	−11.9±0.1
Pergafast 201	Estradiol 17-beta-dehydrogenase 1	DHB1	3HB5	−11.9±0.3
Pergafast 201	Cytochrome P450 2D6	CP2D6	3TBG	−11.9±0.1
Pergafast 201	Nuclear receptor ROR-alpha	RORA	1N83	−11.8±0.0
Pergafast 201	NAD(P)H dehydrogenase [quinone] 1	NQO1	1D4A	−11.8±0.3
Pergafast 201	Breast cancer type 1 susceptibility protein	BRCA1	3COJ	−11.7±0.0

Note: Molecular docking calculations were performed with Autodock Vina 1.1.1. BDP, bisphenol A bis(diphenyl phosphate); BTUM, 4,4′-bis(N-carbamoyl-4-methylbenzensulfonamide)diphenylmethane; NAD(P)H, nicotinamide adenine dinucleotide phosphate; PDB ID, Protein Data Bank identification; ROR, retinoic acid receptor-related orphan receptor; SD, standard deviation.

Cervical Cancer

BPs were predicted to interact with multiple proteins related to cervical cancer. The results of the molecular docking between these compounds and proteins associated with cervical cancer are shown in Excel Table S3. The best affinity score (the highest in silico affinity absolute value) in this group was obtained for the BDP/NAD-dependent protein deacetylase sirtuin-2 (SIR2) complex (; Figure 4). According to the vHTS, this BPA derivative bound the extended C-site of the protein (Rumpf et al. 2015) and protruded into the acetyl–lysine channel and the selectivity pocket through its interaction with the following contact residues: PHE96A, PHE119A, TYR139A, ALA135A, PHE143A, LEU206A, PHE190A, ILE232A, PHE235A, and PHE131A. The top 10 complexes that obtained the best affinity scores (the highest in silico affinity absolute values) in this group are presented in Table 5. The protein SIR2 seems to be a common target in this group for BPs (Table 5). However, other proteins such as the proto-oncogene tyrosine-protein kinase Src (SRC), methylenetetrahydrofolate reductase (MTHR), complement factor I (CFAI), phosphoinositide-3-kinase catalytic alpha polypeptide (PK3CA), nonreceptor tyrosine-protein kinase TYK2 (TYK2), aurora kinase B (AURKB), eyes absent homolog 2 (EYA2), NAD-dependent protein deacetylase sirtuin-7 (SIR7), and alpha-albumin (AFAM) also showed good affinity scores () for these compounds (Figure S5).

Figure 4.

(A) Three-dimensional view of the bisphenol A bis(diphenyl phosphate)/NAD-dependent protein deacetylase sirtuin-2 (BDP/SIR2) complex, showing (B) the binding site and interactions predicted by LigandScout 3.1. Contact residues: PHE143A, LEU206A, PHE190A, PHE96A, ILE232A, PHE235A, PHE131A, PHE119A, TYR104A, LEU103A, LEU134A, ALA135A, LEU138A, TYR139A, and ILE93A. The blue double-sided arrows represent aromatic ring interactions. Note: NAD, nicotinamide adenine dinucleotide.

Table 5

Top 10 docking affinity scores of bisphenol A (BPA) analogs and derivatives (BPs) targeting proteins related to cervical cancer.

BPA analogs and derivatives	Proteins	Short names	PDB ID	Affinity [(kcal/mol) mean±SD]
BDP	NAD-dependent protein deacetylase sirtuin-2	SIR2	4RMH	−13.1±0.3
BPP	NAD-dependent protein deacetylase sirtuin-2	SIR2	4RMH	−12.3±0.0
BTUM	NAD-dependent protein deacetylase sirtuin-2	SIR2	4RMH	−12.3±0.2
Pergafast 201	NAD-dependent protein deacetylase sirtuin-2	SIR2	4RMH	−12.1±0.2
BTUM	Eyes absent homolog 2	EYA2	4EGC	−11.7±0.0
BPPH	NAD-dependent protein deacetylase sirtuin-2	SIR2	4RMH	−11.2±0.0
BDP	Eyes absent homolog 2	EYA2	4EGC	−11.0±0.1
BTUM	Proto-oncogene tyrosine-protein kinase Src	SRC	1FMK	−11.0±0.1
BPM	NAD-dependent protein deacetylase sirtuin-2	SIR2	4RMH	−10.9±0.0
BPPH	NAD-dependent protein deacetylase sirtuin-7	SIR7	5IQZ	−10.9±0.0

Note: Molecular docking calculations were performed with Autodock Vina 1.1.1. BDP, bisphenol A bis(diphenyl phosphate); BPM, Bisphenol M; BPPH, bisphenol PH; BTUM, 4,4′-bis(N-carbamoyl-4-methylbenzensulfonamide)diphenylmethane; NAD, nicotinamide adenine dinucleotide; ROR, retinoic acid receptor-related orphan receptor; PDB ID, Protein Data Bank identification; SD, standard deviation.

Prostate Cancer

BPs were predicted to interact with several proteins associated with prostate cancer. The molecular docking affinities resultant from the vHTS are shown in Excel Table S3. The complex with the best docking affinity score (the highest in silico affinity absolute value) in this group was BPPH/poly [ADP-ribose] polymerase 1 (PARP1) with (Figure 5). Complexes formed by other derivatives and analogs of the plasticizer BPA also presented high affinity scores (absolute values) (Table 6) and included Pergafast 201, BDP, BPFL, and BTUM.

Figure 5.

(A) Three-dimensional view of the bisphenol PH/poly [ADP-ribose] polymerase 1 (BPPH/PARP1) complex, showing (B) the binding site and interactions predicted by LigandScout 3.1. Contact residues: TYR907B, ILE872B, GLY863B, ALA898B, TYR896B, ALA880B, and TYR889B. The black arrows with circles represent hydrogen-bond donor features. Note: ADP, adenine diphosphate.

Table 6

Top 10 docking affinity scores of bisphenol A (BPA) analogs and derivatives (BPs) targeting proteins related to prostate cancer.

BPA analogs and derivatives	Proteins	Short names	PDB ID	Affinity [(kcal/mol) mean±SD]
BPPH	Poly [ADP-ribose] polymerase 1a	PARP1	5WS1	−12.6±0.0
Pergafast 201	Poly [ADP-ribose] polymerase 1a	PARP1	5WS1	−11.6±0.4
BDP	SRSF protein kinase 1	SRPK1	5MY8	−11.2±0.2
Pergafast 201	Flavin-containing amine oxidase domain-containing protein 2	KDM1A	2DW4	−11.1±0.2
BPFL	Poly [ADP-ribose] polymerase 1a	PARP1	5WS1	−11.0±0.0
Pergafast 201	Peroxiredoxin-1b	PRDX1	4XCS	−11.0±0.3
Pergafast 201	Sister chromatid cohesion protein PDS5 homolog B	PDS5B	5HDT	−11.0±0.2
BTUM	Prostatic acid phosphatase	PPAP	1ND6	−10.9±0.1
BDP	Flavin-containing amine oxidase domain-containing protein 2	KDM1A	2DW4	−10.9±0.3
BDP	Poly [ADP-ribose] polymerase 1a	PARP1	5WS1	−10.9±0.5

Note: Molecular docking calculations were performed with Autodock Vina 1.1.1. ADP, adenine diphosphate; BDP, bisphenol A bis(diphenyl phosphate); BPFL, bisphenol FL; BPPH, bisphenol PH; BTUM, 4,4′-bis(N-carbamoyl-4-methylbenzensulfonamide)diphenylmethane; PDB ID, Protein Data Bank identification; SD, standard deviation; SRPK, serine–arginine protein kinase.

Proteins also associated with lung cancer, breast cancer, cervical cancer, and diabetes.

Proteins also associated with lung cancer, breast cancer, and diabetes.

A heatmap representing the molecular docking affinity scores is shown in Figure S6. Numerous BPs presented high binding affinities (absolute values) for proteins associated with prostate cancer, such as serine-arginine protein kinase 1 (SRPK1), polycomb protein embryonic ectoderm development (EED), and PARP1. Furthermore, several of these small molecules exhibited a multi-target behavior with strong affinities for numerous proteins, among them were the BPs BTUM and Pergafast 201.

Diabetes

The results of the vHTS between BPs with proteins related to diabetes are presented in Excel Table S3. The best affinity score (the highest in silico affinity absolute value) in this group was obtained for the BTUM/bile salt-activated lipase [or carboxyl ester lipase (CEL)] complex (affinity: ; Figure 6). According to the vHTS, BTUM interacted with the catalytic site of CEL (Touvrey et al. 2019) (PDB ID: 6H0T) through the following contact residues: TYR123A, ALA108A, and VAL285A. However, other protein–ligand pairs such as those presented in Table 7 also exhibited strong affinities in silico.

Figure 6.

(A) Three-dimensional view of the 4,4′-bis(N-carbamoyl-4-methylbenzensulfonamide)diphenylmethane/bile salt-activated lipase (BTUM/CEL) complex, showing (B) the binding site and interactions predicted by LigandScout 3.1. Contact residues: LEU124A, THR140A, TYR123A, TYR105A, ASN84A, PHE60A, ASN142A, VAL145A, ALA108A, VAL288A, and VAL285A. The red arrows represent hydrogen-bond acceptor features, black arrows with circles represent hydrogen-bond donor features, and blue double-sided arrows symbolize aromatic ring interactions.

Table 7

Top 10 docking affinity scores of bisphenol A (BPA) analogs and derivatives (BPs) targeting proteins related to diabetes.

BPA analogs and derivatives	Proteins	Short names	PDB ID	Affinity [(kcal/mol) mean±SD]
BTUM	Bile salt-activated lipase	CEL	6H0T	−13.2±0.0
Pergafast 201	Aldo-keto reductase family 1 member B1	ALDR	1ADS	−11.7±0.4
BTUM	Angiopoietin-related protein 3a	ANGL3	6EUA	−11.6±0.3
BPPH	Angiopoietin-related protein 3a	ANGL3	6EUA	−11.6±0.0
Pergafast 201	Peroxisome proliferator-activated receptor gamma	PPARG	1ZGY	−11.5±0.0
BTUM	Aldo-keto reductase family 1 member B1	ALDR	1ADS	−11.4±0.3
BTUM	Cytochrome b5 reductase 4	NB5R4	6MV2	−11.0±0.1
Pergafast 201	Angiopoietin-related protein 3a	ANGL3	6EUA	−11.0±0.3
Pergafast 201	Bile salt-activated lipase	CEL	6H0T	−10.9±0.1
Pergafast 201	Fructose-1,6-bisphosphatase 1	F16P1	1FTA	−10.8±0.1

Note: Molecular docking calculations were performed with Autodock Vina 1.1.1. BPPH, bisphenol PH; BTUM, 4,4′-bis(N-carbamoyl-4-methylbenzensulfonamide)diphenylmethane; PDB ID, Protein Data Bank identification; SD, standard deviation.

Proteins also associated with cardiovascular disease.

According to the heatmap and clustering trees (Figure S7), several BPs presented a promiscuous behavior targeting numerous proteins with strong affinities. Among them are BDP, BPPH, and Pergafast 201.

Thyroid Disorders

The docking affinity scores resultant of the vHTS between BPs with proteins related to thyroid disorders are presented in Excel Table S3. The best affinity score (the highest in silico affinity absolute value) in this group was obtained for BTUM/NAD-dependent protein deacetylase sirtuin-6 (SIR6) complex (; Figure 7). However, BPs had the potential to target many other proteins with high affinity (Table 8), such as the kelch-like erythroid cell-derived protein with CNC homology (ECH)-associated protein 1 (KEAP1), mitogen-activated protein kinase 3 (MK03), merlin (MERL), thyroid hormone receptor alpha (THA), THB, and E3 ubiquitin–protein ligase TRIM33 (TRI33). Affinity scores obtained for BPs with proteins related to thyroid disorders are better visualized in the heatmap with clustering trees presented in Figure S8.

Figure 7.

(A) Three-dimensional view of the 4,4′-bis(N-carbamoyl-4-methylbenzensulfonamide)diphenylmethane/NAD-dependent protein deacetylase sirtuin-6 (BTUM/SIR6) complex, showing (B) the binding site and interactions predicted by LigandScout 3.1. Contact residues: ASN240A, VAL258A, LEU241A, TYR257A, GLN242A, ALA58A, THR57A, PHE64A, ALA53A, ILE219A, LEU186A, and TRP188A. The red arrows represent hydrogen-bond acceptor features. Note: NAD, nicotinamide adenine dinucleotide.

Table 8

Top 10 docking affinity scores of bisphenol A (BPA) analogs and derivatives (BPs) targeting proteins related to thyroid disorders.

BPA analogs and derivatives	Proteins	Short names	PDB ID	Affinity [(kcal/mol) mean±SD]
BTUM	NAD-dependent protein deacetylase sirtuin-6	SIR6	5MF6	−12.1±0.3
BTUM	Kelch-like ECH-associated protein 1a	KEAP1	1ZGK	−11.4±0.3
BDP	Mitogen-activated protein kinase 3b	MK03	4QTB	−11.1±0.2
Pergafast 201	Kelch-like ECH-associated protein 1a	KEAP1	1ZGK	−11.0±0.3
BTUM	Merlina	MERL	1H4R	−10.9±0.1
BPM	Thyroid hormone receptor alphac	THA	3ILZ	−10.9±0.1
Pergafast 201	Mitogen-activated protein kinase 3b	MK03	4QTB	−10.9±0.3
BTUM	Thyroid hormone receptor betad	THB	1N46	−10.8±0.1
BTUM	Mitogen-activated protein kinase 3b	MK03	4QTB	−10.8±0.2
BTUM	E3 ubiquitin–protein ligase TRIM33	TRI33	3U5N	−10.7±0.1

Note: Molecular docking calculations were performed with Autodock Vina 1.1.1. BDP, bisphenol A bis(diphenyl phosphate); BPM, Bisphenol M; BTUM, 4,4′-bis(N-carbamoyl-4-methylbenzensulfonamide)diphenylmethane; ECH, erythroid cell-derived protein with CNC homology; NAD, nicotinamide adenine dinucleotide; PDB ID, Protein Data Bank identification; SD, standard deviation.

Proteins also associated with cardiovascular disease, lung cancer, breast cancer, and diabetes.

Proteins also associated with cardiovascular disease, lung cancer, breast cancer, cervical cancer, diabetes, and prostate cancer.

Proteins also associated with diabetes.

Proteins also associated with breast cancer and diabetes.

The protein–protein interaction network developed using STRING 11.0. (Szklarczyk et al. 2019), based on their functional associations, showed that numerous theoretical targets of BPs involved in NCDs (Tables 2–8) were interrelated. Twenty-one of them presented protein–protein associations in the cluster analysis. Furthermore, some of these proteins, such as epidermal growth factor receptor (EGFR), SIR2, PARP1, PPARG, and SRC represented hubs in this network. The two targets with the highest number of protein–protein interactions were EGRF and SIR2, both with five direct interactors, followed by PPARG and SRC with four protein–protein associations for each one. Most of the protein–protein associations have been confirmed experimentally according to the results of STRING 11.0. (Szklarczyk et al. 2019) (Figure S9).

A two-step validation was carried out to evaluate the association between the calculated protein–ligand affinity scores and in vitro data, as well as to assess the accuracy of the docking pose prediction compared with the corresponding crystallographic structures. A data set of 55 protein–ligand complexes with experimental values obtained from PubChem BioAssay (Table 9) were used to perform a correlation with calculated docking affinity values for BPs interacting with ESR1, ESR2, PPARG, PPARD, VDR, and TSHR. The Pearson’s correlation coefficient between in silico docking affinities and the experimental values was with a (Figure 8A).

Figure 8.

(A) Calculated docking affinity scores of bisphenol A (BPA) analogs and derivatives (BPs) with several proteins associated with noncommunicable diseases (NCDs) vs. their experimental half-maximal activity concentration () values obtained from PubChem Bioassay (Wang et al. 2017c). The proteins related to NCDs used for validation purposes were estrogen receptor alpha (ESR1), estrogen receptor beta (ESR2), peroxisome proliferator-activated receptor gamma (PPARG), peroxisome proliferator-activated receptor delta (PPARD), vitamin D3 receptor (VDR), and thyroid-stimulating hormone receptor (TSHR). Superposition of the crystallographic structures and the binding poses resultant from molecular docking of the complexes: (B) bisphenol B/estrogen-related receptor gamma (BPB/ERR3; PDB ID: 1I61) and (C) Bisphenol E/Estrogen-related receptor gamma (BPE/ERR3; PDB ID: 6I64). Crystallographic structures are presented in gray.

Similarly, the structural validation showed that BPs resultant from the molecular docking simulations were located in the correct binding sites and exhibited the same pose reported in the crystallographic structures with ERR3. Both, BPB (PDB ID: 6I61) and BPE (PDB ID: 6I64) obtained RMSD values , and affinity values of and , respectively (Figure 8B,C). Besides, the molecular docking affinity scores of the set of co-crystalized protein–ligand complexes obtained from PDB ranged from to (Table S1). The MD simulation confirmed the stability of the complex that exhibited the best affinity score (the highest in silico affinity absolute value) in the vHTS, BTUM/DHB1, with an average RMSD of of the atomic positions for the dynamics and static models (Figure S10).

Discussion

Fifty BPs with limited or no available toxicological information were identified through a literature search of articles reporting the association of these molecules as analogs, substituents, or derivatives of BPA. This search was conducted using different databases and search engines, such as Science Direct (http://www.sciencedirect.com/), PubMed (https://www.ncbi.nlm.nih.gov/pubmed/), and Google Scholar (https://scholar.google.com/). The generated list of BPs may contribute in the development of further studies to advance in the characterization, evaluation of the potential toxicological effects, and monitoring of these BPs. On the other hand, the identified proteins associated with NCDs aided to evaluate in silico their potential interaction with these emerging contaminants.

The vHTS performed in this research revealed that BPs may have the potential to bind proteins related to NCDs with strong affinity. BTUM, BDP, BPPH, and Pergafast 201 presented the best affinity scores with numerous proteins associated with cardiovascular diseases, lung cancer, breast cancer, cervical cancer, prostate cancer, diabetes, and thyroid disorders (top 10 affinity scores for each disease). According to the U.S. EPA, the compounds BTUM, BPPH (also referred as BisOPP-A), and Pergafast 201 have been used as BPA alternatives in thermal paper (U.S. EPA 2014) and polymers (Zühlke et al. 2020). On the other hand, BDP is used as flame retardant (He et al. 2009; Jing et al. 2013). Interestingly, BPPH has been found in aquatic systems in China (Catenza et al. 2021), and BDP has been detected house dust from Norway and the United Kingdom (Kademoglou et al. 2017). However, more research is needed to determine the occurrence and concentrations of BPs in products for human consumption, as well as in the environment.

According to the World Health Organization, cardiovascular diseases are the leading cause of death globally due to NCDs (WHO 2017). In this study, we found that several BPs present in thermal paper, such as BTUM, Pergafast 201, BPPH, and D-90 (U.S. EPA 2014); the flame retardant BDP (He et al. 2009; Jing et al. 2013); and numerous precursors of polycarbonate plastics and epoxy resins, such as BPFL (Anderson 2020), DGEBP (Zhang and Vyazovkin 2006), Bisphenol M (BPM) (Kricheldorf et al. 2005), 3,5-bis(trifluoromethyl)phenylhydroquinone (BTFMHQ) (Jiang et al. 2018), and Bisphenol P (BPP) (Kricheldorf et al. 2005; Owczarek et al. 2018), exhibited good affinity scores () with proteins involved in cardiovascular diseases. The predicted interactions among these BPs and proteins associated with cardiovascular diseases represent an emerging health concern given that these targets have been reported to mediate the development of cardiovascular pathologies. FGFR3 has been suggested to play a role in cardiac signaling in vitro and in vivo (Touchberry et al. 2013), possibly through its interaction with fibroblast growth factor 23 (FGF23), which has been identified as a stimulator of left ventricular hypertrophy and a marker for cardiovascular risk in human observational studies (Pöss et al. 2013; Stöhr et al. 2018; Udell et al. 2012). CETP has been implied in the regulation of high-density lipoprotein (HDL) and low-density lipoprotein (LDL) cholesterol levels in humans (Brousseau et al. 2004). PPARG mediates the differentiation of adipocytes and insulin signaling (Kintscher and Law 2005; Leonardini et al. 2009; Sauma et al. 2006). In addition, FOLH1 regulates the absorption of folates (Guo et al. 2013; Hiraoka and Kagawa 2017), which have been found to contribute to the prevention of atherosclerotic cardiovascular disease (Verhaar et al. 2002).

The BPs with the highest in silico affinity absolute values for proteins associated with lung cancer were the compounds found in thermal paper: BTUM, Pergafast 201, D-90, and BPPH (U.S. EPA 2014); the flame retardant BDP (He et al. 2009; Jing et al. 2013); and BPFL, employed as a copolymer in the production of polycarbonate plastics (Anderson 2020). These molecules interacted in silico with numerous proteins related to lung cancer, one of the most common cancers worldwide (WHO 2018a). These proteins have been associated with lung cancer progression and metastasis in human observational studies and in in vitro analysis with cancer cells. Among them, LDHA (Koukourakis et al. 2003; Yang et al. 2014), fructose-bisphosphate aldolase A (ALDOA) (Chang et al. 2017), receptor of activated protein C kinase 1 (RACK1) (Shi et al. 2012), Arf-GAP with SH3 domain, ANK repeat, and PH domain-containing protein 3 (ASAP3) (Fan et al. 2014), and EGFR (Bethune et al. 2010). The best affinity score (the highest in silico affinity absolute value), in the group of proteins related to lung cancer, was obtained for Pergafast 201 with LDHA, which was predicted to bind in a different pocket than the reported for pyrazole-based inhibitors (Rai et al. 2017) with the same protein structure (PDB ID: 5W9J). The inhibition of LDHA by small molecules has been considered as an anticancer target in drug design (Feng et al. 2018; Wang et al. 2017b; Xian et al. 2015).

On the other hand, many BPs exhibited estrogenic activity in MCF-7 cells and were suggested to act as EDCs (Rivas et al. 2002). In a previous article, we found that some EDCs, including some BPs, such as BPM, BPB, and BPAF were predicted to bind proteins related to breast cancer in silico (Montes-Grajales et al. 2016). In the present study, the greater number of BPs assessed by vHTS allowed the identification of new targets and compounds to be prioritized for in vitro and in vivo evaluation. Among them, the compounds described as replacements for BPA in thermal paper: BTUM, Pergafast 201, D-90, and BPPH (U.S. EPA 2014), as well as the fire retardant BDP (He et al. 2009; Jing et al. 2013); and some BPA substituents used in the production of polycarbonate plastics and epoxy resins, such as BPFL (Niu et al. 2017), TMBPF (Maffini and Canatsey 2020), BPBP (Česen et al. 2018), DGEBP (Zhang and Vyazovkin 2006), BPP (Wang et al. 2019a), BPAP (Xiao et al. 2018), and BPM (Niu et al. 2017).

These findings suggest that some BPs theoretically bind proteins associated with breast carcinogenesis or hormone imbalance. Among them are DHB1, which has been found to catalyze the synthesis of the hormone estradiol which elicits a role in the growth and proliferation of malignant breast tumors in humans (Mazumdar et al. 2009; Pasqualini et al. 1997); the nuclear receptors retinoic acid receptor-related orphan receptor–alpha (RORA) and ROR-gamma (RORG), which were shown to be associated with breast cancer tumor suppression (Du and Xu 2012) and the regulation of metastatic potential of breast cancer cells (Oh et al. 2016), respectively; as well as breast cancer type 1 susceptibility protein (BRCA1) and serine/threonine-protein kinase Chk2 (CHK2), which were reported to mediate breast cancer tumorigenesis in vivo (McPherson et al. 2004). The interactions between BTUM and the protein DHB1, the complex that obtained the best affinity score in the group of proteins related to breast cancer, showed that the ligand is located in the same binding site reported for an inhibitor of DHB1 that has been proposed as an alternative for cancer therapy (Mazumdar et al. 2009). Some of the shared contact residues were VAL225X, PHE259X, LEU149X, VAL143X, and GLY144X and the hydrogen bond with SER142X.

Several BPs also presented good theoretical affinity () with proteins related to cervical cancer. These included the compounds detected in thermal paper: BTUM, Pergafast 201, BPPH, BPS-MPE (U.S. EPA 2014), and BPZ (Björnsdotter et al. 2017); the flame retardant BDP (He et al. 2009; Jing et al. 2013); and several BPA substituents used in the production of polycarbonate plastics and epoxy resins, such as BPP (Kricheldorf et al. 2005; Owczarek et al. 2018), BPM (Kricheldorf et al. 2005), DGEBP (Zhang and Vyazovkin 2006), and BTFMHQ (Jiang et al. 2018). Some of the protein targets of these molecules were the sirtuins SIR2 and SIR7, related to epigenetic and metabolic regulation in cancer cells (Kalmath et al. 2014; Yu and Guo 2010), as well as EYA2 and SRC involved in cell transformation, migration and metastasis in different cancer types, including cervical cancer (Bierkens et al. 2013; Hou et al. 2013; Krueger et al. 2014; Li et al. 2017).

In addition, BPs used in thermal paper, flame retardants, polycarbonate polymers, and epoxy resins (Badrinarayanan et al. 2008; U.S. EPA 2014; He et al. 2009; Jing et al. 2013; Zühlke et al. 2020) were found to bind proteins related to prostate cancer with strong affinity (top 10 affinity scores for each disease). Among them are BTUM, BPPH, Pergafast 201, BDP, BPP (Kricheldorf et al. 2005; Owczarek et al. 2018), BPM (Kricheldorf et al. 2005), BPFL (Anderson 2020), DGEBP (Zhang and Vyazovkin 2006), and BPTMC (Chenet et al. 2020). The predicted protein targets of these BPs included the DNA repair protein PARP1 (Ko and Ren 2012; Schiewer and Knudsen 2014); SRPK1, which was reported to exhibit a higher expression in samples of human prostate tumor tissue compared with benign tissue (Bullock et al. 2016); flavin-containing amine oxidase domain-containing protein 2 (KDM1A), involved in the epigenetic regulation of diverse cancers in humans (Ismail et al. 2018; Zhang et al. 2019), including prostate tumorigenesis (Ismail et al. 2018); peroxiredoxin-1 (PRDX1), which has been reported to promote survival and growth of prostate cancer cells in vitro (Dasari et al. 2019); the sister chromatid cohesion protein PDS5 homolog B (PDS5B), which is a target of miR-27a and which was able to enhance androgen-stimulated pancreatic cell viability (Takayama et al. 2017); and the prostatic acid phosphatase (PPAP). Low levels of the latter protein were associated with poor prognosis of human prostate cancer by targeted proteomics in urinary samples (Sequeiros et al. 2017).

The complex with the best affinity score (the highest in silico affinity absolute value) in the group of proteins associated with prostate cancer was obtained for BPPH with PARP1. This BPA substituent used in thermal paper occupied the same binding site reported for PARP1 inhibitors. The shared contact residues of these interactions were TYR907B, ILE872B, TYR896B, ALA880B, and ALA989B and a hydrogen bond with GLY863B (Salmas et al. 2016). Therefore, structural moieties of BPPH could be useful for the development of new and safe inhibitors of this protein.

BTUM, Pergafast 201, BPPH, BPS-MPE, D-90, BPP, BDP, BPM, and BPFL, among others, exhibited good theoretical affinities () for proteins associated with diabetes. These compounds are used in the production of thermal paper, flame retardants, polycarbonate plastics, and epoxy resins (Anderson 2020; Björnsdotter et al. 2017; U.S. EPA 2014; He et al. 2009; Jing et al. 2013; Kricheldorf et al. 2005; Niu et al. 2017; Zühlke et al. 2020). According to the vHTS, these BPs were predicted to target several proteins that have been associated with the regulation of glucose metabolism, pancreatic exocrine function, or energy production by human observational studies, and in vitro experiments. Some of them included CEL (Ræder et al. 2014), aldo-keto reductase family 1 member B1 (ALDR; also referred as AKR1B1) (Donaghue et al. 2005), PPARG (Bermúdez et al. 2010; Wang et al. 2019b), fructose-1,6-bisphosphatase 1 (F16P1; also known as FBP1) (Wang et al. 2019b), as well as the cytochrome b5 reductase 4 (NB5R4), which has been reported to mediate the development of hyperglycemia in mice (Wang et al. 2011).

The BPs that presented the highest in silico affinity absolute values for proteins associated with thyroid disorders were BTUM, Pergafast 201, BPPH, and BPS-MPE, which are used in thermal paper (Björnsdotter et al. 2017; U.S. EPA 2014); the flame retardant BDP (He et al. 2009; Jing et al. 2013); as well as BPM and DGEBP, which have been used in the production of polycarbonate plastics and epoxy resins (Kricheldorf et al. 2005; Zhang and Vyazovkin 2006). The protein targets of these BPs were mainly associated with thyroid cancers. Among these are SIR6, which has been found to enhance cell aggressiveness in human papillary thyroid cancer (Qu et al. 2017); KEAP1, which has been related to critical thyroid cancer in humans through a mechanism predominantly mediated by its hypermethylation and subsequent inactivation (Martinez et al. 2013); MK03 (also known as MAPK3 or ERK), which has been reported to participate in the signaling of tumorigenesis in humans (Kohno and Pouyssegur 2006), and MERL, which has been described as a negative regulator of cancer cell growth and proliferation in vitro and in vivo (Stamenkovic and Yu 2010). Other targets include thyroid hormone receptors, such as THA and THB, which regulate the balance of thyroid hormones in humans (Bochukova et al. 2012). Furthermore, the study of the interactions of the complex with the best affinity score (the highest in silico affinity absolute value), in the group of proteins related to thyroid disorders, revealed that BTUM is predicted to bind SIR6 in a slightly different pocket than the described as the binding site for molecular activators of this protein, by sharing only two contact residues PHE64A and TRP188A (You et al. 2017).

Several targets of BPs identified by vHTS, belonging to the protein–protein interaction network generated by STRING 11.0 (Szklarczyk et al. 2019) (Figure S9), were predicted to participate in pathways associated to different processes. NAD(P)H dehydrogenase [quinone] 1 (NQO1), FGFR3, EGFR, SRC, PPARG, BRCA1, and KEAP1 have been related to cancer signaling (KEGG pathway IDs: hsa05200, hsa05219, hsa05206, and hsa05225). ALDOA, LDHA, EGFR, and FGFR3 have been found to participate in central carbon metabolism, gluconeogenesis, and glycolysis (KEGG pathway IDs: hsa05230, hsa00010). RORA and RORG were linked to circadian rhythm and Th17 cell differentiation (KEGG pathway IDs: hsa04710, and hsa04659). The proteins EGFR, FGFR3, and SRC belong to the ErbB, PI3K-Akt and endocytosis signaling pathways (KEGG pathway IDs: hsa04012, hsa04151, and hsa04144), which were related to cell proliferation, differentiation and survival, as well as cell cycle regulation and signaling. Furthermore, these proteins were predicted to participate in other pathways related to EGFR tyrosine kinase inhibitor resistance, regulation of actin cytoskeleton, adherens junctions, gap junctions, endocrine resistance, as well as GnRH and RAP1 signaling (KEGG IDs: hsa01521, hsa04810, hsa04520, hsa04540, hsa01522, hsa04912, and hsa04015), among others (Kanehisa et al. 2019; Ogata et al. 1999).

Proteins highly interrelated in the protein–protein interaction network, such as EGFR, PPARG, SRC, and PPAR1, were also associated with key pathways such as cancer signaling, cell proliferation, and carbon metabolism by STRING 11.0. (Szklarczyk et al. 2019). Therefore, the mechanism of action of BPs with strong theoretical affinities for these proteins could be mediated by protein–ligand interactions with effects in more than one health condition. Some of these compounds were BTUM, BDP, BPP, BPM, BPFL, Pergafast 201, and BPPH. However, due to the limitations of our protocol and the lack of toxicological information, further in vitro and in vivo studies are needed to gain a better understanding of the effects and mechanisms of action of these emerging pollutants.

According to the validation process, the protocol used for the vHTS showed a good predictability in terms of the correct identification of the binding site and ligand pose compared with the corresponding crystallographic structures. As well, there was a good correlation between calculated binding affinities and experimental values (). The values were obtained from quantitative high-throughput screening experiments for the identification of agonistic or antagonistic activity with the protein targets. However, these came from different experiments and should be interpreted carefully given that they may contain a considerable variability due to changes in experimental conditions.

Furthermore, numerous protein–ligand complexes formed by BPs with proteins involved in NCDs (identified through the vHTS) presented better affinity values than the protein–ligand complexes with crystallographic structures ( and ) tested as part of the validation process. This suggests that the proposed BP–protein complexes with the highest in silico affinity absolute values are promising candidates to be prioritized for in vitro and in vivo testing. In addition, the MD simulation showed a good stability of the protein–ligand complex with the strongest affinity score (the highest in silico affinity absolute value) resultant from the vHTS.

Conclusions

This study reveals that BPs have the potential to target proteins associated with NCDs of high prevalence worldwide, as well as their related pathways. Therefore, BPs may be eliciting their effects in NCDs by interacting with proteins involved in crucial pathways of cancer, cell cycle regulation, signaling, and metabolism, among others. Some of them have better in silico binding affinity than BPA. This represents an emerging public health concern because the toxicological information about them is very limited, and humans can be highly exposed to them through everyday products such as thermal paper receipts, electronic devices, or plastic containers. Computer-aided approaches like those used in the present study contribute to improving the speed of the assessment of emerging pollutants, such as BPs. However, the leading compounds should be further examined in vitro and in vivo to confirm their possible interaction, and to elucidate their mechanisms of action. Therefore, BPs that presented strong in silico affinity for proteins involved in NCDs, such as BTUM, BDP, BPPH, and Pergafast 201, are proposed as high priority compounds to be further assessed by using in vitro and in vivo models of these diseases.

Click here for additional data file.

Click here for additional data file.

Click here for additional data file.