Comprehensive Analysis and Large-Scale Screening of Binding Interactions Between PFAS and Their Mixtures with Nuclear Receptors

1. Introduction

Per- and polyfluoroalkyl substances (PFAS), synthetic for-ever-chemicals used for 60 years, are exposed to 99% of the human population every day. PFAS are predominantly encountered in the environment, mainly through food and water consumption. They are integral to various consumer goods, providing waterproof, grease proof, and nonstick properties (Pelch et al., 2019). Drinking water, in particular, has been identified as a major source of PFAS exposure to humans (Domingo et al., 2019). Beyond their widespread use, PFAS are increasingly recognized for their adverse health effects, including cancer, metabolic issues, and reproductive issues. Higher PFAS exposure has been linked to irregularities in menstrual cycles and increased cycle length (Ding et al., 2020). While PFAS offer functional benefits in products, their persistence and bioaccumulation raise considerable health concerns, including cancer risk. The carcinogenic potential of PFAS is attributed to factors such as action-metabolism disruption, endocrine system interference, and epigenetic changes (Boyd et al., 2022). Studies have shown associations between hepatocellular carcinoma (HCC) and PFAS, where PFOS was regulated and the study yielded when exposed to greater amounts of PFOS, increased signs of HCC occurred (Goodrich et al, 2022). Outside the scope of cancer, PFAS continues to harm humans through the impact on reproductive systems, developmental stages, and other organs in the body leading to kidney or liver disease (Fenton et al, 2021).

More alarmingly, there is growing evidence of PFAS pro-carcinogenic effects, especially concerning kidney, breast, ovarian, liver, prostate and testicular cancers. In the context of liver cancer, research has revealed a correlation between hepatocellular carcinoma (HCC) and PFAS exposure. Shearer et al showed the correlations between increased concentrations of PFOA and other PFAS types and renal cell carcinoma (RCC) (Shearer et al 2021). Goodrich et al. further reported studies regarding PFOS exposure and the causation of HCC (Goodrich et al. 2022). Associations between PFAS exposure and breast cancer are becoming prevalent as PFAS interfere with mammary gland development and hormonal functions, potentially elevating breast cancer risk (Lopez-Espinosa et al., 2021). The endocrine-disrupting potential is highest for PFAS, as they have been shown to increase risk of breast cancer at higher concentrations (Barry et al. 2023). Jones et al. (2022) reported that individuals with higher PFAS exposure showed a greater incidence of prostate cancer, due to PFAS-induced hormonal imbalances. Importantly, there is evidence of associations between PFAS exposure and kidney, testicular, ovarian, thyroid, endometrial cancers (Shearer et al 2021) (Rhee Chang et al 2023) (Purdue et al 2023) (Chang et al 2023) (Rhee Barry et al 2023) (Madrigal et al 2023) (Jones et al 2023).

The physicochemical properties of PFAS significantly influence their binding interactions with nuclear receptors. The hydrophobicity and stability of PFAS contribute to their persistence and potential bioaccumulation, further influencing their ability to interact with and disrupt the normal functioning of nuclear receptors. Such interactions can lead to significant biological effects, given the role of nuclear receptors in regulating gene expression related to development, metabolism, homeostasis and cancer. Understanding these interactions is essential for assessing the toxicological impacts of PFAS on human health. PFAS consist of a carbon backbone attached to fluorines, which results in a nonpolar hydrophobic nature, similar to that of the endogenous ligands of nuclear receptors like VDR, PPARγ, PXR, and ERα. This nonpolarity allows the chemicals to bind to nonpolar ligand binding domains, composed of mainly nonpolar amino acids. Long and short PFAS chains have been shown to disrupt the function of nuclear receptors such as PXR causing the overactivation of the receptor, and leading to endocrine disruption, oxidative stress, and hepatic steatosis (Lai et al 2020). It has been shown that PFOA competitively binds at the VDR (Singam et al 2023). Beyond their carcinogenic implications, PFAS continue to pose significant health risks, impacting reproductive health, developmental stages, and causing damage to vital organs, leading to conditions like kidney or liver disease (Fenton et al., 2021). For example, when PFAS are bound to the VDR, the VDR cannot properly regulate genes that are involved in calcium and phosphate homeostasis, immune response, and cellular proliferation and differentiation (Wang et al, 2012). PFOA, when bound to PPARγ, has shown anomalies in bone development and adipose tissue homeostasis (Kirk et al 2021).

Previous studies showed correlations with greater PFAS exposure and increased rate of thyroid cancer (van Gerwen et al, 2023). Colorectal cancer susceptibility increases when the VDR is dysregulated, a condition that can be induced by PFAS through their binding to the VDR, as well as to other nuclear receptors (Issa et al 2022). The VDR has been shown to regulate blood pressure, work against organ failure, and play as a factor for autoimmunity, angiogenesis, inflammation and vascular cell activity (Lin et al, 2019). Given the regulatory functions of the VDR, when the competitive binders of calcitriol like PFOA bind instead, structural flexibility is decreased leading to an altered response of the vitamin D responsive genes (Di Nisio et al, 2020). In the VDR, 1,25-dihydroxyvitamin D (1,25(OH)2 D) is an agent that works against age-related osteoporosis (Yang et al, 2019). As more of the PFAS are able to bind the VDR, the function of 1,25-dihydroxyvitamin D is limited and that leads toward a higher likelihood of developing osteoporosis. There are many common PFAS molecules that could trigger the effects of osteoporosis, like PFOA or PFOS as these are the most abundant in relation to humans (Liang et al, 2022).

The PXR is involved in the regulation of intermediate metabolism through trans-activation and trans-repression of genes controlling glucose, lipid, cholesterol, bile acid, and bilirubin homeostasis (Pavek 2016). Like any other member of the NR superfamily, PXR is composed of a DNA binding domain (DBD), involved in receptor dimerization and binding of specific DNA sequences, an H region, a flexible domain connecting the DBD with the LBD, and a C-terminal LBD (Pavek 2016). PXR was also originally characterized as the key transcription factor that activates hepatic genes, econding drug-metabolizing enzymes and drug efflux transporters, and protects the body from harmful foreign toxicants or endogenous toxic substances (Pavek 2016). Studies have been conducted, detailing the interactions between PFAS and the human PXR receptor, with binding typically leading to an overactivation of the receptor (Lai et al 2020). Overactivation of PXR is linked with potential endocrine disruption, oxidative stress, hepatic steatosis, and other adverse drug interactions (Lai et al 2020).

The PPARγ plays a major role in cell differentiation, lipid metabolism, and glucose homeostasis (Villapol 2019). Thiazolidinedione (TZD) are agonists that bind to the PPARγ, which are prescribed for the treatment of type II diabetes. When associated with PPARγ, PFAS have been shown to have effects on adipose tissue and bone differentiation, which leads to significant cancerous effects (Kirk et al 2021). Some adipose tissue related effects like obesity have been associated with younger children that had a higher exposure to specific PFAS like PFOA and PFHxS (Kirk et al 2021). Inclusive to young children, prenatal babies also experience negative effects when exposed to PFAS. Correlations have shown with increased PFAS concentration, there is an increased count of adipose tissue, leading to a greater BMI and fat deposition.

The Estrogen receptors (ERs), are ligand-activated nuclear hormone receptors (NHRs), responsible for gene expression in the nucleus. The native ligand that activates ER𝛼 is estrogen, the primary female sex hormone. There are four types of estrogen including estrone, estradiol, estriol, and estetrol, however, estradiol is the most dominant and prevalent during primary and secondary reproductive years, making it the reference ligand of choice in our study. As a steroid, estradiol diffuses into the cell and binds to ER𝛼 in the cytoplasm. ER𝛼 is structurally homologous to most nuclear receptors, composed of six domains, A/B, C, D, and E/F, for 595 amino acids. PFAS have shown to bind to human ER𝛼 at greater rates, resulting in an increased activation of the estrogen receptor and resulting in overexpression of estrogen (Qiu et al 2020). Additionally, there are two regulators of transcriptional activity, activation factor-1 (AF-1) and activation factor-2 (AF-2) which are hormonally activated. The ligand binding domain (LBD) is located within domain E, the largest domain of the protein. In activation, estradiol (E2) binds to the hydrophobic pocket. When E2 binds to the pocket, the hydroxyl on the A ring forms a hydrogen bond to Glu353 on helix 3 and Arg394 from helix 5, while the terminal hydroxyl of E2 binds with His524 on helix 11 (Bafna et al 2020). Estrogen bound ER𝛼 undergoes a conformational change moving alpha helix 12 to an open position. The movement of H12 opens up the Activation Factor 2 (AF-2) site for subsequent coactivator recruitment and leads to ER𝛼 dimerization. In this study, we exploited orthosteric and allosteric sites containing ER𝛼 LBD and AF-2 sites for PFAS docking simulations.

Given the hazardous effects and ability to disrupt NR functionality, in this study we performed PFAS binding interactions with NRs at the orthosteric and allosteric sites. We report the results of the specific nuclear receptors: VDR, PPARγ, PXR, and ER𝛼. Distinguishing itself through a comprehensive screening process, this study evaluates 9,507 PFAS molecules against these nuclear receptors. By employing docking analysis, we provided a comprehensive understanding of the molecular interactions between PFAS in a standalone and mixed context at the orthosteric and allosteric sites of the nuclear receptors. Our findings demonstrated a large range of PFAS that had an ability to bind to the VDR, PPARγ, PXR, and ERα, with binding affinities comparable to or greater than those of endogenous ligands, at orthosteric and allosteric sites.

2. Methods and Materials

2.1. Overall Methodology Workflow

Our docking process analyzes the binding interactions between receptor and ligands to understand how different PFAS chemicals can compete with the nuclear receptors at their orthosteric and allosteric sites that cause dysregulation to nuclear receptors. We first retrieved all the required protein files from the rcsb.org protein data bank; two files per protein, one for allosteric and the other for orthosteric sites. Along with the protein files, we obtained the endogenous ligand files from the same protein data bank website. Next, we retrieved all the PFAS files from the Environmental Protection Agency (EPA) CompTox chemical dashboard (https://www.epa.gov/comptox-tools/comptox-chemicals-dashboard). We ran Autodock Vina on the endogenous ligands with the corresponding protein files, and placed both files in UCSF ChimeraX to visualize the superimposition. Once Autodock vina proved to have an accurate docking for the endogenous ligands, we performed the same docking procedure for the PFAS chemical screening, using the python code (Supplementary file Python Code). The output is a large array of PFAS, ranked by their binding affinity in kcal/mol, with their proper conformation with respect to the specific nuclear receptors’ binding sites. Once we had the ranked PFAS at the respective binding sites for NR proteins, we analyzed the high energy binding PFAS structures to understand which PFAS will have higher likelihood of competing with the NRs activity (Supplementary Figures - Suppl 1). Finally, we developed several PFAS mixed model scenarios where we displayed combinations of PFAS and endogenous ligands that can bind and interfere NRs activity and induce dysregulatory effects.

2.2. Receptor and Native Ligand Preparation and Docking Simulations

We retrieved orthosteric and allosteric structures from the Protein Data Bank: one containing the Vitamin D Receptor and 1,25-dihydroxyvitamin D3 (calcitriol) (PDB:1IE9), and the other consisting of the Vitamin D Receptor and (3α,5β)-3-hydroxy-cholan-24-oic acid (lithocholic acid - LCA) (PDB:4Q0A). The Pregnane X Receptor (PXR) was acquired, one for orthosteric site docking (PDB:5X0R) and other for allosteric AF-2 site docking (PDB:7AXB). In the case of PPARγ, two individual structures of the ligand-binding domain (LBD) of PPARγ with known synthetic ligands were obtained (PDB: 3ET3, and 5GTO). Docking analysis was conducted for the ERα at the orthosteric (PDB: 1ERE) and allosteric sites (PDB: 3UUD). Using ChimeraX (Pettersen et al., 2021; Goddard et al., 2018), all water molecules and ligands were removed from the receptor. Subsequently, these receptor files were processed with Autodock 4/4.2 Tools (Morris et al 2009), adding polar hydrogens to facilitate intermolecular bonding with ligands, merging of nonpolar hydrogens, the assignment of Kollman charges, and setting the pH to 7. The Kollman charges were selected; they are calculated using the AMBER force field and are based on a molecule's molecular orbital electron density, which makes them more suitable for accurately representing proteins and other large biological molecules. Finally, we assigned the AD4 type forcefield. Subsequently, the protein was exported in .PDBQT format, for Autodock Vina.

All ligands were reprocessed using a software Gypsum DL (Ropp et al, 2019). In Gypsum dl, we input smiles (.smi) version of the ligands, and 3D .sdf output files of the ligands were obtained. We created a python code file to convert the 3D .sdf files into 3D .pdbqt files that was optimized for the best conformer structure (Supplementary file Python Code). In preparation for docking, we used Autodock tools to locate the allosteric and orthosteric sites. Next, we created a grid box with specific dimensions, for the VDR the center was 10.839, 21.703, 33.959 (x, y, z), with a size of 38, 36, 18 (x, y, z) grid box units at the orthosteric site. For the VDR allosteric site the docking was positioned on the center at 6.395, 46.188, 25.575 (x, y, z) with dimensions of 16, 40, 18 (x, y, z) grid box units. For the PPARγ orthosteric binding site, the grid box was centered on the amino acid TYR327 and encompassed other known amino acid residues present in the active site such as HIS323, TYR472, and HIS449. For the PPARγ allosteric binding site, the grid box was centered on the amino acid LEU270 and encompassed other known amino acid residues present in the allosteric site, and the grid box size was 30Å in each dimension (x, y, z). For the PXR orthosteric site, we defined the grid box parameters with 56, 60, 52 (x, y, z), with center coordinates of 36.163, 16.732, 26.889 (x, y, z). These specific parameters were chosen with the aim of encompassing the Ser247, Gln285, His407, and Arg410 polar residues, which play a crucial role in defining the ligand binding position in PXR. For the ER𝛼 orthosteric site we chose the coordinates of residue ARG-394 (45.437, 55.270, 44.046). To cover the size of the ligand binding, the grid box size was set to 26, 20, 20 (x, y, z). Following this, the parameters were exported as a .GPF file, and Autogrid4 was executed to generate a .GLG file for the respective docking procedure. These dimensions were then shifted from Autodock Tools to a .txt file where the configuration was kept so the docking would be consistent to the specific binding domain.

The .PDBQT file of the protein was set as a macromolecule when preparing the orthosteric and allosteric docking file. Subsequently, we selected the .PDBQT ligand file as the ligand and ensured the accuracy of the atom types. For the docking process, a genetic algorithm with a long run containing a maximum of 25,000,000 evaluations, and 20 ligand conformations was employed. These parameters were saved as a .DPF file, utilizing the Lamarckian GA algorithm due to its improved efficiency and enhanced performance compared to the genetic algorithm alone. Finally, we executed Autodock4 to generate a .DLG file, which was used for further docking analysis. The docking was restricted within the orthosteric and allosteric sites of the VDR, PXR, PPARγ, ER𝛼 to obtain a greater success rate (Mukherjee, Balius, Rizzo 2010). To identify the accuracy of the docking,and for the controlled simulations, we redocked the endogenous ligands in their respective binding domains. For the VDR we redocked the calcitriol and lithocholic acid molecules onto the orthosteric and allosteric site respectively, using Autodock Vina (Eberhardt et al, 2021; Trott et al, 2010). Similarly, we carried out redocking of endogenous ligands for PPARγ, PXR and ER𝛼.

2.3. PFAS Preparation and Docking

The chemical structures for the PFAS were downloaded from the Environmental Protection Agency (EPA) CompTox chemical dashboard in June 2023 (https://www.epa.gov/comptox-tools/ comptox-chemicals-dashboard). Starting with 12,034 chemicals, we filtered out the PFAS that had no structures, multiple structures, and isotopes, leaving us with 9,903 chemicals. We used the DTX SID as a way of differentiating the PFAS chemicals, as it is the chemical identifying method for CompTox. The chemical structures of PFAS were downloaded as SDF files. Using OpenBabel 2.3.2, the SDF file was converted to a SMILES file containing all the PFAS chemicals with SMILES codes. The PFAS chemicals without proper SMILES codes were discarded and removed from the dataset. Total of 9,870 PFAS chemicals were prepared using Gypsum-DL by converting the 2D PFAS structures into accurate 3D compounds. In addition to generating the 3D structures, the compounds were neutralized at a pH of 7. The 9,870 compounds successfully prepared using Gypsum-DL were energy minimized using obminimize from OpenBabel 2.3.3, and using the general amber force field (GAFF) with default parameters were converted into PDBQT files. We developed a Python code (Supplementary file Python Code) capable of optimizing the conformations of all ligands using the AMBER force field and the Obminimize software. This optimization process resulted in the generation of .PDBQT output files for each ligand to be used in the Autodock Vina software, and in turn sped up the processing of the PFAS files.

2.4. Molecular Docking

Molecular docking simulations were conducted using Autodock Vina employing the AMBER force field. The simulations involved interactions between proteins, native ligands, and PFAS. Specifically, the docking simulations were focused on the orthosteric and the allosteric sites with the corresponding endogenous ligand. Grid box coordinates are created as described above in section 2.2. Autodock Vina based virtual screening was executed for 9,870 PFAS molecules onto the orthosteric site. Following completion of the orthosteric site, the grid box parameters were adjusted to target the allosteric site where the endogenous ligand binds, and all PFAS molecules were subsequently docked onto this binding site. The docking protocol and the accuracy of Autodock Vina predictions were verified using cross docking of native ligands at both the orthosteric and allosteric sites for the VDR, PXR, ERα, and PPARγ.

3. Results and Discussion

3.1. Docking Verification: Cross-Docking of Known Ligands and Non-Native Ligands and Comparison with Experimental Binding Affinity

We examined the interactions of PFAS, native ligands, and non-native ligands to the VDR, PXR, PPARγ, and ERα. Our approach involved validating Autodock Vina docking protocol by docking known ligands and confirming the accuracy of native ligands by cross-docking of non-native ligands for each nuclear receptor. These steps are important to ensure the reliability of our docking study protocol. Our docking results demonstrate that native ligands are docked correctly to their respective nuclear receptors, as seen by low RMSD values from the X-ray binding ligand conformation, highlighting the precision of our Autodock Vina protocol. Our examination of non-native ligands docking results further validates our Autodock Vina docking protocol.

3.1.1. Cross-Docking of Known Ligands: VDR

We started with verifying the accuracy of Autodock Vina predictions using cross docking of native ligands at the orthosteric and allosteric sites of the VDR (Figure 1). We docked two native ligands to verify the molecular docking procedure. The docking of the native ligands, calcitriol (PDB ID: 1IE9, resolution: 1.40 Å) and lithocholic acid (PDB ID: 4Q0A, resolution: 1.90 Å), were compared to the co-crystalized structures. The superimposed structures for calcitriol and lithocholic acid are shown in Figure 1. The docking at both the orthosteric and allosteric sites yielded similar accuracies, with calcitriol displaying a RMSD of 1.31 Å and lithocholic acid exhibiting an RMSD of 1.306 Å. The amino acid residues within 5 Å of the orthosteric site include TYR 143, TYR 147, PHE 150, SER 237, LEU 233, ARG 274, CYS 288, GLU 277, SER 278, SER 275, TYR 295, ILE 271, TRP 286, ASN 276, PHE 279, LEU 313, MET 272, VAL 300, LEU 309, HIS 305, HIS 397, ILE 268, ALA 303, LEU 404, TYR 401, LEU 230, VAL 234, ALA 231, LEU 227, VAL 418, and LEU 414. In the co-crystalline structure, calcitriol forms hydrogen bonds with HIS 397 and TYR 143, and in the redocked structure, it forms additional hydrogen bonds with ARG 274, SER 278, and HIS 305. The amino acid residues surrounding the allosteric site within 5 Å are ARG 184, ASP 181, ASP 260, TYR 264, SER 263, LYS 268, LEU 443, PHE 185, THR 441, VAL 444, and GLN 267. Although the co-crystalline structure of lithocholic acid did not show hydrogen bonding with the VDR, the redocked ligand formed a hydrogen bond with SER 263, demonstrating the flexibility of the structure and its capacity to optimize binding with VDR. The docking score for calcitriol is -12.3 kcal/mol, indicating strong interactions with the VDR LBD, whereas the docking score for lithocholic acid is -6.6 kcal/mol, suggesting much weaker interactions.

3.1.2. Cross-Docking of Known Ligands: PPARγ

To confirm the accuracy of the docking protocol for the PPARγ orthosteric binding site, the binding poses of native ligand, ET1, was analyzed. The pose obtained via molecular docking was compared with the co-crystallized structure of ET1 (PDB ID: 3ET3), resulting in a RMSD of 1.017 Å. Further, the docked ET1 at the orthosteric site formed a hydrogen bond with T473.The docking protocol was similarly validated for the allosteric binding site. The pose of PPARγ allosteric site native ligand, T35, was compared with the co-crystallized structure of T35 (PDB ID: 5GTO), yielding an RMSD of 0.758 Å. Additionally, the docked T35 formed hydrogen bonds with residues L270 and Q282. These findings validate the molecular docking protocol efficacy in accurately placing ligands within their respective binding pockets in the correct orientation, confirming that the parameters for docking small molecules to the LBD of PPARγ successfully reproduce the experimentally observed binding poses of known ligands.

3.1.3. Cross-Docking of Known Ligands: PXR

To confirm the docking protocol for the orthosteric site of PXR, the native ligand 4WH was re-docked onto the orthosteric site of the PXR. We used this approach to validate the accuracy of the docking protocol by comparing the computationally predicted binding mode to the experimentally determined position and orientation of the ligand in the PXR. Figure 2 a) displayed the superimposed conformation of the re-docked ligand and the co-crystallized ligand on the PXR. The RMSD value was calculated for the superimposition, resulting in a value of 0.793 Å. This value falls within the RMSD criterion of less than 2 Å, indicating a high level of docking accuracy. Our protocol also placed the re-docked reference ligand 4WH in the correct orientation, further confirming the accuracy of our docking protocol for the PXR. To verify the docking protocol for the allosteric AF-2 Site of PXR, we re-docked the glycerol ligand onto the allosteric site of the PXR. Figure 2 b) shows the superimposed conformation of the re-docked ligand and the co-crystallized ligand on the PXR. The RMSD value was calculated for the superimposition, resulting in a value of 0.098 Å. This value falls within the RMSD criterion of less than 2 Å, indicating a high level of docking accuracy. Our protocol also placed the re-docked reference ligand glycerol in the correct orientation, further confirming the accuracy of our docking protocol for the PXR.

3.1.4. Cross-Docking of Known Ligands: ER𝛼

To validate the docking protocol, we redocked endogenous ligands at the orthosteric and allosteric site of the ER𝛼, to produce superimposed conformations. To begin with, at the ER𝛼 LBD, we focused on the endogenous ligand estradiol. We superimposed estradiol onto the 1ERE protein-ligand complex, as depicted in Figure 3. The analysis revealed a docking score of -11.1 kcal/mol for estradiol. Further, the arrangement of the docked ligand closely superimposed the native ligand of the 1ERE complex structure, providing strong validation for our docking methodology. At the AF-2 site, we conducted the same protocol for the steroid receptor coactivator-1 (SRC-1) to show accuracy for docking at the allosteric site. Through computational docking, we noticed the binding of SRC-1 to specific residues, particularly SER-433, within the AF-2 site. The calculations yielded a binding affinity of -9.9 kcal/mol for this interaction, displaying an accurate docking protocol, and the redocked ligand is displayed in Figure 3. The protocol for docking is deemed to be accurate for both protein files for the ER𝛼 due to the redocked ligands superimposing the native ligands from the protein files.

3.1.5. Cross Docking of Non-native Ligands: VDR

To test the accuracy of Autodock Vina’s results, we have also docked other proteins’ ligands, like Chenodeoxycholic, Estradiol, Tretinoin, 15-deoxy-Δ-12,14-Prostaglandin, Triiodothyronine, into the orthosteric and allosteric site of the VDR. The results supported Autodock Vina’s accuracy as it held calcitriol as the greatest docking score among other native ligands for other proteins in their respective orthosteric site, as shown in Table 1. At the orthosteric site, calcitriol has a docking score of -12.3 kcal/mol whereas the highest docking score for the crossed docked ligands is -9.4 kcal/mol by chenodeoxycholic and estradiol. As depicted in Table 1, Autodock Vina’s accuracy has been further established as the same molecules were docked in the allosteric site of the VDR, and still lithocholic acid had the greatest docking score. In comparison to the crossed docked ligands, LCA doesn’t have a far greater docking score as it is -6.6 kcal/mol and estradiol has a docking score of -6.5 kcal/mol. Calcitriol yielded to be far greater binder to the VDR at the orthosteric site than lithocholic acid is to the allosteric site. Resulting in more potential competition at the allosteric site than orthosteric site, and further proving a precise docking protocol.

3.1.6. Cross Docking of Nonnative Ligands: PPARγ

In addition to verifying the protocol's binding pose accuracy, the docking scores' accuracy was also verified. This was done by docking the native ligands of other nuclear receptors to the orthosteric PPARγ binding site and comparing the obtained docking scores to known native and synthetic ligands of PPARγ. In total, one native and one synthetic ligand of PPARγ, as well as 6 non-endogenous ligands, were docked into the LBD of PPARγ (PDB ID: 3ET3) (Supplementary Tables - Suppl 1). The protocol gave the two highest docking scores to the known ligands while the 6 ligands not endogenous to PPARγ were given lower docking scores. These results thus confirm that the parameters for obtaining the docking scores are accurate and can be relayed to the docking of PFAS as depicted in Table 2. Docking scores for known PPARγ ligands as well as known ligands for other nuclear receptors. All ligands are bound to the orthosteric binding site and are within 5Å of known amino acid residues including HIS323, TYR327, and LEU330.

3.1.7. Cross Docking of Nonnative Ligands: PXR

We performed docking experiments with various native ligands of PXR onto the Orthosteric and Allosteric AF-2 Site to evaluate the accuracy of our docking protocol. To reinforce our findings, we compared the binding affinities of these native ligands with those of external receptor ligands (Supplementary Tables - Suppl 2). Vitamin K2 displayed the highest binding affinity of -9.7 kcal/mol among PXR native ligands at the Orthosteric LBD site and Lovastatin displayed the highest binding affinity of -5.6 kcal/mol at the Allosteric AF-2 site. Importantly, all of PXR's native ligands exhibited superior binding affinities to the orthosteric Site compared to the native ligands collected from external receptors, with the exception of PXR’s native ligand Schisandrin A. However, literature studies have shown that this compound activates human PXR with a similar efficacy and potency as rifampicin, which is known to be a relatively weak PXR binder (Ngan et al 2009) (Chang 2009). Similarly, at the AF-2 site, most of PXR’s native ligands displayed higher binding affinities to the allosteric binding site compared to the native ligands collected from external receptors, with the exception of 5-beta-pregnane-3-20-dione and 3-keto-lithocholic acid. However, these ligands had binding affinity values of -4.6 kcal/mol, relatively close to the native ligand with the highest binding affinity that belonged to an external receptor, which was -4.7 kcal/mol. This minor difference could stem from the larger size and intricacy of the former ligands, making their precise docking at the allosteric site more challenging. Given the remarkably close resemblance of these outlier values to the optimal binding affinity of an external receptor’s native ligand, our docking methodology evidently retains its accuracy.

3.1.8. Verification of Experimental Binding Affinity Values for Docking Accuracy

To further validate the accuracy of our docking protocol for the orthosteric site, we collected and analyzed the IC₅₀, EC₅₀, Ki, and Kd values of the native ligands. These experimental binding affinity values will provide empirical support and confirmation of our computational findings. The 4WH reference ligand’s EC₅₀ value of 880 nM was chosen as the primary experimental data point because it surpassed the IC₅₀ value of 340 nM for the reference ligand. Further, this EC₅₀ binding affinity was translated into a value of -8.3 kcal/mol, which closely matched the binding affinity of the re-docked 4WH ligand in its most similar conformation to the co-crystallized ligand within the PXR protein. For the VDR, the Kd value for calcitriol is 0.06 nM, while the Kd value for lithocholic acid is 330 nM. At the PPARγ, T35 had an EC₅₀ value of 183nM. In the ER𝛼, Estradiol had max Ki value of 100nM, max Kd value of 100nM, max IC₅₀ value of 46 nM, and max EC₅₀ value of 10 nM. We were unable to identify the experimental values of the Glycerol, ET1, and SRC-1 reference ligands to validate the accuracy of our docking protocol. Glycerol plays a function when binding to the allosteric site of the PXR, while ET1 binds to the orthosteric site of the PPARγ, and SRC-1 interacts with the allosteric site of the ER𝛼. However, we were able to confirm that the Ile255, Ile273, Leu276, Leu424, Glu427, and Leu428 amino acid residues, which are well documented in the literature as a part of the Allosteric AF-2 binding site for the PXR, were all within a 5 Å distance of the Glycerol ligand, providing strong evidence that our docking protocol for the allosteric site was accurate (Kamaraj et al 2022).

3.2. Virtual Screening of PFAS Molecules against VDR, PXR, PPARγ, and ER𝛼

We processed a total of 9,903 PFAS to be tested onto the VDR, PXR, PPARγ, and ER𝛼 at their orthosteric and allosteric sites. Of the 9,903 chemicals, we excluded 396 ligands from the study as Autodock Vina wasn’t able to dock them onto the protein. 9,507 ligands have been completely processed and docked onto the orthosteric and allosteric site of the VDR, PXR, PPARγ, and ER𝛼. There was similarity in structure within the top PFAS among the individual receptor binding sites. For example, the majority of the top PFAS at the VDR orthosteric site displayed many similarities in their structures. However, top PFAS at differing NRs had differences in structure due to the nonidentical amino acid residues occupying the binding sites. In terms of PFAS structure, fluorinated carbon rings and chains were common in the structure of various PFAS, with some also incorporating elements such as oxygen, nitrogen, hydrogen, and sulfur, like sulfonic acid. Despite these variations, the PFAS displayed structural likeness. Examining the structural similarities and differences among prominent PFAS bindings to NR is key for identifying the properties that contribute to high binding energy and NR activation. We discuss below the PFAS screening results against the VDR, PXR, PPARγ, and ER𝛼.

3.2.1. PFAS Virtual Screening: VDR Orthosteric Site

There were 130 PFAS that bound to the orthosteric site of the VDR at a greater docking score than calcitriol. The range of PFAS docking scores with greater binding affinity than calcitriol were from -12.4 kcal/mol to -14.7 kcal/mol, as the top 15 PFAS are shown in Table 3. Additionally, at the orthosteric site, there are 21 PFAS that had similar docking scores as calcitriol, all with a docking score of -12.3 kcal/mol. Among the remaining PFAS, DTXSID101034031 stood out the most, as it had a docking score of +2.5 kcal/mol at the orthosteric site, which was the only positive binder among the 9,356 PFAS that had a lower docking score than calcitriol. The positive binding energy observed is due to that PFAS is very large, with a molecular weight of 2290.668 g/mol, and does not properly fit into the VDR (Supplementary Figures - Suppl 2). Similarly, other large PFAS have been observed through our docking study in PPARγ, PXR, and ER𝛼 displaying as positive binders. The inability of these large PFAS molecules to adequately accommodate within the binding site of the VDR could result in suboptimal interactions, leading to less favorable binding energy (Supplementary Figures - Suppl 2). The docking of the top 5 PFAS in comparison to the orthosteric endogenous ligand, calcitriol, are displayed in Figure 4. Among the top 5 orthosteric PFAS binding ligands, 4 happened to have only carbons and fluorines, and only one had a sulfonic acid, but it still was mainly of fluorines attached to a carbon chain (Supplementary Figures - Suppl 3). PFAS that had similar binding energy to calcitriol were different from the top ligands in the sense that they also contained other atoms that are not as commonly found in the top PFAS, such as nitrogen, iodine, more oxygen, sulfur, phosphate, and sulfonic acid. A major similarity between the top 5 PFAS binding ligands is the non-polarity, as there are carbon fluorine bonds decorate throughout the molecule, negating the dipole-dipole charge and neutralizing it. DTXSID40881335 (Fluoropolymer - FP), DTXSID60881337 (Fluoropolymer - FP), and DTXSID101023399 (Fluoropolymer - FP) all contain rings with fluorines attached around it causing structures that are nonpolar and stable enough to interact with the VDR, therefore giving a greater binding affinity. DTXSID701036930 (Perfluoroalkyl Sulfonic Acid- PFSA) and Hexatetracontafluorodocosane (Fluoropolymer - FP) have more linear type structures, allowing for more stereochemistry about the structures, but their main non-polarity gives a greater binding affinity when docked into the orthosteric site.

3.2.2. PFAS Virtual Screening: VDR Allosteric Site

2,229 PFAS chemicals bound to the allosteric site of the VDR with greater binding affinity than lithocholic acid. The scores of those PFAS ranged from -6.7 kcal/mol to -10.4 kcal/mol as shown in Table 4. 412 PFAS had a similar docking score to that of lithocholic acid, with a docking score of -6.6 kcal/mol. Figure 5 shows the top 5 PFAS binding to the allosteric site of the VDR adjacent to the native ligand, lithocholic acid. The top 5 allosteric binding PFAS by docking score is displayed. 4 of the 5 top PFAS have only carbons and fluorines, while one contains some oxygens, nitrogens, and hydrogen (Supplementary Figures - Suppl 4). PFAS that had similar binding energy to lithocholic acid didn’t have much structural similarities to LCA other than consisting of carbon chains and rings. DTXSID70597457 (Fluoropolymer - FP), DTXSID60881337 (Fluoropolymer - FP), DTXSID40881335 (Fluoropolymer - FP), DTXSID90984683 (Fluoropolymer - FP), and DTXSID501041016 (Perfluoroalkyl Acids - PFAAs) consist of carbon rings with fluorines attached. Similar to the orthosteric site, the allosteric site prefers nonpolar molecules.

3.2.3. PFAS Virtual Screening: PPARγ Orthosteric Site

For the orthosteric binding site, PFAS chemicals with a docking score greater than -9 kcal/mol were shortlisted given that the minimum docking score for a known ligand was -9 kcal/mol. Of the 9,505 docked PFAS compounds into the orthosteric site, 1,823 molecules have the binding affinity stronger than the endogenous ligand, ET1. The top 15 PFAS with the greatest docking scores were then compiled in Table 2. The top 15 PFAS are bound to the PPARγ orthosteric site with binding affinities ranging from -11.3 to -12.2 kcal/mol, displaying the top ranked PFAS in Figure 6 and second top ranked in Figure 7. The top ranked orthosteric PFAS had many structural similarities that were slightly different than that of the allosteric site of the PPARγ. The orthosteric PFAS included more functional groups that made the overall compound more polar, like amine groups, ethers, nitrogens, and oxygens. Some compounds contained stronger electron withdrawing groups like sulfonic acid and carbonyls, some with a mix of different functional groups, whereas one of the PFAS had no functional groups resulting in no polarity. Therefore, the inference can be drawn that the structure of the orthosteric site has more polarity than the structure of the allosteric site of PPARγ, which explains why there are more functional groups, causing polarity, present in PFAS that bind at high affinities at the orthosteric site in comparison to that at the allosteric site.

3.2.4. PFAS Virtual Screening: PPARγ Allosteric Site

For the allosteric binding site, PFAS compounds with a docking score greater than -10 kcal/mol were shortlisted given that the minimum docking score for a known ligand into the allosteric binding site was -10 kcal/mol. Of the 9505 docked PFAS compounds, 693 compounds have the binding affinity greater than -10 kcal/mol, which are stronger binders than the endogenous ligand, T35. The structure and binding interactions of the top PFAS compounds were then analyzed, and the top PFAS at the allosteric site is displayed in Figure 8. These top PFAS molecules formed polar interactions with the following amino acid residues: HIS323, TYR473, SER289, CYS285, and ARG288. Previous research indicates that HIS323, TYR473, and SER289 form H bonds with known PPARγ agonists and are crucial for PPAR activation (Gim et al 2018). The top 15 PFAS have many similarities, as they all contain some sort of cyclic features like rings and aromatic rings. Majority of the top 15 PFAS binding ligands have nonpolar substituents, however there are some with ethers, esters, alcohol, and amine groups, but come out to be nonpolar compounds.

3.2.4. PFAS Virtual Screening: PXR Orthosteric Site

The docking outcomes of PFAS molecules onto the orthosteric site of PXR were organized based on their descending order of binding affinities. Among these, 650 PFAS molecules exhibited binding affinities surpassing that of the reference ligand 4WH, which was -10.2 kcal/mol. The top ligands demonstrate the highest binding affinities shown in Table 5. Additionally, the superimposed conformations of the top ligands with the highest binding affinities shown in Figure 9. These molecules share several structural similarities that may have contributed to their notable elevation in docking. In our study, orthosteric site top ranked PFAS chemicals 1, 2, 4, 5, 9, 10, and 11 namely DTXSID501041016, DTXSID401026885, DTXSID701026989, DTXSID70896735, DTXSID801026875, DTXSID60896264, and DTXSID40881335, exhibit indications of polycyclic aromaticity, exhibiting their ability for intramolecular resonance engagement as described below (Supplementary Tables – Suppl. Table 3). Upon analysis, we noticed that π-stacking interactions involving aromatic rings play a crucial role, significantly contributing to the strength of ligand binding (Brylinski 2018). These interactions, characterized by traits like planarity and a discernible π-electron cloud above and below the rings, provide diverse interaction possibilities (Brylinski 2018). PFAS molecular binding (Figure 9 and Table 5) such as π–π interactions, cation−π, amide−π, halogen−π, and hydrogen-bond interactions facilitated by heteroatoms are encompassed by these interactions (Brylinski 2018). It has been established that π–π stacking is essential for favorable electron correlation (Rimac at al 2021). The capacity of PFAS ligands to display resonance characteristics and distribute electrons uniformly across the entire molecule is crucial for strong binding affinities, as PFAS ligand-receptor interactions are influenced by the extent of electron cloud overlap. These findings from our study emphasize the importance of aromatic π-stacking in influencing PFAS ligand-receptor interactions, providing valuable insights into the factors that contribute to strong binding affinities.

The PFAS fluorine atom’s unique characteristics enhance interactions between molecules. Its electronegativity forms strong bonds with the receptor, optimizing binding and stabilizing the ligand in the binding site. Further, the small size of fluorine atoms reduces clashes and increases PFAS ligand flexibility, leading to a more favorable alignment with the receptor, resulting in improved docking scores and potentially higher binding affinities. The significance of fluorine atoms in establishing high binding affinities between proteins and ligands finds empirical support in existing literature. This can be supported by the binding affinity of PFAS to PXR (Figure 9 and Table 5). Notably, the analysis revealed that among the top 14 PFAS exhibiting the strongest binding affinities to the PXR receptor, fluorine substituents were attached to all carbon atoms constituting the ligand's backbone, excluding the carbon atom present in the carboxylic acid group (Lai et al 2020). Of particular interest, a PFAS known for its substantial binding affinity, namely FTCA (5:3 Fluorotelomer carboxylic acid), featured a carboxylic acid functional group (Lai et al 2020). This specific ligand shared a similar chemical structure with the PFAS denoted as DTXSID701026989 (Perfluoroalkyl Acids - PFAAs). This ligand played a prominent role in our investigation and demonstrated a distinctly high binding affinity to the orthosteric site within the PXR receptor. Further, all of the top ranked PFAS that were docked at the orthosteric site of PXR had at least one fluorine substituent present, seeing a trend with ligands of higher binding affinities containing more fluorine substituents.

3.2.5. PFAS Virtual Screening: PXR Allosteric Site

The docking outcomes of PFAS molecules onto the allosteric AF-2 site of PXR were rank ordered based on their descending order of binding affinities. Among these, 9,148 PFAS molecules exhibited binding affinities surpassing that of the reference ligand glycerol, which was -2.3 kcal/mol. Our top binding ligands had affinities ranging from -7.6 kcal/mol to -9.3 kcal/mol, and the PFAS ligand with the highest binding affinity to the allosteric AF-2 site of the PXR molecule was cobalt (Supplementary Tables – Suppl. Table 4). Many similarities in regards to the structure of these molecules could be seen when analyzing the binding affinities, with the primary ones being the presence of fluorine substituents and cyclic structures in 10 and 11 namely DTXSID40204383 and DTXSID301027559 of the top 13 binding ligands respectively. As stated earlier, fluorine's distinct properties amplify molecular interactions by creating robust bonds with PXR, enhancing binding of the PFAS ligand within the site. Additionally, the small size of fluorine atoms minimizes collisions and boosts ligand flexibility, promoting a stable cooperative alignment with the receptor. This alignment potentially elevates docking scores and enhances binding affinities. Existing literature highlighted the crucial role of fluorine atoms contributing to strong binding affinities between proteins and ligands in PXR orthosteric sites, but minimal research has been done of the effects of fluorine on the allosteric site of PXR. However, for the purpose of this study, we infer that the effects will be rather similar. Further, the presence of cyclic structures has a significant effect on the ability of PFAS ligands to bind to the AF-2 allosteric site of PXR, as the AF-2's ligand-dependent groove is made up of helices 3, 4, 5, and 12, which are primarily hydrophobic regions (Chang 2009). The cyclic structures found in PFAS ligands are made up of adjacent cyclohexane molecules that are completely surrounded by Fluorine molecules whose electronegativities oppose each other, resulting in an overall nonpolar state of the PFAS ligands. Because the top ligands are nonpolar, they are more likely to exhibit a strong attraction to the AF-2 site due to the hydrophobic interaction between the ligand and receptor.

3.2.6. PFAS Virtual Screening: ER𝛼 Orthosteric Site

At the ER𝛼, only 20 PFAS bind at a greater binding affinity than the endogenous ligand estradiol, ranging from -11.2 to -14.3 (Supplementary Tables – Suppl. Table 5). For the most part about the top 10 PFAS ligands, there were very different structural features, as some were more linear, others contained more cyclic features, and few had other functional groups. The binding energy range for the top 10 was -12.2 kcal/mol to -14.3 kcal/mol, which is significantly stronger binder than that of the endogenous ligand with binding energy -11.1 kcal/mol. The majority of the PFAS had nonpolar features, showing minimal other elements other than carbons and fluorines, but took apart on many different conformations, like bicyclic features.

3.2.7. PFAS Virtual Screening: ER𝛼 Allosteric Site

In comparison to the orthosteric site of the ER𝛼, the PFAS when bound to the allosteric site exhibited a greater binding affinity. The range for the top 10 PFAS at the allosteric site was -12.3 kcal/mol to -14.8 kcal/mol (Supplementary Tables – Suppl. Table 6). Majority of the top PFAS display more structural similarities than those at the orthosteric site, as most contain cyclic features. The top 10 ligands don’t have many differing functional groups and just have carbon and fluorine bonds for the most part. The allosteric site seems to be more nonpolar than the orthosteric site, as there are more variable functional groups in the PFAS that bind to the orthosteric site rather than those that bind to the allosteric site.

3.3. Commercially Relevant PFAS Molecules Virtual Screening: VDR

Virtual screening is performed on the commercially relevant and commonly exposed PFAS only for the VDR as a test case study to understand their binding potential and impacts. PFOA and PFOS are commonly exposed PFAS class at higher concentrations to human compared to other PFAS (Wu et al 2023). PFOA and PFOS have been linked with the suppression of nuclear receptors like HNF4α (Fragki et al, 2021). Childhood diseases and diminished antibody response has been correlated with greater PFOA and PFOS exposure (von Holst et al, 2021). PFOA and PFOS have been associated with cancer through mechanisms such as DNA/RNA damage, lipid peroxidation, and cell death (Pierozan et al 2023). More specifically, PFOS has been shown to interact with lung surfactant in effects to result in lung collapse and failure, whereas PFOA has been associated with human cell apoptosis, and increasing studies are relating it to human carcinogenicity (Tsuda, 2016, Tarapore Ouyang 2021). Our analysis shows that PFOS and PFOA can display a similar or higher binding affinity to the VDR orthosteric site, regardless of its lower than calcitriol individual binding affinity. In the allosteric site, high concentration and binding strength play a factor for VDR overactivation or dysregulation caused by PFOA and PFOS. Other than PFOA and PFOS, other commercially relevant PFAS include Perfluorononanoic acid (PFNA), Tridecafluorohexane-1-sulfonic acid (PFHxS), Perfluoroheptanoic acid (PFHpA), and Perfluorodecanoic acid (PFDA). PFNA, PFDA, PFHpA, and PFHxS are similar to PFOA and PFOS, have similar effects to the orthosteric and allosteric sites of the NRs. Previous studies on PFDA reported the associated effects with increasing gastric cell proliferation and the suppression of cell senescence (Dong et al 2017). In higher concentrations, PFNA has shown to increase rates of melanoma in women and chances of uterine cancer. Germ cell tumor responses have a relationship with increased PFHxS concentration (Lin et al 2020).

3.3.1. Commercially Relevant PFAS Virtual Screening: VDR Orthosteric Site

We analyzed the binding strength of the commercially relevant PFOA and PFOS. At the orthosteric site, PFOA showed to have a docking score of -8.6 kcal/mol while PFOS is -9.6 kcal/mol. By docking score, PFOA and PFOS are not competitive to the orthosteric site of the VDR as calcitriol has a docking score of -12.3 kcal/mol. Previous research indicates PFOA and PFOS exposure doesn’t lead to dysregulation to the fullest extent on the VDR. However, that is not limited to other PFAS, as other PFAS types show greater binding affinity to the VDR than PFOS and PFOA (Singam et al 2023). Since PFOA and PFOS do not display a greater affinity for the VDR in comparison to calcitriol, the contributing factor for modified receptor function is due to its high concentration exposure. Alongside PFOA and PFOS, other common PFAS that showed significance to humans are PFDA, PFNA, PFHxS, and PFHpA. PFDA had the greatest binding affinity as the docking score was -9.7 kcal/mol, the strongest binding among the relevant PFAS (Supplementary Tables – Suppl. Table 7). Similarly, PFNA was relatively strong compared to the other relevant PFAS with a docking score of -9.6 kcal/mol, which is the same as PFOS. PFHxS and PFHpA had the lowest docking scores with -8.3 kcal/mol and -8.2 kcal/mol respectively. With lower docking scores, commercially relevant PFAS have shown to affect receptor functionality due to their high concentrations. PFDA having the greatest binding affinity along with the higher exposure possibility to humans that increases its chance to bind to the VDR. However, it is important to note that there are 130 different PFAS chemical classes that can bind to the VDR's orthosteric site, and they can also found in commercial products, though not as prevalent as the ones listed above. While the commercially significant PFAS have lower docking scores than calcitriol at -12.3 kcal/mol, the prevalence and higher exposure concentration of these relevant PFAS can induce carcinogenic effects through VDR.

3.3.2. Commercially Relevant PFAS Virtual Screening: VDR Allosteric Site

At the allosteric site, PFOS and PFOA both displayed a docking score of -6.7 kcal/mol, greater than the docking score of lithocholic acid at -6.6 kcal/mol. Lithocholic acid plays a protective role for TNF-α-induced injury of intestinal barrier by binding to the VDR (Yao et al, 2019). LCA has proven to bind at both the LBD and allosteric site, but when bind to the allosteric site, it provides more structural balance for the VDR. Our earlier study targeting the allosteric site, we showed that LCA binds to Helix 12, therefore stabilizing the active agonist conformation (Issa et al, 2022). Given the elevated concentrations of PFOA and PFOS, it can bind to the VDR allosteric site as well, altering the activity of the receptor. PFDA, PFNA, PFHxS, and PFHpA were also bind in the allosteric site, and showed more significant binding affinity than in the orthosteric site (Supplementary Tables - Suppl Table 7). PFDA and PFNA had the greatest binding affinity among these four relevant PFAS at -6.8 kcal/mol and -6.7 kcal/mol, respectively. They both displayed greater binding affinity than lithocholic acid, which had a docking score of -6.6 kcal/mol. This is significant as the activity of the VDR can be altered as LCA plays a role in regulating the VDR functions. PFDA and PFNA in high exposure can cause VDR activation or dysregulation once bound to the allosteric site, due to its greater binding affinity and concentration. PFHxS and PFHpA have docking scores of -6.4 kcal/mol and -6.5 kcal/mol, which is just under the docking score of LCA. Due to PFHxS and PFHpA having binding affinities that are only lower than LCA by 0.1-0.2 kcal/mol, their individual binding affinity is not significantly weaker than LCA. This means they can be competitive binders to the allosteric site. In combination with the high concentration of PFHxS and PFHpA, the allosteric site of the VDR is very susceptible to being regulated by commercially relevant PFAS.

3.4. Modeling of PFAS Mixture

Within the context of nuclear receptors, it is well-established that hydrophobic compounds beyond endogenous ligands have the potential to interact with the orthosteric and allosteric sites of the nuclear receptors (NR). Such interactions often result in NR activation or dysregulation, leading to disease outcomes such as cancer. NR overactivation appears to be a consequence of when nonendogenous hydrophobic chemicals engage with nuclear receptors (Issa et al, 2022). Distinct combinations of these nonendogenous chemicals binding to nuclear receptors have been identified, and recognizing there are many scenarios, we have reviewed seven combinations that are deemed significant. Due to different PFAS classes, different effects are induced, resulting in many possible PFAS binding scenarios. Firstly, a possible situation is when the endogenous ligand occupies the orthosteric site concurrently with another compound binding to the allosteric site. Another scenario occurs when an endogenous ligand binds at the allosteric site simultaneously with other compound binding to the orthosteric site. Lastly, a situation arises when both the orthosteric and allosteric sites of the receptor are engaged by nonnative exogenous ligands. Under this category of activation, there are many subcategories as different PFAS classes alter their induced effects. This scenario supports that NR overactivation can occur when two exogenous ligands simultaneously bind to the receptor, as the hydrophobicity brings stabilization, and the effect causes a super-activated receptor (Issa et al, 2022).

We developed several PFAS mixture models to represent the different possibilities of PFAS binding to the orthosteric and allosteric sites of nuclear receptors. In this study, we report only models of the VDR, and similar models can be made for other NRs. Figure 10 depicted various mixed models with seven significant combinations of PFAS and endogenous or native ligands, each contributing to the activation or altered function of the receptor. Similar to our previous environmental chemical motor fuel oil Salpn-VDR and Vitamin D3 standalone and mixed model association studies (Issa et al, 2022), both in-silico and in-vitro, we discuss several of the following PFAS-Ligand-VDR mixed model scenarios. PFAS may have toxic effects when bound simultaneously to the allosteric site in combination with the calcitriol at the orthosteric site, potentially causing an overactivation of the VDR as shown in Figure 10 panel a). In this model, LCA can be replaced by the PFAS DTXSID70597457 (Fluoropolymer - FP) at the allosteric site, while the orthosteric site is occupied by calcitriol. This form of receptor activation requires PFAS to bind to the allosteric site, which has a lower binding affinity and is more accessible because it is more superficial than the orthosteric site on the protein. This accessibility allows various PFAS to bind and potentially cause receptor overactivation.

The combination of PFAS at the allosteric site with calcitriol at the orthosteric site is a potential hazard, as the PFAS stabilizes the VDR while calcitriol performs its native function, potentially resulting in an overactivated VDR complex. PFAS can bind to nuclear receptors at the orthosteric site alongside the endogenous ligand at the allosteric site, potentially inducing carcinogenic effects. For example, as shown in Figure 10, panel b, LCA binds to the allosteric site while a high-energy binding PFAS, DTXSID40881335 (Fluoropolymer - FP), occupies the orthosteric site. These combinations are common due to the hydrophobic nature of PFAS and the abundance of endogenous ligands like calcitriol and LCA. Consequently, one endogenous ligand will bind to the receptor while the other site is occupied by a PFAS, leading to the overactivation of the VDR. In scenarios of high PFAS exposure, two PFAS molecules could bind to both the orthosteric and allosteric sites, as illustrated in Figure 10, panel c, with two high-affinity binding PFAS, including DTXSID40881335 at the orthosteric site and DTXSID70597457 at the allosteric site. Both PFAS molecules, consisting only of carbon-fluorine rings and chains, enhance the hydrophobic interaction, which is highly carcinogenic due to the non-endogenous nature of these chemicals, resulting in significant protein stabilization and super-activation. Particularly concerning are PFAS commonly found in humans, such as PFOA and PFOS, which are often found at high concentrations. As depicted in Figure 10, panels d through g, these chemicals can occupy both the orthosteric and allosteric sites of the VDR, leading to overactivation that could lead to carcinogenesis. This is significant given the prevalent rates of PFOA and PFOS in humans, which could lead to VDR overactivation and potential cancer risks. These super activation mechanistic scenarios are similarly applicable to other nuclear receptors such as ERα, PXR, and PPARγ.

Our previous investigation demonstrated the super activation of the VDR by the motor fuel oil compound Salpn, which binds to the allosteric site, enhancing stability and inducing excessive agonistic activity when combined with calcitriol (Issa et al 2022). Similarly, PFAS molecules, due to their hydrophobic nature similar to endogenous hormone ligands, exhibit strong affinities for both allosteric and orthosteric sites of nuclear receptors. This interaction has been implicated in the overproduction of hormones, particularly in the case of ERα, potentially leading to carcinogenic effects similar to those observed with exposure to carcinogens (Qui et al 2020). PFAS have shown effects like altering receptor function when bound with calcitriol to the VDR creating a super-agonistic effect (Fenton et al, 2021) (Issa et al 2022). The LBD orthosteric site has a large volume that can accommodate both calcitriol and PFAS, presenting another possibility of overactivation where both ligands occupy the orthosteric site; however, this scenario is not explored in the context of this study. Nevertheless, the structural similarities between the orthosteric and allosteric sites of nuclear receptors allow for the prediction that PFAS could have similar overactivation effects on other nuclear receptors as observed in the VDR by Salpn (Issa et al 2022) (Fenton et al, 2021).

4. Strengths and Limitations

Our PFAS virtual screening study is distinct from others as it investigates over 9,000 PFAS chemical interactions with PPARγ, PXR, VDR, and ERα, both at their orthosteric and allosteric sites. The interactions between PFAS and the allosteric sites of these nuclear receptors (NRs) are particularly a strength, as they include more possibilities for receptor activation. Further, the inclusion of allosteric sites is significant, as our study demonstrates that PFAS can readily bind these sites due to easier access and lower binding potential. In NRs, the allosteric site plays a major role in maintaining receptor stability, contributing to homeostatic responses. However, when occupied by exogenous ligands, hyperactivation can occur, underscoring the importance of studying the allosteric site (Ishizawa et al., 2018) (Issa et al., 2022) (Fenton et al, 2021). We developed mixture models of PFAS in combination with other PFAS and endogenous ligands to illustrate how NRs can be overactivated. Binding of PFAS to either the orthosteric or allosteric sites of nuclear receptors can result in hyperactivation, similar to the effects observed with carcinogenic compounds. Different classes of PFAS may also alter the activation status of the nuclear receptor. For example, in the context of the VDR, simultaneous binding of PFAS and endogenous ligands can lead to receptor stability and hyperactivation, reflecting carcinogenic pathways. In our earlier Salpn-VDR study, Salpn has been associated with the overactivation of other nuclear receptors, including PPARγ, androgen receptor (AR), ERα, and thyroid hormone receptor alpha (TRα) (Issa et al., 2022). The hydrophobicity of Salpn, similar to that of PFAS, supports the conclusion that PFAS can have analogous effects on nuclear receptors when bound. It is noteworthy that some of these mixed model scenarios are hypothetical and require further validation. However, these mixture models provide valuable representations of how the combinatorial effects of PFAS can lead to different types of functional alterations in the receptor. In this study, though we included local binding site flexibility, however, future investigations are needed to include full flexibility of protein and ligand complex. Fully flexible docking will offer more precise predictions and provide additional insights into how PFAS can bind with NR complexes. These full flexible docking studies are our future focus observing interactions between PFAS and nuclear receptors (NRs) under flexible conditions. Further improvements to this study would involve employing MM-GBSA or MM-PBSA methods to assess the full-fledged binding free energy of top PFAS ligands to the NRs. This approach would enhance the accuracy of the computational model. Additionally, expanding the scope of the study beyond PFAS to understand how other chemicals, both standalone and in combination with PFAS, can induce carcinogenic responses is important. This broader examination is important for understanding the full range of toxicological impacts and potential health risks posed by chemical interactions with nuclear receptors.

5. Future Directions

Our ongoing studies are to understand how PFAS and other endocrine disrupting chemicals (EDCs) in combination, can cause overactivation of NRs. This can be performed by the mixture modeling, as we performed a case study for the VDR and earlier mixed model study (Issa et al 2022), for EDCs and PFAS which will display how they can cause overactivation simultaneously. EDCs such as polybrominated diphenyl ethers (PBDEs), microplastics (MPs), nanoplastics (NPs), and pesticides including PFAS pesticides, and nonpersistent EDCs like parabens, phenols, and phthalates have shown carcinogenic effects (Wan et al 2021). Analyzing the combinatorial effect of PFAS and EDCs are relevant due to the carcinogenic nature of both chemicals. Similar to PFAS, EDCs and nonpersistent EDCs have been correlated with exposure to humans through consumption, and causing nuclear receptor dysregulation (Yilmaz et al 2020). Studying PFAS alongside PFAS pesticides and EDCs holds importance for future research, given their common presence in humans, according to the latest report by the Pesticide Action Network (PAN). Their potential to induce cancerous effects emphasizes the significance of the analysis of their combinatorial effect. Studies has shown EDCs effect on nuclear receptors such as PPARγ, ERα, and PXR, to have detrimental effects at low concentrations (Toporova and Balaguer 2020), raising the concern when combined with PFAS, with disrupting capabilities to nuclear receptors, resulting in super activation and carcinogenesis. Therefore, analyzing the effects PFAS, PFAS pesticides, and EDCs induce together is important for understanding how they can cause carcinogenic effects.

6. Conclusion

We performed computational docking simulations followed by binding affinity analyses to evaluate the binding strength of PFAS with nuclear receptors such as VDR, PXR, PPARγ, and ERα. The binding affinities were quantified in kcal/mol, and PFAS were subsequently ranked based on these values. Numerous PFAS exhibited binding energies higher than corresponding endogenous or native ligands. Binding strength analysis reveals that 130 PFAS bind strongly than the orthosteric ligand calcitriol to VDR and 2,229 outcompeting lithocholic acid (LCA) at the allosteric site. In the case of PPARγ, 1,863 PFAS bind with greater affinity than orthosteric ligand ET1 and 693 binds stronger than the allosteric site ligand T35. In the case of PXR, 650 PFAS had affinity greater than the orthosteric ligand 4WH, while 9,148 had higher affinity than allosteric ligand glycerol. Finally, for the ERα, 40 PFAS had a greater binding affinity than the native ligand estradiol, while 8 had higher affinity than the allosteric ligand SRC-1. These stronger PFAS binding to the nuclear receptors, influenced by inherent hydrophobic interactions between PFAS and the ligand binding site of the nuclear receptors. Despite diverse interactions of PFAS across different receptors due to varied amino acid residues within the binding domains, commonalities persist in the hydrophobic nature of these interactions. Notably, the top-binding PFAS at orthosteric sites of PXR, VDR, and ERα were predominantly hydrophobic, consisting mainly of carbon and fluorine, whereas top binders to PPARγ included additional elements, reflecting less hydrophobicity in the ligand-binding domain (LBD) of PPARγ compared to other receptors. Unlike other NRs studied here, PPARγ had a smaller count of stronger binding PFAS in comparison to the endogenous ligand at the allosteric site than the orthosteric site.

Further investigation on the binding affinity of commercially relevant PFOA and PFOS, towards the VDR, demonstrate that PFOA, PFOS do not effectively displace or compete with calcitriol, at the orthosteric site, a primary binding site on the VDR. In contrast, PFAS can compete with the allosteric ligand lithocholic acid (LCA), a bile acid, and can easily displace it from its binding position due to its relatively lower affinity. Additionally, PFAS structurally mimic the bile acid and exhibit highly similar physicochemical properties, forming stronger binding interactions with the VDR at the allosteric site. This displacement allows the stable PFAS to bind more readily to secondary, allosteric sites on the VDR. Such allosteric binding does not hinder the primary activation of the receptor by calcitriol but contributes to the additional affinity and influence receptor conformation leading to the activation. Moreover, the higher concentrations of these PFAS significantly enhance their likelihood of interacting with the VDR. This concentration-dependent binding suggests that in environments with elevated PFAS levels, such as in industrial or contaminated sites, there could be an increased potential for these substances to interact with the VDR or other nuclear receptors, potentially leading to altered receptor activity and downstream biological effects. We are currently expanding our analysis to other receptors, given that these PFAS are frequently exposed at higher concentrations with potential for activating proteins and disrupting their function. Given the association of nuclear receptor or other protein dysregulation with a range of adverse health outcomes, including cancer, osteoporosis, obesity, and metabolic disorders, highlight the importance of our findings. These results support the necessity for direct in vitro assessments and further in-depth analyses to better understand the implications of PFAS exposure on human health.

Through mixture modeling, we developed several mixture combinations of PFAS, calcitriol, and LCA against the nuclear receptors to analyze simultaneous binding potential at their orthosteric and allosteric sites. In the case of VDR, the mixed models suggest that some PFAS can hyperactivate the VDR by the synergistic effects through simultaneous binding at their orthosteric and allosteric sites, potentially leading to carcinogenic effects. Further, these models not only highlight potential mechanistic insights of carcinogenic potential but also shed light on various other functional alterations induced by PFAS. As both endocrine disrupting chemicals (EDCs) and other commercially relevant chemicals co-exposed with PFAS through dietary and environmental exposure such as water sources, understanding their combined effects is essential for developing strategies to mitigate their impact. The continued studies into the combinatorial effects of PFAS are ongoing. This is important, specifically, given the persistent, and rising environmental concentrations of these substances and their hazardous effects to humans.