Combined Effects of Multiple Per- and Polyfluoroalkyl Substances Exposure on Allostatic Load Using Bayesian Kernel Machine Regression

This study aims to investigate the combined effects of per- and polyfluoroalkyl substances (PFAS) on allostatic load, an index of chronic stress that is linked to several chronic diseases, including cardiovascular disease and cancer. Using data from the National Health and Nutrition Examination Survey (NHANES) 2007–2014, this study examines the relationship between six PFAS variables (PFDE, PFNA, PFOS, PFUA, PFOA, and PFHS) and allostatic load using Bayesian Kernel Machine Regression (BKMR) analysis. The study also investigates the impact of individual and combined PFAS exposure on allostatic load using various exposure-response relationships, such as univariate, bivariate, or multivariate models. The analysis reveals that the combined exposure to PFDE, PFNA, and PFUA had the most significant positive trend with allostatic load when it was modeled as a binary variable, while PFDE, PFOS, and PFNA had the most significant positive trend with allostatic load when modeled as a continuous variable. These findings provide valuable insight into the consequences of cumulative exposure to multiple PFAS on allostatic load, which can help public health practitioners identify the dangers associated with potential combined exposure to select PFAS of interest. In summary, this study highlights the critical role of PFAS exposure in chronic stress-related diseases and emphasizes the need for effective strategies to minimize exposure to these chemicals to reduce the risk of chronic diseases. It underscores the importance of considering the combined effects of PFAS when assessing their impact on human health and offers valuable information for policymakers and regulators to develop strategies to protect public health.


Introduction
Exposure to multiple contaminants concurrently is a fact of life [1]. Exposure to individual toxic substances, such as Perfluorooctanesulfonic acid (PFOS), has been shown to have adverse health outcomes. Yet, little is understood about how exposure to multiple contaminants affects health [2]. This is especially the case for the exposure to multiple contaminants on the stress response system, which potentially potentiates the development of chronic diseases. It is thus critical to explore the effect of multiple contaminants such as Per-and polyfluoroalkyl substances (PFAS) on stress.
PFAS have been utilized extensively in consumer goods for more than 70 years. The utilization of PFAS in commercial and industrial settings is due to their special thermal stability and surface activity characteristics, including hydro-and lipophobicity [3].
PFAS have been used as lubricants in industrial processes, as additives in pesticides and pharmaceuticals, and as lubricants in food packaging, firefighting foams, and coatings to produce nonstick and stain-resistant qualities [4]. The negative consequences of exposure to multiple PFAS could be deleterious and may have long-term effects on the life course of affected peoples [5,6]. exposures in epidemiological and toxicological studies, making BKMR a valuable tool for researchers [13,18].
Through the use of BKMR, this study hypothesizes that there exists a significant association between exposure to multiple PFAS simultaneously and high levels of allostatic load. The primary objectives of this study are to determine the extent of association between exposure to multiple PFAS and allostatic load, identify the specific PFAS that contribute the most to the observed association, and investigate these factors by modeling allostatic load as a categorical and continuous variable.
The findings of this study provide critical insights into the effects of multiple PFAS exposure on human health and contribute to the public health literature aimed at understanding and reducing the risk of PFAS exposure. Furthermore, identifying specific PFAS that contribute the most to the observed association can help prioritize regulatory actions and interventions aimed at reducing exposure to these harmful chemicals.

Description of Cohort
This investigation utilized data from the National Health and Nutrition Examination Survey (NHANES), administered by the Centers for Disease Control and Prevention (CDC). The NHANES collected biological specimens and survey data on demographics, social, and clinical factors for the civilian, non-institutionalized population residing in any of the 50 states of the United States of America and the District of Columbia (DC) using a complex, multistage probability design. This study used data for adults aged 20 years and older.
For this study, data were assessed over four cycles: 2007-2008, 2009-2010, 2011-2012, and 2013-2014. The final sample size for our study was n = 23,482 individuals, consisting of 11,414 males and 12,068 females. The study aimed to investigate several variables of interest, including exposure to specific perfluoroalkyl substances (PFAS), such as Perfluorodecanoic acid (PFDE), Perfluorononanoic acid (PFNA), Perfluorooctanesulfonic acid (PFOS), Perfluoroundecanoic acid (PFUA), Perfluorooctanoic acid (PFOA), and Perfluorohexanesulfonic acid (PFHS). Additionally, the study also examined the allostatic load index, which was operationalized using the following clinical and biomarkers: Systolic Blood Pressure (SBP), Diastolic Blood Pressure (DBP), triglycerides, HDL cholesterol, total cholesterol, C-reactive protein (CRP), Body Mass Index (BMI), hemoglobin A1C, albumin, and creatinine clearance. Overall, this study utilized a large and diverse sample size, as well as a comprehensive set of variables of interest, to investigate the relationships between PFAS exposure and allostatic load.

Blood Serum Measurements
As part of the NHANES, participants visited a mobile examination center (MEC) where they provided blood samples for laboratory analysis. The blood samples were used to assess a wide range of health indicators, including markers of exposure to environmental contaminants such as per-and polyfluoroalkyl substances (PFAS).
To ensure that laboratory measurements were conducted under comparable circumstances, the MEC's controlled environment was used. Blood samples were collected in polypropylene and polyethylene containers, and depending on the participant's age, at least 0.5 mL of serum was processed. Samples were then frozen or refrigerated before being sent to labs across the country.
For this study, participants who were scheduled for a morning session were required to fast for nine hours before the blood draw. After the initial blood draw, they were instructed to ingest 75 g of dextrose (10 oz of glucose solution) within ten minutes. Two hours later, a second blood sample was collected.

PFAS Extraction and Quantitation
To extract and quantify PFAS from human serum, a 50 L aliquot was first diluted with formic acid and injected into a commercial column switching system for concentration on a solid-phase extraction column. This method of extraction was used to separate the analytes from other serum constituents. High-performance liquid chromatography was then employed to separate the analytes from each other. Detection and quantification of the PFAS was achieved using tandem mass spectrometry, which converts liquid phase ions into gas-phase ions. To detect and quantify the PFAS, an electrospray ionization source was modified to produce negative-ion Turbo Ion Spray (TIS) ions. This method has a low limit of detection (LOD) in the low parts per billion range, which enables the rapid detection of these PFAS in human serum. Overall, this method provides a reliable and sensitive approach for the detection and quantification of PFAS in human serum.
All PFAS analytes in the dataset had a constant detection limit of 0.10 ng/mL. For each analyte, two variable names were provided in NHANES: a value of "0" indicated that the result was at or above the detection limit, while a value of "1" indicated that the result was below the detection limit. In cases where the analytic result was below the lower limit of detection, an imputed fill value was inserted into the analyte results field. This value was calculated as the lower limit of detection divided by the square root of 2 (LLOD/sqrt(2)), which was 0.10/ √ 2 = 0.07. Consequently, the LOD for each PFAS was either 0.10 or 0.07 [19,20].

Quantifying Allostatic Load
The markers of interest that went into the allostatic load were: Systolic blood pressure-SBP, diastolic blood pressure-DBP, triglycerides, high-density lipoprotein (HDL) cholesterol, and total cholesterol, inflammatory (C-reactive protein-CRP), and the metabolic systems (body mass index-BMI, hemoglobin A1C, albumin, and creatinine clearance and they were quantified as follows. The NHANES procedure for obtaining blood pressure involved collecting a maximum of 3 brachial systolic and diastolic BP readings for each participant using appropriate cuff sizes based on midarm circumference. Trained medical personnel followed a standard protocol and used a Baumanometer true gravity mercury wall model at the MEC. Triglycerides and HDL-C were measured using the Roche modular P chemistry analyzer (Roche Diagnostics, Indianapolis, IN, USA). CRP was measured via latex enhanced nephelometry and quantified using anti-CRP antibodies with a hydrophilic shell and polystyrene core of CRP particles. BMI was calculated as weight (kg) divided by the square of standing height (m 2 ), while glycohemoglobin was measured using the A1c G7 HPLC Glycohemoglobin Analyzer (Tosoh Medics, Inc., San Francisco, CA, USA). The Roche/Hitachi Modular P Chemistry Analyzer was used to measure urine creatinine (Roche Diagnostics, Indianapolis, IN, USA). Finally, a non-competitive, double-antibody method was employed for measuring human urinary albumin, involving covalently attaching an antibody to human albumin to polyacrylamide beads, reacting the solid-phase antibody with a urine specimen, and measuring the fluorescence of the resulting complex using a fluorometer.

Operationalizing Allostatic Load
Influenced by earlier research [2,[21][22][23][24][25], the study operationalized allostatic load as a cumulative index of physiological dysfunction in the cardiovascular, inflammatory, and metabolic systems. This index included multiple markers, such as triglycerides, highdensity lipoprotein (HDL) cholesterol, total cholesterol, and systolic and diastolic blood pressure, along with other measures like body mass index, hemoglobin A1C, albumin, and creatinine clearance. Quartiles were used to divide the markers based on their distribution within the dataset. Specifically, the top quartile was considered high-risk for all markers except for albumin, creatinine clearance, and HDL cholesterol, for which the lowest 25% of the distribution was deemed to be the highest risk based on prior literature [26][27][28][29][30][31][32]. Each study participant was given a value of 1 if they were in the high-risk category or a value of 0 if they were in the low-risk category for each marker which was then added up to get a total allostatic load value out of 10. Using a high-score criterion, we divided the allostatic load scores into two categories [33]. An allostatic load greater than or equal to 3 was required for a score to be considered high; otherwise, it was considered low.

Statistical Analysis
This study aimed to investigate the impact of six per-and polyfluoroalkyl substances (PFAS)-PFDE, PFNA, PFOS, PFUA, PFOA, and PFHS-on allostatic load, which was the dependent variable. Allostatic load, a measure of physiological dysregulation, was assessed using established biomarkers. The six PFAS were considered the independent variables, and we aimed to model their combined effect on allostatic load through statistical analysis. This study employed BKMR analysis to explore the relationship between the independent and dependent variables and assess the magnitude and significance of their associations.
2.6.1. Statistical Results Analysis Using BKMR

BKMR Analysis
We used BKMR with the hierarchical variable selection method due to highly correlated variables and collinearity of key variables within the dataset. In this study we built BKMR models in the R program via the R package (bkmr). BKMR models evaluated the role of mixtures or multipollutant exposures (e.g., PFOA, PFOS, etc.) on allostatic load using the kmbayes function.
2.6.2. BKMR Modeling for the Binary Outcome (Allostatic Load as 0 and 1) Using BKMR we modeled the effect of multiple PFAS on allostatic load (as a categorical variable) adjusting for covariates "Gender", "Age", "Cigarette smoking", "Physical activity", "Ethnicity", "Occupation", "Income", "Alcohol Consumption", "Education", "Birthplace", and "Time in the US" in our primary analysis and then in subsequent analysis adding lead (Pb) to the covariates in order to account for the demonstrated association of blood lead levels with allostatic load.
Binary outcomes were included in the BKMR package using the probit model via a general linear model for efficiency and convenience of computation and to overcome some of the issues that may arise in the dataset, such as collinearity under Bayesian inference [14,34]. Finally, we modeled allostatic load as a continuous variable to see if the associations held.

Critical PFAS Values Elevating Allostatic Load Is US Adults
The overall analytical plan for the first part of this study was to assess the combined effect of PFDE, PFNA, PFOS, PFUA, PFHS, and PFNA on allostatic load (as a categorical variable) while adjusting for the following covariates: Gender, Age, Cigarette smoking, Physical activity, Ethnicity, Occupation, Income, Alcohol Consumption, Education, Birthplace, and Time in the US. Figure 1 and Table 1 below present the PIPs which measure the percentage of the data that backs the inclusion variables in the model. In other words, it quantifies a variable's importance to be included within the model. The exposures that were included in our model based on their posterior inclusion probabilities (PIPs) values were PFDE, PFNA, and PFUA.
The posterior inclusion probabilities (PIPs) of selected variables in the model range from 0 to 1, where a higher value indicates a stronger association with the outcome variable, which is allostatic load in this study. PIPs close to 1 suggest that the variables are important predictors of allostatic load and play a significant role in the model. Conversely, variables with PIPs closer to 0 are less relevant to the outcome of interest and do not impact allostatic load significantly. These results are presented in Table 1 5 is PIP (posterior inclusion probability), which quantifies the importance of the variable in the variable selection. The posterior inclusion probabilities (PIPs) of selected variables in the model range from 0 to 1, where a higher value indicates a stronger association with the outcome variable, which is allostatic load in this study. PIPs close to 1 suggest that the variables are important predictors of allostatic load and play a significant role in the model. Conversely,  The study firstly compares different models (BKMR, Oracle, Linear, and True) and finds that the BKMR model is suitable for analyzing multipollutant exposures. The BKMR model provides a better fit than other models for both individual variables and the overall model, as shown in Table 1. Figure 2 displays the univariate relationship between the allostatic load variable and each individual PFAS included in the model, while holding constant the other exposures at their median values and the covariates at fixed levels. The analysis reveals that some of the PFAS variables did not show significant associations with the outcome allostatic load.
Notably, the results suggest that exposures to PFDE, PFNA, and PFUA are significantly associated with allostatic load, while other variables did not show significant associations.
The study firstly compares different models (BKMR, Oracle, Linear, and True) and finds that the BKMR model is suitable for analyzing multipollutant exposures. The BKMR model provides a better fit than other models for both individual variables and the overall model, as shown in Table 1. Figure 2 displays the univariate relationship between the allostatic load variable and each individual PFAS included in the model, while holding constant the other exposures at their median values and the covariates at fixed levels. The analysis reveals that some of the PFAS variables did not show significant associations with the outcome allostatic load. Notably, the results suggest that exposures to PFDE, PFNA, and PFUA are significantly associated with allostatic load, while other variables did not show significant associations.   Table 2 demonstrate the relationships between exposures and the outcome variable (allostatic load) after adjusting for covariates (previously described as confounders). These associations reveal that certain exposures, such as PFDE, PFNA, and PFUA, are linked to allostatic load, with varying levels of impact indicated by the steepness of the incline in the figure. The uphill and downhill sections of the figure represent high and low levels of exposure, respectively. For example, PFDE has a mean value of 2.40 for high allostatic load and −4.00 for low allostatic load, while PFHS has a mean value of 1.10 for high allostatic load and −1.03 for low allostatic load. These values increase or decrease based on the exposure levels and the impact of cofactors that affect the exposure.   Table 2 demonstrate the relationships between exposures and the outcome variable (allostatic load) after adjusting for covariates (previously described as confounders). These associations reveal that certain exposures, such as PFDE, PFNA, and PFUA, are linked to allostatic load, with varying levels of impact indicated by the steepness of the incline in the figure. The uphill and downhill sections of the figure represent high and low levels of exposure, respectively. For example, PFDE has a mean value of 2.40 for high allostatic load and −4.00 for low allostatic load, while PFHS has a mean value of 1.10 for high allostatic load and −1.03 for low allostatic load. These values increase or decrease based on the exposure levels and the impact of cofactors that affect the exposure.  Figure 3 below displays the results of probit regression modeling with risk difference (RD) inferred from the output. Here we calculated the conditional RD comparing the probability of allostatic load when the second exposure was set to its 75th percentile versus its 50th percentile (with other exposures at their median value), for exposures fixed at their 25th, 50th, and 75th percentiles. Overall, the estimation (est) point and true values were 0.012 and 0.024 for the 25th, 0.014 and 0.030 for the 50th, and finally, 0.017 and 0.031 for the 75th percentile. Overall, the data suggests that higher exposure to multiple PFAS increases the risk of higher allostatic load.
(RD) inferred from the output. Here we calculated the conditional RD comparing the probability of allostatic load when the second exposure was set to its 75th percentile versus its 50th percentile (with other exposures at their median value), for exposures fixed at their 25th, 50th, and 75th percentiles. Overall, the estimation (est) point and true values were 0.012 and 0.024 for the 25th, 0.014 and 0.030 for the 50th, and finally, 0.017 and 0.031 for the 75th percentile. Overall, the data suggests that higher exposure to multiple PFAS increases the risk of higher allostatic load.

BKMR Modeling for the Binary Outcome (Allostatic Load as 0 and 1) with Lead (within the Covariates)
The modeling below was repeated with lead (Pb) added as a covariate to see if the association between allostatic load and blood lead levels significantly altered the outcomes of the study. Figure 4 and Table 3 below present the PIPs for the exposures. The exposures to be included in the model based on the resulting PIP values are PFDE, PFNA, and PFUA. These PFAS are the same variables that were included in the prior models which did not have lead (Pb) as a covariate.

BKMR Modeling for the Binary Outcome (Allostatic Load as 0 and 1) with Lead (within the Covariates)
The modeling below was repeated with lead (Pb) added as a covariate to see if the association between allostatic load and blood lead levels significantly altered the outcomes of the study. Figure 4 and Table 3 below present the PIPs for the exposures. The exposures to be included in the model based on the resulting PIP values are PFDE, PFNA, and PFUA. These PFAS are the same variables that were included in the prior models which did not have lead (Pb) as a covariate.        Table 4 provides further details on these associations, showing the mean values of each PFAS for high and low allostatic load. The relationships between the exposures and allostatic load vary in terms of their steepness, indicating different levels of exposure. The values can either increase or decrease depending on the amount of exposure and the impact of cofactors. The uphill and downhill sections of Figure 5 and Table 4 represent high and low levels of exposure, respectively, and illustrate these relationships. Overall, these results highlight the complex associations between PFAS exposures and allostatic load, which can be further explored using multivariate models. and downhill sections of Figure 5 and Table 4 represent high and low levels of exposure, respectively, and illustrate these relationships. Overall, these results highlight the complex associations between PFAS exposures and allostatic load, which can be further explored using multivariate models.

BKMR Modeling for the Continuous Outcome (Allostatic Load as Countinuous Variable)
To investigate the relationship between multiple PFAS and allostatic load, we modeled allostatic load as a continuous variable rather than a binary variable, as the latter approach could have missed important nuances in the relationship. We employed the Bayesian Kernel Machine Regression method to evaluate the impact of multiple PFAS exposures on allostatic load, while controlling for covariates such as gender. This approach allowed us to comprehensively analyze the associations between PFAS and allostatic load and examine potential nonlinear and interactive effects. The models were adjusted for relevant covariates to ensure the accuracy and reliability of the results.
The following three figures (Figures 7-9) illustrate the variations of parameters/values throughout the sample runs. Figure 7 displays the beta1 (first exposure) value, which is equal to −0.18. Beta1 was selected to represent the sample relationship, and the figure demonstrates how the relationships between study variables change throughout the simulation runs, from zero to the end (5000 iterations in this case). Figure 8 depicts an average epsilon (ε)-(sigsq.eps) value of 1.7, which measures the error

BKMR Modeling for the Continuous Outcome (Allostatic Load as Countinuous Variable)
To investigate the relationship between multiple PFAS and allostatic load, we modeled allostatic load as a continuous variable rather than a binary variable, as the latter approach could have missed important nuances in the relationship. We employed the Bayesian Kernel Machine Regression method to evaluate the impact of multiple PFAS exposures on allostatic load, while controlling for covariates such as gender. This approach allowed us to comprehensively analyze the associations between PFAS and allostatic load and examine potential nonlinear and interactive effects. The models were adjusted for relevant covariates to ensure the accuracy and reliability of the results.
The following three figures (Figures 7-9) illustrate the variations of parameters/values throughout the sample runs. Figure 7 displays the beta1 (first exposure) value, which is equal to −0.18. Beta1 was selected to represent the sample relationship, and the figure demonstrates how the relationships between study variables change throughout the simulation runs, from zero to the end (5000 iterations in this case). Figure 8 depicts an average epsilon (ε)-(sigsq.eps) value of 1.7, which measures the error of the true regression line. Lastly, Figure 9 displays the correlation value (r1), which is 0.71. r1 reflects the relationship between PFNA and allostatic load and how the parameters change when the model is executed. Table 5 shows all the r1 values from the models.       Table 6 and Figure 10 below show how th interact, highlighting the effects of single-variable outcomes when ot fixed at different percentiles. The values demonstrate the associatio exposure and allostatic load, indicating varying effects at different per column shows the overall risk values of changes in a single pollutan    Table 6 and Figure 10 below show how the study variables interact, highlighting the effects of single-variable outcomes when other variables are fixed at different percentiles. The values demonstrate the association between PFAS exposure and allostatic load, indicating varying effects at different percentiles. The last column shows the overall risk values of changes in a single pollutant when all other pollutants are fixed and compared from their 25th to 75th percentiles. Note: a1 PIP is posterior inclusive probability for the exposures. a2 is univariable or single-independent variable health risk on the response when each of the exposures is fixed at the 25th, 50th and 75th percentile or the change effect value. a3 single-variable health risks when on the response variable by comparing the changes on each exposure for its 25th to its 75th percentiles or the result of subtracting the value of 25th percentiles from 75th percentiles. Note: a1 PIP is posterior inclusive probability for the exposures. a2 is univariable independent variable health risk on the response when each of the exposures is fixed at 50th and 75th percentile or the change effect value. a3 single-variable health risks wh response variable by comparing the changes on each exposure for its 25th to its 75th perc the result of subtracting the value of 25th percentiles from 75th percentiles.     For example, in 2007-2008, the lower quartile, median, and upper quartile concentrations of PFDE were 0.070, 0.090, and 0.200, respectively. The interquartile range was 0.130, indicating that the middle 50% of the data was between 0.025 and 0.155. The same pattern can be observed for the other PFAS variables in each cycle.
Overall, Table 7 provides valuable information on the distribution of the PFAS data, which is important for understanding the study results and drawing conclusions. Figure 11 depicts the univariate exposure-response relationship, where each individual PFAS and its association with allostatic load are shown while fixing the remaining exposures to their median (50th percentile), as presented in Table 6 above. The covariate, gender, is held constant in this analysis. The figure illustrates the associations between exposures and the allostatic load after adjusting the model for gender as a covariate. For example, the PFDE exposure is associated with a steeper increase in allostatic load, indicating varying levels of exposure, and the uphill trend in the figure represents a higher level of exposure. This suggests that the values of allostatic load increase with increasing levels of PFDE exposure. Figure 12 displays the bivariate relationship between 2-way PFAS exposures and allostatic load, focusing on pairs of PFAS-allostatic load interactions while fixing the remaining exposures at their 50th percentile (median). This approach was taken based on the association of PFAS with allostatic load and the use of PIPs. In this analysis we excluded any exposure not associated with the outcome in this bivariate relationship. Figure 13 depicts an investigation of bivariate or 2-way interactions, which examines the relationship between two exposures and their corresponding response variables. This analysis examines the effect of one exposure on the response variable while holding the other exposure fixed at different percentiles (namely, the 10th, 50th, and 90th percentiles). Figure 14 below shows the visualization of 2-exposures-response or three-way interactions by fixing the third variable (exposure) at the 10th, 50th, and 90th percentiles.    Figure 13 depicts an investigation of bivariate or 2-way interactions, which examines the relationship between two exposures and their corresponding response variables. This analysis examines the effect of one exposure on the response variable while holding the other exposure fixed at different percentiles (namely, the 10th, 50th, and 90th percentiles).      Figure 15 measure the total effect of all exposures or mixtures. The exposures are fixed at different quantities starting from the 25th percentile to 75th percentile at increments of 5 using the 50th percentile (median) to compare the exposures. The estimation for all exposures at the 50th percentile shown at zero (red dashed line) demonstrates that when all exposures are fixed at the 75th percentile, the allostatic load decreases by 0.035 units, and for the 40th percentile, the allostatic load increases by 0.002   Figure 15 measure the total effect of all exposures or mixtures. The exposures are fixed at different quantities starting from the 25th percentile to 75th percentile at increments of 5 using the 50th percentile (median) to compare the exposures. The estimation for all exposures at the 50th percentile shown at zero (red dashed line) demonstrates that when all exposures are fixed at the 75th percentile, the allostatic load decreases by 0.035 units, and for the 40th percentile, the allostatic load increases by 0.002 units.  Figure 15. Summary of overall health effects of the exposures (multimixers) on the outcome depe on various percentiles (form 25th to 75th percentiles). Note: "est" is the estimation value multivariable effects on allostatic load (e.g., the estimation of overall effects of all exposures allostatic load is 0.009 when are fixed on the 60th percentile as compared to when all of them ar the 50th percentile (median)).  Figure 15. Summary of overall health effects of the exposures (multimixers) on the outcome depends on various percentiles (form 25th to 75th percentiles). Note: "est" is the estimation value of multivariable effects on allostatic load (e.g., the estimation of overall effects of all exposures on allostatic load is 0.009 when are fixed on the 60th percentile as compared to when all of them are at the 50th percentile (median)). Figure 15 demonstrates the overall summary of the effect of multiple PFAS on allostatic load. The collective degrees across the PFAS and allostatic load both increase and decrease simultaneously. In general, the overall summary of the exposure effects demonstrates that allostatic load decreases both with increases and decreases of the percentiles in and around the 50th percentile, with some differences in a direction depending on the degrees of rising or diminishing.

Discussion
In this study, we used BKMR to assess the relationships between allostatic load and six PFAS combinations in a sample of US adults from the NHANES. The concentration of various PFAS were relatively homogenous across the four cycles of data making the data ideal for our study. Specifically, the findings from the study reveal important insights into the concentrations and time trends of the PFAS of interest in the study cohort. The results indicate that the concentrations of PFAS were relatively stable over time, with only small variations observed.
BKMR modeling is a statistical technique that is employed to evaluate the significance of variables based on their posterior inclusion probability, and the multivariable relationships between the predictor variables and outcome variables, including univariate, bivariate, and multivariate exposure-response relationships, as well as their effects on one another. This method was used to identify the most important PFAS and their interactions in producing high allostatic load, while accounting for potential confounders and covariates. By considering the relationships between multiple variables simultaneously, BKMR modeling enabled a more comprehensive understanding of the Figure 16. Univariate exposure of health effects at 95% confidence intervals comparing the changes in each exposure from the 25th, 50th, and 75th percentiles when the remaining exposures are fixed at the similar percentiles. Note: "est" is the estimation value of the univariate effects of PFAS on allostatic load.

Discussion
In this study, we used BKMR to assess the relationships between allostatic load and six PFAS combinations in a sample of US adults from the NHANES. The concentration of various PFAS were relatively homogenous across the four cycles of data making the data ideal for our study. Specifically, the findings from the study reveal important insights into the concentrations and time trends of the PFAS of interest in the study cohort. The results indicate that the concentrations of PFAS were relatively stable over time, with only small variations observed.
BKMR modeling is a statistical technique that is employed to evaluate the significance of variables based on their posterior inclusion probability, and the multivariable relationships between the predictor variables and outcome variables, including univariate, bivariate, and multivariate exposure-response relationships, as well as their effects on one another. This method was used to identify the most important PFAS and their interactions in producing high allostatic load, while accounting for potential confounders and covariates. By considering the relationships between multiple variables simultaneously, BKMR modeling enabled a more comprehensive understanding of the complex interplay between the exposures and outcomes. Among the PFAS examined in this study, BKMR analysis revealed that higher serum concentrations of PFDE, PFNA, and PFUA tended to be associated with a higher allostatic load with these PFAS being identified as the most important in the relationship between PFAS and allostatic load. When allostatic load was modeled as a continuous variable the most critical PFAS were PFDE, PFOS, and PFNA. The fact that PFOS became more important in the analysis suggests that the nature of the data and the underlying relationship between allostatic load and PFAS is complex. To best capture nonlinear interactions between the variables, modeling allostatic load as continuous may paint a more accurate picture of the effects of multiple PFAS. This is so because allostatic load modeled as a continuous variable has more granularity to detect subtle changes, and as a categorical variable allostatic load only provides information on the category to which each observation belongs, and not on specific values within the category. Finally, as a continuous variable, allostatic load has increased sensitivity to changes or variations in data.
This study is the first of its kind to model PFAS with allostatic load using BKMR. That said, other studies may offer additional insight. Ding et al. conducted a study to examine the association between PFAS exposure and adipocytokine concentrations and found significant correlations between PFAS mixtures and the health-related outcomes of leptin and adiponectin in their models. Specifically, the PFAS chemicals PFDA, PFOS, PFNA, and PFHpS were identified as important contributors to these associations [35]. Tian et al. studied the link between PFAS exposure and folate levels using BKMR. They found that PFOA interacted with PFOS, PFNA, and PFDA to negatively affect RBC folate, with the negative impact increasing as PFOA levels rose. The mixture of five PFAS was inversely associated with RBC folate overall [36]. Other studies using BKMR analysis have shown adverse effects of multiple PFAS and liver injury [37], cardiovascular disease [38] and kidney dysfunction (glomerular filtration rate) [39] among others. These findings together speak to the plausibility of our hypothesis in the context of the results we obtained.
In our study, using BKMR's feature of being able to estimate hierarchical importance, PFNA was found to have a larger influence on allostatic load, followed by PFDE and PFUA when modeling allostatic load as a binary variable. Such insight can be critical for public health officials working in environments where the PFAS exposure is vast and contains various PFAS. By considering all exposures in a multivariable space and connecting the exposure with the health outcome via the Gaussian kernel function, the BKMR model allows for an excellent exploration of the relationship between PFAS and allostatic load.
The BKMR model's conclusions about the direct relationship between PFNA, PFDE, and PFOS with allostatic load when modeled as a continuous variable are interesting in that the most critical variables in order from high to low were PFOS, PFDE, and PFNA. Thus, when allostatic load is viewed as a continuous variable the traditional PFAS, PFOS, becomes the most important variable. This is especially insightful since BKMR modeling correct for the confounding effect, which is brought on by the strong correlations between PFAS when compared to traditional multiple linear regression models [2]. Specifically, BKMR deals with either high correlations among multiple PFAS or non-linear relationships and thus captures associations which may have been missed in traditional models.
PFAS compounds in combination with one another may interact and produce a cumulative effect that is either more or less than the total of the fractions of the individual PFAS because people are typically exposed to numerous PFAS at once from various sources [40]. The mechanisms by which the interactions occur between PFAS to produce allostatic load need to be further investigated to confirm our findings. Despite not performing mechanistic studies, we hypothesized that combined exposure to PFAS may increase allostatic load through promoting oxidative stress and diminishing antioxidant systems i.e., via GSH de-pletion, increasing blood pressure and cholesterol, and promoting an obesogenic state. This hypothesis is based on research that has occurred using animal models and in vitro [40][41][42].
The strengths of this study are as follows. The study used mixture analysis to evaluate the combined exposure of multiple PFAS on allostatic load. In addition, we used a nationally representative sample of adults. Finally, the study adjusted for the relevant confounders which lowered the bias that may have occurred. The cross-sectional study design, one of the study's drawbacks, makes it challenging to conclude whether a causal relationship between PFAS exposure and allostatic load exists. Finally, we did not employ the NHANES survey weights in BKMR models, comparable to other earlier studies using NHANES data, because weighting algorithms were not included in the R packages for BKMR.

Conclusions
In the current investigation, we found direct correlations between exposure to selected PFAS and allostatic load. The PFAS mixtures showed statistically significant overall direct correlations with allostatic load. Our study demonstrated that exposure to particular PFAS in combination might adversely promote the body's stress response.
Based on the results of the current study, there are significant associations between allostatic load and exposure to PFDE, PFNA, PFOS, and PFUA. Our findings suggest that an accumulation of PFAS exposure could have negative effects on the body's ability to cope with stress.