Bayesian Estimation of Human Population Toxicokinetics of PFOA, PFOS, PFHxS, and PFNA from Studies of Contaminated Drinking Water

The goal of this data collection was to find published PFAS concentrations in human serum and drinking water levels to predict those serum levels. We located human serum concentrations and corresponding data on DWCs by searching databases such as PubMed, ScienceDirect, and Google Scholar for English language journal articles and reports that presented human serum data for PFOA, PFOS, PFNA, or PFHxS levels. We also conducted a tree search on identified review articles and reports. We did not apply constraints on publication date, but the literature review to identify studies was concluded in 2019. We focused on populations exposed through contaminated drinking water, as defined by study authors. When serum data for additional PFAS were identified in studies where contamination of one PFAS was indicated prior to conducting the study, we did also collect these serum concentrations and include them if we could also identify corresponding drinking water information. In some cases, when PFAS levels in drinking water were reported as below the minimum reporting level (MRL) for the U.S. EPA’s third Unregulated Contaminant Monitoring Rule (UCMR3), but identified in serum of the study participants, we estimated the DWC as below the MRL as part of the analysis. We excluded studies focused on occupational cohorts, and where possible, we excluded individuals in community studies also known to have likely occupational exposures, such as employees of some chemical companies or firefighters, because our goal was to develop a tool that is predictive of community PFAS levels. For individual-level data, this information was only available for the Decatur, Alabama, data set, and so other data sets may include individuals with significant occupational exposure. Although several studies collected some information about occupation, only Bartell et al. (2010) and Emmett et al. (2006) specifically detailed how they excluded those with known or likely occupational exposure to the PFAS. All data were obtained from previous publications, and this analysis did not require institutional review board approval.

To be included in our model, a study needed to provide sufficient detail so that we could understand the timing of the PFAS exposure and map the serum levels to relevant DWCs for the study participants if the study did not explicitly give these concentrations. We reviewed the background, methods, results, figures, and tables of each study for reporting of community- and individual-level human serum data. Data extraction included recording the individual serum level if available, as well as the mean, geometric mean (GM), median, minimum, maximum, standard deviation (SD), sample size, dates of sampling, geographic location, and relevant demographic information, such as participant age and sex. When they were provided, we recorded dates of PFAS water contamination and remediation to better estimate participant exposure. If a study provided water consumption or finished water concentration data along with serum levels, we also collected this information. We aimed to develop as complete a picture as possible of DWCs over time prior to the serum measurements. When participant geographic information such as ZIP code, city, or water district service area was reported, we searched PFAS water concentration data to match them to serum data by location and dates. We obtained finished water concentration data for matching locations from journal articles, service districts’ water quality reports, and data collected under the U.S. EPA’s UCMR3. We used the Environmental Working Group’s (EWG’s) National Tap Water Database (U.S. EPA 2017; EWG 2019a) to identify additional sources of drinking water concentration data but verified data contained in the EWG database with the original sources. We did not include source water concentrations, and we included finished water data for all available time periods, including during participant sampling, before remediation, or after remediation, if available. Data recorded included water collection data, water type, arithmetic mean, GM, median, minimum, maximum, SD, and number of samples. The authors of several human serum studies reported the PFAS concentration that participants were exposed to in their drinking water, and in those cases, we recorded the PFAS concentration for the location directly from the serum-level studies.

To be included for analysis, studies with time-course serum data had to have available water concentration data during the time period between the first and last serum collection. For studies with only a single serum collection, water concentration data had to precede the serum collection date. If we could not reliably match study participants with reliable water concentrations, we excluded those serum levels. In addition, owing to limited data, no data sets focusing exclusively on children were included, and for a single study with both children and adults $>18\text{\hspace{0.17em}years}$ of age, only the adult values were used in estimation. We excluded individuals with a known potential for occupational PFAS exposure, those who could not be linked to specific water concentrations, those who reported filtering their water, and those with missing values for important variables such as weight and sex. We did include populations drinking bottled water in the testing data sets. We ultimately included data from nine sites for parameter estimation. Study data sets, citations, and location-specific water concentration handling are listed in Table 1 and, in further detail, in Tables S1 and S2. Table S1 includes details on why a location was excluded, which was typically due to a lack of sufficient information on water concentrations or serum concentrations.

In keeping with previous and current PBPK modeling by the U.S. EPA in adult (nonpregnant) humans (Lorber and Egeghy 2011; U.S. EPA 2016), we used a one-compartment toxicokinetic model for each PFAS. Moreover, the parameters of the one-compartment model are also included in human PBPK models that include pregnancy and lactation (MDH 2021; NJDEP 2021). Thus, this approach will enable direct use of our modeling results in a wide range of PFAS public health assessments. The parameters were defined as follows: $C_{bgd}$ is the background serum concentration (in micrograms per liter); $C_0$ is the initial serum concentration (in micrograms per liter) at $t=0$; $DW{I}_{BW}$ is the drinking water intake on a body weight-adjusted basis (e.g., as liters per kilogram per day); $k=\text{ln}\left(2\right)/T_{1/2}$ is the rate constant per day; and $V_d$ is the V_d (in liters per kilogram).

All of these parameters are assumed to be constant for the duration of the simulation, except for the DWC, $DWC\left(t\right)$ (in micrograms per liter), which may vary over time. The ordinary differential equation (ODE) for the concentration of PFAS in the body is as follows:

As shown in Figure 1, this model is equivalent to the usual one-compartment model used for pharmaceuticals, but with the usual bolus dose (typically modeled as either an initial condition or an exponential input) replaced by a time-varying input that includes both a drinking water component (first term) and a background exposure component. Moreover, the background dose rate (e.g., $D_{bgd}$, in micrograms per kilogram per day) has been reparameterized in terms of the background serum concentration, $C_{bgd}=D_{bgd}/\left(k{V}_d\right)$. This reparameterization improves model fitting by reducing parameter correlations. In addition, because clearance is used as the basis for both chemical-specific interspecies extrapolation and interindividual variability factors (WHO et al. 2005; U.S. EPA 2014), CL (in liters per kilogram per day) $=k\times V_d$ was also calculated.

For a constant water concentration, Equation 2 has an analytic solution for the serum concentration as a function of time:where $C_0$ is the concentration at time 0, and $C_{ss}$, the steady-state serum concentration (due to drinking water alone, without background), is given by Equation 3:

The first term in Equation 2 is simply the background concentration, the second term describes the time-dependence from time 0 to background in the absence of water contamination, and the third term describes the time-dependence of the transition to steady state. If data are available only at a single time point, then steady state is assumed, and DWCs are assumed to be constant:

In some cases, the data at a single time point involve an intervention in which the source of contamination was removed at a certain time interval, $\mathrm{\Delta }t$, in the past. This simplification of intervention is necessary when no water concentrations for after intervention are available. In this case, the serum concentration takes the following form:

We employed a Bayesian population approach to estimate model parameters. This approach involves two statistical models: one for population variability and one for the likelihood of the observed data (serum concentrations) given a set of model parameters (such as T_½ and V_d). With respect to the population variability model, each individual $i$ is assumed to have their own (unknown) model parameters ${\mathrm{\theta }}_i$ and serum concentration data $C_i$ at time $t_j$. The individual parameters (log-transformed) are drawn from a normal population distribution:

Here, ${\mathrm{\alpha }}_{\mathrm{\theta },i}$ represents the $z$-score of the individual in the population for parameter $\mathrm{\theta }$. For some parameters, the population mean, ${\mathrm{\mu }}_{\mathrm{\theta }}$, and SD, ${\mathrm{\Sigma }}_{\mathrm{\theta }}$, may be fixed (e.g., background, drinking water intake), whereas for others, ${\mathrm{\mu }}_{\mathrm{\theta }}$ and ${\mathrm{\Sigma }}_{\mathrm{\theta }}$ may have distributions that are due to uncertainty (e.g., T_½, V_d).

With respect to the statistical model for the likelihood of the observed data, we assume lognormal errors with an additional parameter to specify the SD of the error distribution, so the likelihood function for an observation, $C_{i,obs}\left(t_j\right)$, in individual $i$ at time point $t_j$ would be defined as follows:where $C_{pred}$ is the model prediction for serum concentration as a function of time as described above. Different values for the error of the geometric standard deviation ($GS{D}_{err})$ are used for individual time-course data, individual steady-state data, and summary-level steady-state data.

For individual data, there are two common situations. In the first, time-course data are available for the individual. In this case, the model sets $t=0$ at the initial sampling time point, where $C\left(0\right)=C_0$. We then evaluate subsequent time points at times after the initial sample point using Equation 1 for time-varying water concentrations and Equation 2 for constant water concentrations. This approach avoids the uncertainty regarding what the background and water concentrations were before the first sampling point but requires data on the water concentrations between the first and last serum sampling point. In the second situation, serum data are available only at a single time point, but historical water concentration data are available. In this case, we use Equation 5, where steady state is assumed until the point of intervention (if any), after which the water concentration level is assumed to be negligible. This simplification allowed us to include data from study locations with no postintervention water concentration data.

If only summary data are available, then the only statistic that the model can directly predict at the population level is the population arithmetic mean because of the presence of the background term. For instance, at steady state, taking the mean of Equation 4, we have the following:The $C_{bgd}$ term has a population mean as follows:

Because each of the parameters—$DW{I}_{BW}$, $DWC$, $k$, and $V_d$—that make up $C_{ss}$ are assumed to be lognormally distributed in the population, their product (and quotient) is also lognormally distributed. Thus, the population mean of $C_{ss}$ is given as follows:where the population GM $\left(e^{\mu }\right)$ and geometric standard deviation (GSD) $\left(e^{\Sigma }\right)$ of $C_{ss}$ are related to those of the other parameters as follows:

In the case of an intervention at time $\mathrm{\Delta }t$ prior to the serum measurement, taking the population mean of Equation 5 gives the following:Because the term $e^{-k\mathrm{\Delta }t}$ cannot be expressed exactly in a formula, we applied an approximation. When appearing in the exponential, $k$ is assumed to be normally, rather than lognormally, distributed, with the mean and variance matched to the actual distribution of $k$. This implies that $e^{-k\mathrm{\Delta }t}$ is lognormally distributed, with GM $e^{-\langle k\rangle \mathrm{\Delta }t}$ and GSD $e^{\mathrm{\Delta }t\sqrt{VAR\left(k\right)}}$. Then, we make an additional approximation of independence so that each term above consists of a product of independent lognormal distributions. Thus, the product $Y=C_{ss}{e}^{kt}$ has a population mean as follows:

A similar approach can be used for summary data consisting of multiple time points as long as the DWC is constant. Taking the arithmetic mean of Equation 2 and separating the terms gives the following:

For the terms $\langle C_0{e}^{-kt}\rangle$ and $\langle C_{bgd}{e}^{-kt}\rangle$, we can use the same approximation from Equations 14–16, replacing the ${\mu }_{C_{ss}}$ and ${\text{Σ}}_{C_{ss}}^2$ with the corresponding values for $C_0$ and $C_{bgd}$. We checked the approximations in Equations 14–17 with simulations and found them to have errors mostly of $<10\%$ for $\mathrm{\Delta }t$ up to about twice the T_½.

Prior distributions for each model parameter are summarized in Table 2 (for additional detail, see Table S3). For the background serum concentrations, $C_{bgd}$, NHANES data collected closest to the year(s) during which the study data were collected were used as the prior central estimate (Table S3). However, 80% of the appropriate NHANES level was used as $C_{bgd}$, corresponding to a relative source contribution from drinking water of 20%, so as to avoid double-counting background exposures with exposures due to the measured DWCs. This assumption is consistent with U.S. EPA drinking water guidance, which as a default assumes that 20% of exposure to a contaminant is from drinking water in the absence of chemical-specific data (U.S. EPA 2018). For time-course data, the prior central estimate for the initial concentration $C_0$ was set to the reported initial value for each individual. For drinking water intake, the population distribution was fixed based on community water source intake data from the Food Commodity Intake Database “What We Eat in America” (FoodRisk 2020) for the years 2005–2010 for those 16–81 years of age (consumers only). For residual error, each data type (individual time-course, individual steady-state, and summary steady-state) involves different assumptions, so each is assumed to have a different residual error.

For elimination rates, we assigned informative prior distributions for each PFAS. We reviewed the available body of literature on T_½ for the four PFAS of interest. We focused on those studies that estimated T_½ from consuming contaminated drinking water. We focused exclusively on nonoccupational populations exposed through contaminated drinking water given that this is a common exposure pathway of relevance to public health and eliminates any potential differences in toxicokinetics associated with high occupational exposure levels. Therefore, we excluded occupationally exposed cohorts without drinking water measurements because ongoing exposures in the community can confound T_½ estimates (e.g., Costa et al. 2009; Olsen et al. 2007). We also excluded studies that estimated elimination from urine alone because total elimination may include other pathways for some PFAS (e.g., Harada et al. 2005; Zhang et al. 2013) and sufficient data were available for estimation without including those studies. However, the results of these studies were considered in setting the bounds on the prior distributions for T_½.

For PFOA, prior distributions were based on the estimates from Bartell et al. (2010) and Seals et al. (2011), the smallest and largest T_½ values identified, respectively. The population mean was centered on the mean log elimination rate [corresponding to a median T_½ of 4.6 y with uncertainty $\text{GSD}=1.5$, so the 95% confidence interval (CI) of 2.1 to 10.2 covers both studies]. Population variation was based on the reported 95% interval of individual T_½ from Bartell et al. (2010) $\left(\log \text{\hspace{0.5em}SD}=0.286\right)$. Assuming the precision ${\mathrm{\Sigma }}^{-2}$ has a prior uncertainty coefficient of variation $\left(\text{CV}\right)={\mathrm{\alpha }}^{-\mathrm{½}}=33\%$, this implies a shape parameter $\mathrm{\alpha }=9$. The rate parameter $\mathrm{\beta }$ is then derived by matching the prior $\text{mean}=\mathrm{\alpha }/\mathrm{\beta }={0.286}^{-2}\approx 12$. For PFOS, this is based on the range of several studies, centered on a GM of 4.8 from Olsen et al. (2007) and the same GSD as PFOA. Based on this range, we dropped a study that found a T_½ of 27 y for PFOS in adults $>50\text{\hspace{0.17em}years}$ of age (Zhang et al. 2013) given that it is an outlier compared with the other studies. For population variation, because PFOA has the most data, we used the posterior from PFOA as the prior for PFOS and the other PFAS. For PFNA and PFHxS, the prior central estimates for the T_½ were 4.3 y (from Zhang et al. 2013) and 5.3 y (from Li et al. 2018), respectively. Limited T_½ data were available for PFNA; therefore, we used the data from Zhang et al. (2013) to center our prior distribution even though these data were based on blood–urine pairs. We chose to use the higher estimate (4.3 in all males and females $>50\text{\hspace{0.17em}years}$ of age) vs. the lower (2.5 y in younger females) as our central estimate for the prior distribution given that the elimination rate is hypothesized to be longer for the longer-chained PFAS (PFNA is a 9-carbon chain PFAS) (Graber et al. 2019). For PFHxS, we chose to center our estimate on the study by Li et al. (2018) that the MDH used in their PBPK model for PFHxS. T_½ for PFHxS also ranged widely in the literature, from 4.7 y (females only) to 15.5 y (Li et al. 2018; Worley et al. 2017a), when excluding the urine-based outlier of 35 y from Zhang et al. (2013).

We based V_d priors on the same literature search as was conducted for T_½. The prior distribution for the population GM was centered on $0.17\mathrm{L}/\mathrm{kg}$ for PFOA and $0.23\mathrm{L}/\mathrm{kg}$ for PFOS, based on Thompson et al. (2010) with uncertainty $\text{GSD}=1.3$ for both PFOA and PFOS. The prior distribution for population variation was based on a weakly informative half-normal prior (per a recommendation by Gelman 2006) with mean of 0.16. For PFHxS, the prior central estimate was $0.25\mathrm{L}/\mathrm{kg}$, based on Sundström et al. (2012), which found a V_d range of $0.2\text{\hspace{0.17em}to\hspace{0.17em}}0.3\mathrm{L}/\mathrm{kg}$ informed by rat, mice, and monkey data, and Koponen et al. (2018), which assumed the same V_d for PFHxS as for PFOS. Owing to a lack of data, the prior central estimate for PFNA was assumed to be the same as for PFOA ($0.17\mathrm{L}/\mathrm{kg}$), consistent with the assumption made in Koponen et al. (2018), which was informed by rodent data.

For some studies, the DWCs are reported to be “below the minimum reporting level.” In this case, a uniform distribution uncertainty between zero and the reporting level is used as a prior for the actual DWC.

Data were separated into “training” and “testing” data sets, with only training data used for calibration (i.e., part of the likelihood function). For studies with individual data, training data consisted of half the individuals of each sex, randomly selected, with the remaining individuals treated as testing data sets. For studies with summary data only, half of the studies were selected for training and half for testing. In both cases, if an odd number of individuals or studies were present, the additional case was randomly assigned to either testing or training.

We implemented the model in the open-source software MCSim (version 6.1.0; GNU MCSim; https://www.gnu.org/software/mcsim/), and all analyses were performed in R within RStudio (R Development Core Team; version 3.6.1; Rstudio Team). Plotting and summarization were performed with tidyverse packages (Wickham et al. 2019). Four independent Markov chain Monte Carlo (MCMC) chains were run for each PFAS. We assessed convergence using the potential scale reduction factor ($\widehat{R}$), which approaches 1.0 with convergence, and for which a value of $\le 1.2$ is proposed as acceptable (Gelman et al. 2013). As mentioned above, PFOA was run first, because it has the most individual data. The posterior distribution for the population variance of the elimination rate for PFOA was used as the prior distribution for this parameter for the other PFAS (i.e., replacing the inverse gamma distribution prior used for PFOA). We fit the PFOA posterior for this variance parameter to a lognormal distribution. R codes for all modeling and analyses are included in supplementary materials and available at https://github.com/wachiuphd/2022-Bayes-PFAS-PK.

We performed multiple local sensitivity analyses with respect to our modeling assumptions and approach to validate the parameter estimates produced by the model. These included a) changing the relative source contribution from drinking water in background serum concentration from 20% to either 0% or 80%, b) shifting the prior distributions for V_d to 20% above or below the primary estimate, and c) evaluating the effect of using alternative training vs. testing data sets. Specifically, we analyzed the effect of using only individual-level serum concentrations, using only population-level aggregates, and switching the training and testing data sets from the primary analysis. Finally, we performed an analysis focused only on parameter estimation, where all data were used for training. In all cases, we reran the entire model calibration and evaluation process using the alternative parameter values, prior distributions, or data sets.

Comprehensive details on the identified populations and studies are provided in Table 1 and in Table S1. Briefly, for the Decatur population ($n=37$ after excluding eight individuals owing to occupational history), the Agency for Toxic Substances and Disease Registry (ATSDR) supplied the individual serum data. For the Arnsberg population, the study by Hölzer (2008) provided individual serum PFOA levels graphically for 2006 and 2007 ($n=151$). We derived and estimated those values from the line graph provided in the study to use them as individual serum data inputs in the model at two time periods. For the Minnesota population, the study by Johnson et al. (2017) provided individual-level serum data for PFOS and PFOA correlated with individual drinking water exposure ($n=98$). These were also derived from digitizing figures with WebPlotDigitizer (version 4.2) because tabular data were not provided. We obtained all other data from text or tables from study publications. Individual water consumption information was not available in any of the identified studies.

A number of identified populations were excluded because the publications were missing critical information. We excluded the North Wales, Pennsylvania, cohort because the water data from before the serum measurements were not available and because some of the water was purchased from a neighboring town (North Wales Water Authority 2018). We excluded the Ronneby cohort because serum concentrations were not reported at individual time points (Li et al. 2018). The Uppsala County cohort was excluded because it contained only a cumulative estimate of months of exposure for children. For the Arnsberg cohort, we included only data for adults, given that it was the only study that separated children’s serum levels from adults’.

The study populations we used for modeling, along with summary statistics of the serum and water concentrations for each PFAS, are summarized in Table 1. In some cases, there was only a single water or serum level for an individual or a population. In other cases, there were multiple water levels preceding the serum level, or ideally, more than one serum level with corresponding water concentrations over time at the individual level.

Parameter estimates in all four models achieved excellent convergence ($\widehat{R}<1.05$) with a reasonable number of MCMC iterations (4 chains; 20,000 iterations per chain). The resulting overall model fits comparing posterior median predictions and data are shown in Figure 2. For training data (Figure 2, left panels), the model was able to match the data very tightly. The GSD_err for individual time-course data was $\sim 1.1$ for all four PFAS, indicating that the model fit had a residual error (difference between predictions and training data) of only $\sim 10\%$ when time-course data were available. In cases where data are only available at a single time point per individual, such as for the Minnesota data for PFOA and PFOS, the GSD_err was $\sim 1.5$, indicating $\sim 50\%$ residual error. A larger error such as this one is expected because of the need to approximate steady state in these cases. Moreover, the Minnesota data appeared to tend toward overprediction, which is consistent with the use of a steady-state approximation (i.e., in reality, steady state would not have been reached, so the actual concentration would tend to be lower). Summary time-course data, available only for PFOA, had a GSD_err of $\sim 1.2$, indicating $\sim 20\%$ error. For summary data at a single time point, however, the residual errors were quite a bit larger: $\sim 2$-fold for PFOA and PFOS, 2.3-fold for PFNA, and 2.8-fold for PFHxS. Remarkably, however, the residual error in all cases was $<3$-fold. These results show the importance of individual data for accurately estimating PFAS elimination.

The comparisons with the testing data (Figure 2, right panels) showed a similar trend. For individual data, the model performed best when we used data sets containing individual time-course data, with residual errors well within 50%. Individual data with a single time point and summary steady-state data performed worse, although residual errors were generally within 3-fold. Figure 3 shows the posterior distributions of predictions for the Decatur time-course data, which were available for all four PFAS. For both training and testing data sets, the data were well within the CIs of the posterior predictions. Results were similarly accurate for the individual data sets across all PFAS (Figures S1–S4). However, for summary data, the predictions were less accurate for both the training and testing data sets.

Posterior distributions for the main toxicokinetic parameters of T_½ and V_d, as well as the derived parameter of clearance, are shown in Table 3. These reflect updating of the prior distributions after consideration of the likelihood of the data. In addition to estimates for the population GM and population GSD, the model makes a prediction for a “random individual,” which combines uncertainty and variability. The random individual is relevant to the general public because for any individual, there are two sources of uncertainty: a) the uncertainty in the population distribution (GM and SD), and b) uncertainty in where one is located on the population variability distribution (e.g., does an individual have a higher or lower than typical T_½ for a PFAS?).

Comparisons of prior and posterior distributions for the T_½ are shown in Figure 4. For all four PFAS, the posteriors for the population GMs were noticeably shifted and narrower than the prior distributions. This indicates that the data were informative relative to the prior. For PFOA, the posterior for the population GSD was also shifted and indicated more population variation in the T_½ than under prior assumptions. However, for the other PFAS, the priors and posteriors for the population variation were similar because there were insufficient numbers of individuals to inform this parameter; in other words, the data were not informative relative to the prior. Nonetheless, this indicates that the data are consistent with the amount of population variation in T_½ observed with PFOA. T_½ estimates (95% CIs) for the population GM were 3.14 (2.69, 3.73) y for PFOA, 3.36 (2.52, 4.42) y for PFOS, 2.35 (1.65, 3.16) y for PFNA, and 8.30 (5.38, 13.5) y for PFHxS. For PFHxS, there was a noticeable shift from the prior to the posterior, with the central estimate moving from 5.29 to 8.30 y. Although the uncertainty ranges in Figure 4 look similar on the natural scale, the 95% CI range on the log-scale shrank substantially, from 4.9- to 2.5-fold.

Comparisons of prior and posterior distributions for the V_d are shown in Figure 5. V_d posterior estimates ranged from $0.19\text{\hspace{0.17em}to\hspace{0.17em}}0.428\mathrm{L}/\mathrm{kg}$ across the four PFAS. In general, the posterior estimates for the V_d across the four PFAS are shifted to larger values as compared with the prior distributions, indicating that the data were informative for the population mean of V_d. A larger V_d indicates that for a given intake, the ratio between DWCs and serum concentrations would tend to be larger as compared with previous studies; this result is likely due to accounting for background (nondrinking water sources). There are insufficient data, however, to substantially inform the population variability in the V_d given that the prior and posterior distributions are similar.

With respect to clearance, which is the product of the rate coefficient k and the V_d, posterior medians ranged from $0.025\text{\hspace{0.17em}to\hspace{0.17em}}0.095\mathrm{L}/\mathrm{kg}\text{\hspace{0.17em}per\hspace{0.17em}year}$, with $\text{PFHxS}<\text{\hspace{0.17em}PFNA}<\text{\hspace{0.17em}PFOS}<\text{\hspace{0.17em}PFOA}$. These values tended to be on the higher end of the prior distributions for PFOA, PFOS, and PFNA, and toward the lower end for PFHxS. We also found that posterior distributions for the components of clearance, T_½, and V_d were uncorrelated (all $R^2<0.05$). Posteriors for the population variation were slightly larger than those for T_½, reflecting the relatively smaller contribution of V_d to interindividual variation.

To check for the possibility of systematic biases, we examined the posterior distributions for individual k and V_d parameters and compared them across study cohorts. For instance, if individual posteriors for one location were systematically different from individual posteriors for another location, then that would suggest errors in the model or parameters. As shown in Figure S5, there is no discernable difference across cohorts from different cities/locations, and all individual parameter posterior samples had $z$-scores that were statistically consistent with the expected standard normal distribution.

The posterior estimates for parameters were generally stable across all sensitivity analyses, and in all cases, there was substantial overlap between the 95% CI from each sensitivity analysis and that from the primary analysis (Table S4, Figures S6–S9). Changing relative source contribution from drinking water and nondrinking water sources had little effect on T_½ or V_d estimates. Swapping the data sets used for testing with those used for training produced nearly identical results to the primary analysis. Increasing and decreasing the priors for V_d resulted in nearly identical posterior estimates for all parameters. The largest differences from the primary analysis occurred when removing either the individual- or the population-level data for training. The posteriors for the GSD for variability in T_½ were much lower when we used only population summary data, as would be expected because it is more challenging to estimate variability when only summary data are available.