A Probabilistic Approach to Evaluate the Risk of Decreased Total Triiodothyronine Hormone Levels following Chronic Exposure to PFOS and PFHxS via Contaminated Drinking Water.

BACKGROUND: Extensive exposure to per- and polyfluoroalkyl substances (PFAS) have been observed in many countries. Current deterministic frameworks for risk assessment lack the ability to predict the likelihood of effects and to assess uncertainty. When exposure exceeds tolerable intake levels, these shortcomings hamper risk management and communication.
OBJECTIVE: The integrated probabilistic risk assessment (IPRA) combines dose-response and exposure data to estimate the likelihood of adverse effects. We evaluated the usefulness of the IPRA for risk characterization related to decreased levels of total triiodothyronine (T3) in humans following a real case of high exposure to PFAS via drinking water.
METHODS: PFAS exposure was defined as serum levels from residents of a contaminated area in Ronneby, Sweden. Median levels were 270 ng/mL [perfluorooctane sulfonic acid (PFOS)] and 229 ng/mL [perfluorohexane sulfonic acid (PFHxS)] for individuals who resided in Ronneby 1 y before the exposure termination. This data was integrated with data from a subchronic toxicity study in monkeys exposed daily to PFOS. Benchmark dose modeling was employed to describe separate dose-effect relationship for males and females, and extrapolation factor distributions were used to estimate the corresponding human benchmark dose. The critical effect level was defined as a 10% decrease in total T3.
RESULTS: The median probability of critical exposure, following a combined exposure to PFOS and PFHxS, was estimated to be [2.1% (90% CI: 0.4%-13.1%)]. Gender-based analysis showed that this risk was almost entirely distributed among women, namely [3.9% (90% CI: 0.8%-21.6%)]. DISCUSSION: The IPRA was compared with the traditional deterministic Margin of Exposure (MoE) approach. We conclude that probabilistic risk characterization represents an important step forward in the ability to adequately analyze group-specific health risks. Moreover, quantifying the sources of uncertainty is desirable, as it improves the awareness among stakeholders and will guide future efforts to improve accuracy. https://doi.org/10.1289/EHP6654.

Introduction

Per- and polyfluoroalkyl substances (PFAS) is a class of chemicals used in a multitude of applications (Buck et al. 2011; Lau 2015). Although they have useful technical properties, some of these chemicals have been found worldwide in human, animal, and environmental samples (Banzhaf et al. 2017; Jian et al. 2017). Worldwide, the presence of high levels of these chemicals in drinking water, groundwater, and ecosystems has been associated with industrial production sites or the usage of aqueous film-forming firefighting foams (Banzhaf et al. 2017; Filipovic et al. 2015; Thalheimer et al. 2017). In a situation where exposure is clearly above the level of concern, for example, by exceeding tolerable intake levels, there is a need to characterize the risk in terms of probability for adverse health effects. Probabilistic risk assessment has been proposed to provide decision makers and other stakeholders with information about the magnitude of health risks and, at the same time, present uncertainties in a transparent and precise manner (WHO/ILO/UNEP 2014). This approach has so far not been used to evaluate exposure scenarios with PFAS, and it remains to be elucidated how probabilistic risk estimates differ from the traditional deterministic approaches that calculate a discrete Margin of Exposure (MoE) for the exposed population.

The main aim of the present study was to perform an integrated probabilistic risk assessment (IPRA) analysis and to discuss its usefulness as a potential tool for risk assessment in comparison with the traditionally used deterministic MoE approach. The IPRA method, first described by van der Voet and Slob (2007), combines dose–response data with exposure data in a probabilistic manner in order to generate a distribution of individual margins of exposure (IMoE) and to estimate a probability of critical exposure (PoCE), that is, the likelihood of having an exposure that exceeds a predefined effect level. In addition, the IPRA approach takes data variation into account, which can be divided by variability and uncertainty sources (e.g., intraspecies variability and duration extrapolation uncertainty) and quantifies their contribution to the final estimate (van der Voet and Slob 2007). This uncertainty estimation is a major advantage over traditional deterministic risk assessment strategies such as the MoE approach, which results in a single estimate where the impact of uncertainty remains unknown (van der Voet and Slob 2007). The IPRA approach expresses the estimated risk though a confidence interval (CI) that describes the percentile of the population at risk of having an effect due to the chemical exposure (van der Voet and Slob 2007). Previously, the approach has been used for both human and environmental risk assessment purposes for agents as diverse as nanoparticles (Jacobs et al. 2016), cadmium, mycotoxins, pesticides, and acrylamide (Bokkers et al. 2009). The IPRA approach is most often used in cases of general population exposure and, to the best of our knowledge, it has not been applied to cases where the population has had an extensive exposure to soil and water contaminants.

In order to investigate the usefulness of the probabilistic approach for risk assessment purposes, we performed an analysis focusing on PFAS in drinking water, an emerging issue in many industrialized countries. PFAS represent a large class of chemicals known to be bioaccumulative, persistent and toxic (Jian et al. 2017). Among other effects, PFAS have been described as hepato- and immunotoxic and as developmental toxicants (ATSDR 2018; Wang et al. 2017). A growing body of literature points toward the understanding that thyroid hormone levels can be affected following exposure to perfluorooctane sulfonic acid (PFOS) and perfluorohexane sulfonic acid (PFHxS). Disruption of the thyroid hormone system has been reported in both epidemiological (Dallaire et al. 2009; Knox et al. 2011; Wen et al. 2013) and animal studies (Chang et al. 2008; Curran et al. 2008; Luebker et al. 2005; Martin et al. 2007; Seacat et al. 2002; Thibodeaux et al. 2003). However, not all epidemiological studies have found a negative correlation with PFOS levels in serum (Olsen et al. 2003), nor all animal studies (Chang et al. 2017). Nevertheless, PFOS has been identified as a thyroid hormone disruptor (ATSDR 2018; Coperchini et al. 2017).

The mechanistic aspects of PFOS-induced effects on the total triiodothyronine () hormone are numerous and complex. In epidemiological studies, effects on the thyroid hormone system might be difficult to observe due to a large interindividual variation. Moreover, the hypothalamus–pituitary–thyroid (HPT) axis is regulated by a tight feedback mechanism. Thus, xenobiotics leading to changes in circulating and thyroxine () levels can be considered disruptive to the HPT homeostasis, but that is not necessarily the cause of disease given that this can be compensated for by an increase in thyroid-stimulating hormone (TSH) secretion. It is plausible that PFOS exerts its thyroid-disrupting effects by displacing circulating and hormones from their binding proteins in the bloodstream, namely albumin, thyroxine-binding globulin, and transthyretin (Chen and Guo 2009; Hebert and MacManus-Spencer 2010; Ren et al. 2015, 2016).

For the IPRA analysis performed in the present study, we used serum levels of PFAS measured in habitants of the Swedish municipality of Ronneby. In December 2013, high levels of PFAS were detected in the outgoing water of the Brantafors waterworks, which provides drinking water to approximately one-third of Ronneby’s 28,000 inhabitants (Li et al. 2018; Livsmedelsverket 2013a). The Brantafors waterworks are located about from the firefighting training location of Ronneby’s airfield where PFAS-containing firefighting foams were estimated to have been used since the mid-1980s (Li et al. 2018). However, the exact composition of the firefighting foams that were used, the average annual volume usage, and the frequency of the firefighting training sessions in the Ronneby airfield are unknown. The two main contaminants identified in drinking water delivered from Brantafors were PFOS and PFHxS (Li et al. 2018; Livsmedelsverket 2013a). Based on data from this population, the average half-lives of PFOS and PFHxS were estimated to be 3.4 y and 5.3 y, respectively (Li et al. 2018).

Given the persistency and potential toxic effects of PFOS and PFHxS, both compounds are currently objects of concern to various international agencies (ECHA 2017; U.S. EPA 2016). PFOS has also been listed under Annex B (Restriction) of the Stockholm Convention on Persistent Organic Pollutants (UNEP 2017a). Chemicals proposed for listing under the convention and currently under review include PFHxS and perfluorooctanoic acid (PFOA) (UNEP 2017b).

Dose–response data were obtained from an animal study where cynomolgus monkeys were exposed daily to PFOS (Seacat et al. 2002). Among other effects, the exposure caused a significant decrease of total and free and increased levels of TSH (Seacat et al. 2002) (see Figure S1). In 2008, the European Food Safety Authority (EFSA) selected the paper by Seacat et al. (2002) as a critical study to establish a tolerable daily intake (TDI) value for PFOS (EFSA 2008).

Traditionally, risk is characterized by estimating the MoE, calculated by dividing an experimentally derived point-of-departure (POD), such as a no-observed-adverse-effect-level (NOAEL) or benchmark dose lower bound (BMDL), and an exposure metric from a typical or worst-case scenario. If the exposure is 100 times lower than the POD, there is generally no concern for health effects. However, this deterministic MoE approach does not express the contribution of uncertainty or variability to the estimate, nor does it inform how many individuals are likely to be affected if the exposure exceeds the reference value. For PFAS, there is currently no consensus regarding the threshold for a safe level of exposure and multiple guideline values, for both PFOS and PFHxS, have been established by different agencies (Table 1).

Table 1

A noncomprehensive compilation of guideline values for chronic exposure to PFOS and PFHxS expressed in terms of body weights.

Authority	Year	PFOS limit value	PFHxS limit value	Effect	Reference	MoE
ATSDR	2018	—	MRL=20 ng PFHxS/kg BW per dayCombined UF of 300	Rat, thyroid follicular cell damage	ATSDR 2018	0.59
Swedish National Food Agency (Livsmedelsverket)	2013	—	TDI=5,000 ng PFHxS/kg BW per dayCombined UF of 200	Rat, liver effects	Livsmedelsverket 2013b	147
ATSDR	2018	MRL=2 ng PFOS/kg BW per day Combined UF of 300	—	Rat, delayed eye opening and decreased pup weight	ATSDR 2018	0.018
EFSA	2020a	TWI=8 ng/kg BW per week for the sum of four PFAS (PFOS, PFHxS, PFNA and PFOA) ∼1.1 ng/kg BW per day		Epidemiological studies, lower antibody titers	Unpublished worka	0.008
EFSA	2018	TWI=13 ng PFOS/kg BW per week∼1.86 ng PFOS/kg BW per dayNo UFs used	—	Epidemiological studies, increase of serum cholesterol	EFSA 2018	0.016
EFSA	2008	TDI=150 ng PFOS/kg BW per dayCombined UF of 200	—	Monkeys, lipids and thyroid metabolism	EFSA 2008	1.3
U.S. EPA	2016	RfD=20 ng PFOS/kg BW per dayCombined UF of 30	—	Rat, decreased neonatal body weight	U.S. EPA 2016	0.18

Note: —, no data; ATSDR, Agency for Toxic Substances and Disease Registry; BW, body weight; EFSA, European Food Safety Authority; EPA, Environmental Protection Agency; MoE, Margin of Exposure; MRL, minimal risk level; PFHxS, perfluorohexane sulfonic acid; PFNA, perfluorononanoic acid; PFOA, perfluorooctanoic acid; PFOS, perfluorooctane sulfonic acid; Rfd, reference dose; TWI, tolerable weekly intake; UF, uncertainty factor.

This TWI has not been adopted at the time of publication of this article.

In the present study, we used the case of the PFAS-contaminated drinking water in Ronneby and the potential effects on thyroid hormone levels to evaluate the usefulness of the IPRA method in comparison with the commonly used deterministic approach MoE.

Methods

The IPRA analysis was performed in accordance with the approach described by van der Voet and Slob (2007). This full probabilistic analysis used Monte Carlo simulations, which repeatedly sampled randomly from the input data, to emulate possible exposures and dose–effect scenarios (van der Voet and Slob 2007; WHO/ILO/UNEP 2014). The MoE can be defined as the ratio of a NOAEL or BMDL-value and a point estimate of the exposure (WHO/ILO/UNEP 2014). Analogously, it has been suggested that the IMoE in an integrated probabilistic risk assessment can be defined as follows (van der Voet and Slob 2007): where IEXP is the individual exposure distribution, and IBMD is the probabilistically derived individual benchmark dose distribution. In the present analysis, each IMoE distribution described a scenario of possible combinations that departed from the original measurements in the exposed population. This process was repeated 10,000 times so that the IEXP, IBMD, and IMoE distributions would reflect as many scenarios as possible. The analysis was performed in R (version 3.4.2; R Development Core Team). An overview of the IPRA process is given in Figure 1.

Figure 1.

A schematic overview of the integrated probabilistic risk assessment approach, describing the simulations performed for residents that lived in Ronneby, Sweden, for at least 1 y (). The same process was repeated for residents living for 10 y () and 29 y (). The IEXP distribution size was always 1,000 greater than the size of the respective population from which the draws were originally made. The IBMD distribution size matched the number of IEXPs, keeping the gender proportions. Note: BMD; benchmark dose; EF, extrapolation factor; , GSD for interindividual human variability; IBMD, probabilistically derived individual benchmark dose distribution; IEXP, individual exposure distribution; IMoE, individual margins of exposure; PoCE, probability of critical exposure; Prob, probability; TD, toxicodynamics.

Exposure Assessment

The serum concentrations were determined at the Division of Occupational and Environmental Medicine, Department of Laboratory Medicine, at Lund University, Lund, Sweden. The quantification was performed using liquid chromatography and tandem mass spectrometry after extraction by protein precipitation (Li et al. 2018). The PFAS analyses are part of an interlaboratory control progran coordinated by the University of Erlangen-Nuremberg, Germany). For the present study, we used results from serum samples from Ronneby’s municipality residents, obtained between June 2014 and December 2015, from 1,845 individuals living in households provided with PFAS-contaminated drinking water. Given the high PFAS concentrations found in the drinking water, this exposure route was assumed to outweigh all others when contributing to the PFAS body burden (Table 2; see also Table SI) (Li et al. 2018; Livsmedelsverket 2013b). Individuals years of age were excluded from the analysis due to PFAS exposure during fetal life and via breastfeeding (Mondal et al. 2014). The residence history was collected from all participants. All individuals included in the analysis lived continuously in the areas provided with PFAS-contaminated drinking water for at least 1 y before exposure was terminated in December 2013. To explore the distribution of risk in the population, we stratified the group by gender and duration of residency, that is, the group was divided into those who lived continuously in the households provided with PFAS-contaminated drinking water during the last exposure year (2013, ), the last 10 y (2004–2013, ), or the last 29 y (1985–2013, ) before exposure was terminated in December 2013. Population descriptors for the three groups are presented in Table 3. The 5th, 50th, and 95th percentiles of PFOS and PFHxS serum levels measured in this population are presented in Table 4. Calculations were based on exposure to PFOS only and to a combined exposure to PFOS and PFHxS, assuming an equipotent toxicity on a molar basis. PFOS and PFHxS molecular weights are 500.13 and , respectively.

Table 2

Measured levels of PFAS in drinking water from the Brantafors waterworks in December 2013.

Chemical name	CASN	Concentration (ng/L)a
Perfluorinated sulfonate acids (PFSA)
Perfluorooctane sulfonic acid (PFOS)	1763-23-1	4,000
Perfluoroheptane sulfonic acid (PFHpS)	375-92-8	67
Perfluorohexane sulfonic acid (PFHxS)	355-46-4	1,200
Perfluorobutane sulfonic acid (PFBS)	75-22-4	140
Perfluorinated carboxylic acids (PFCAs)
Perfluordodecanoic acid (PFDoDA)	307-55-1	<10
Perfluroundecanoic acid (PFUnDA)	2058-94-8	<10
Perfluorodecanoic acid (PFDA)	335-76-2	<10
Perfluorononanoic acid (PFNA)	375-95-1	1.2
Perfluorooctanoic acid (PFOA)	335-67-1	130
Perfluoroheptanoic acid (PFHpA)	375-85-9	40
Perfluorohexanoic acid (PFHxA)	307-24-4	340
Perfluoropentanoic acid (PFPeA)	375-85-9	52
∑PFAS		5,970.2

Note: CASN, Chemical Abstracts Services number; PFAS, per- and polyfluoroalkyl substances; , sum of PFAS.

Data from Livsmedelsverket (2013a).

Table 3

Sampled Ronneby, Sweden, population description, gender, and age. Duration exposure regarding those living continuously in the areas provided with PFAS-contaminated drinking water for at least 1 y (2013), 10 y (2004–2013), or 29 y (1985–2013).

Category	≥1y in Ronneby (n=1,845)	≥10y in Ronneby (n=1,176)	≥29y in Ronneby (n=506)
Women (%)	54.3	55.7	55.9
Median age (y)	45	51	64.5
Age [y (%)]
3–18	19.2	11.1	0
19–65	62.1	64.6	54.3
66–94	18.7	24.3	45.7

Note: PFAS, per- and polyfluoroalkyl substances.

Table 4

Median (5th and 95th percentiles) serum PFOS and PFHxS concentrations of the Ronneby population. Duration exposure regarding those who lived continuously in the areas provided with PFAS-contaminated drinking water for atleast 1 y (2013), 10 y (2004–2013), or 29 y (1985–2013).

PFAS serum concentrations (ng/mL)	≥1y in Ronneby (n=1,845)		≥10y in Ronneby (n=1,176)		≥29y in Ronneby (n=506)
	Females	Males	Females	Males	Females	Males
PFOS	270 (58, 834.5)	290.2 (69.3, 830.3)	367 (85.8, 909.4)	379 (136.1, 900.2)	473.6 (153.8, 1,007)	498.4 (167.2, 999.4)
PFHxS	229 (41.1, 819.2)	263 (52.3, 763)	346.4 (63.6, 897.4)	377.3 (112.3, 869.2)	483.6 (124.6, 995.7)	477.7 (169.7, 971.3)

Note: PFAS, per- and polyfluoroalkyl substances; PFHxS, perfluorohexane sulfonic acid; PFOS, perfluorooctane sulfonic acid.

In each iteration of the analysis, an IEXP distribution was created by randomly drawing, with replacement, from the group of exposed individuals in Ronneby. This IEXP distribution was 1,000 times larger than the original exposed population; for example, 1,845,000 random draws, with replacement, were performed from the sampled 1,845 individuals that lived at least 1 y in Ronneby and drank PFAS-contaminated drinking water. For each scenario, a total of 10,000 IEXP distributions were created.

Estimation of the IBMD

Animal data.

The benchmark dose modeling approach, introduced by Crump (1984), is a statistical method used to describe dose–response relationships. It estimates a benchmark dose (BMD), which is the dose most likely to give rise to a prespecified effect. A 90% confidence interval (CI) is derived, composed of the BMDL and BMDU, the lower and upper bound of the CI, respectively. This CI is estimated by taking into account the uncertainty in the data. were calculated based on data on total levels reported by Seacat et al. (2002). Cynomolgus monkeys (Macaca fascicularis) were exposed once a day for 6 months to 0, 0.03, 0.15, or body weight (BW) per day by intragastric intubation (Seacat et al. 2002). There were six animals of each sex at each dose level, except in the group, where there were only four animals. For the purpose of the current study, the total levels measured by AniLytics (Gaithersburg, Maryland) at the end of the study (day 184) were used as the key end point (Seacat et al. 2002). The group mean serum levels of PFOS at Day 183 were used as a measure of internal exposure (Seacat et al. 2002).

Dose–effect modeling.

The Individual Critical Effect Doses (IBMDs), also known as individual critical effect doses (van der Voet and Slob 2007), were calculated as follows: where is the serum concentration of PFOS leading to a 10% decrease in total levels in cynomolgus monkeys (Seacat et al. 2002). A 10% reduction was used instead of EFSA’s 5% default critical effect size (EFSA Scientific Committee 2017) since such a small change is likely to be without biological significance and within analytical variation (Andersen et al. 2003). According to the paper by van der Voet and Slob (2007), the BMD (and not the BMDL) should be used, corresponding in this case to the best estimate of the serum concentration that leads to the defined critical effect size. The uncertainty in the data is taken into account by the distribution for the BMD, that is, the BMD uncertainty.

To obtain the distribution, the levels were fitted to the serum PFOS levels data from the study by Seacat et al. (2002) using the R package Proast (version 65.5) and following EFSA’s BMD modeling guidance (EFSA Scientific Committee 2017). In the study by Seacat et al. (2002), the PFOS serum levels and the levels were reported on the group level [arithmetic (SD)]. As individual data were unavailable, it was assumed that all monkeys within a dose group had the same serum levels of PFOS, but different levels of at Day 184. Data from males and females were fitted separately, using a 10% decrease in total levels as the critical effect size. The model’s critical difference value, using the Akaike information criterion (AIC), was set to 2, according to EFSA’s guidance (EFSA Scientific Committee 2017). Model 3 of the exponential family had the best fit. The estimated was (90% CI: 2.4, ) and the was (90% CI: 29.3, ) (Figure 2). Then, the bootstrap method (as described by Moerbeek et al. 2004) was employed to estimate 3,000,000 BMDs for males and females, respectively, constituting the distribution. For each iteration, an IBMD distribution was obtained by randomly drawing from the bootstrapped and the extrapolation factor (EF) distributions , , and (further explained below). The number of draws in each IBMD distribution corresponded to the size of the IEXP distribution, keeping the gender proportions.

Figure 2.

Benchmark dose–response analysis for (A) female and (B) male monkeys, using serum PFOS concentrations at day 183 (x-axis) and total levels at Day 184 (y-axis), as described in the study by Seacat et al. (2002). Doses are based on median serum concentrations measured in the dose groups. Note: a, background response, according to the fitted model; AIC, Akaike information criterion; b, potency parameter, according to the fitted model; CED, critical effect dose (also known as benchmark dose); CEDL, lower bound of the CED 90% confidence interval; CEDU, upper bound of the CED 90% confidence interval; CES, critical effect size; conv, convergence, denoted by 1 if the fit algorithm converged and 0 if not; d, steepness parameter, according to the fitted model; dtype, data type, 10 for continuous summary data (expressed in ); loglik, loglikelihood of the fitted model; PFOS, perfluorooctane sulfonic acid; , triiodothyronine; var, the within-group variance (related to the natural log-responses).

Interspecies extrapolation.

There is potentially a difference in toxicodynamic sensitivity to PFOS between the average cynomolgus monkey and the average human, and the uncertainty about this difference was accounted for by the interspecies toxicodynamics extrapolation factor (

Intraspecies extrapolation.

In order to account for interindividual differences in sensitivity to the exposure, an extrapolation factor for the intraspecies variability in toxicodynamics was introduced (). It has been shown that the difference in sensitivity to a specific chemical is often approximately lognormally distributed (WHO/ILO/UNEP 2014). In accordance with data compiled by WHO’s International Programme on Chemical Safety (IPCS; WHO/ILO/UNEP 2014), we assumed a GM of 1 for . The GSD for interindividual variability () of the lognormal distribution is chemical-dependent and the specific value for PFOS is unknown. We assumed that was lognormally distributed and that the 5th, 50th, and 95th percentile of the were 0.0776, 0.221, and 0.631, respectively (WHO/ILO/UNEP 2014). A unique draw for based on the drawn for was performed for each iteration (Figure 1). The distribution is simultaneously a source of variability and uncertainty given that toxicodynamics differences between humans can be described by a variability distribution that is, in turn, subject to uncertainty.

Extrapolation between durations.

The length of a chronic study should cover a considerable part of the life span of the tested species, and cynomolgus monkeys have expected life spans of 25–30 y (Choi et al. 2016). Therefore, the 6-month study performed by Seacat et al. (2002) was considered subchronic. An extrapolation factor taking into account the uncertainty in extrapolating between different durations was used (). As suggested by WHO/ILO/UNEP (2014), was characterized as a lognormal distribution with a GM of 2 and a GSD of 0.84, resulting in a distribution with a 50th percentile of 2 and 5th and 95th percentiles of 0.5 and 8, respectively.

Individual Margin of Exposure Distributions

For each iteration, an IMoE distribution was obtained by dividing a randomly drawn IBMD with a randomly drawn IEXP from the respective distributions (van der Voet and Slob 2007). This was done in accordance with gender, so that IBMDs derived from female monkeys were divided with an IEXP value referring to women in the originally exposed population, and vice versa. Each IMoE distribution contained as many individual values as the IEXP and IBMD distributions. The iteration process was repeated so that a total of 10,000 IMoE distributions were obtained. Then, the percentage of was estimated, constituting the PoCE [probability (Prob) ()] (van der Voet and Slob 2007).

Evaluating Sources of Uncertainty and Variability

The quantification of the relative contribution of each uncertainty and variability source to the PoCE and estimated IMoE distributions was performed as described by van der Voet and Slob (2007). Each uncertainty analysis was repeated while ignoring one or more of the sources of uncertainty, that is, exposure, BMD, , , and . The complete set of possible combinations is represented by a full factorial design. In short, the 10,000-iteration analysis were performed times. When , , or were not included, they were replaced with the median of their distributions. When the BMD uncertainty was excluded, the maximum likelihood estimates of the BMDs for males and females in the original data set were used, that is, for the and for the (Figure 2). When the exposure uncertainty was excluded, the original measures from each of the three exposure scenarios were replicated 1,000 times and used as the IEXP distribution. For each of the 32 combinations, the variances of the PoCE for PFOS as well as for a combined exposure of PFOS and PFHxS were calculated, and the contribution of each of the five uncertainty sources was evaluated with multiple linear regression using the nonnegative linear models (NNLM) package in R (version 3.4.2; R Development Core Team) (Lin and Boutros 2016).

Margin of Exposure Approach

PFOS and PFHxS were the two most abundant PFAS chemicals in the outgoing water from the Brantafors waterworks in December 2013, with concentrations of and (Table 2) (Livsmedelsverket 2013a). It was assumed that these levels reflected the average historic concentrations of PFOS and PFHxS. For deterministic risk assessment purposes, we applied the EFSA-recommended values when using the MoE approach, which suggest the use of for human BW and a standard daily liquid intake of per day (EFSA Scientific Committee 2012). These values were used to estimate the average exposure in the adult population and may, therefore, be underprotective of more highly exposed populations. The MoE was calculated as the ratio of the derived tolerable daily intake (TDI) and an exposure metric for PFOS or PFHxS. A MoE where the numerator is the TDI is also known as margin of safety (MoS) or hazard quotient (HQ). Generally, an exposure below the TDI, that is, an , will be interpreted as being without appreciable risk.

Several TDI values have been proposed for PFOS and PFHxS (Table 1). In 2008, EFSA published a scientific opinion on PFOS (EFSA 2008). A TDI of was established, based on the lowest NOAEL of BW per day derived from the same study (Seacat et al. 2002) used in the present study. A total uncertainty factor (UF) of 200 was applied to the NOAEL: a factor of 100 for inter‐ and intraspecies differences and an additional factor of 2 to compensate for uncertainties related to the extrapolation from subchronic to chronic exposure. EFSA’s CONTAM Panel established a TWI, in 2018, of (), departing from epidemiological data indicating increased serum cholesterol levels, decreased birth weight, and response after vaccination as critical end points (EFSA CONTAM Panel 2018). A new scientific opinion by the EFSA CONTAM Panel has been prepared, estimating a TWI of [week ()] for the sum of four PFAS [PFOS, PFHxS, perfluorononanoic acid (PFNA), and PFOA]. However, this TWI has not been adopted at the time of this publication. In addition to the EFSA, the U.S. Agency for Toxic Substances and Disease Registry (ATSDR 2018) has published a draft for public comment where they derive a minimal risk level (MRL) for PFOS of based on laboratory animal data (delayed eye opening and decreased pup weight in rats, and applying a total uncertainty factor of 300). For PFHxS, there is no EFSA-established TDI, but the EFSA CONTAM Panel was asked by the European Commission to write a scientific opinion on the risk to human health related to the presence of other PFAS in food, including PFHxS (EFSA-Q-2017-00549). The ATSDR (2018) and the Swedish National Food Agency (Livsmedelsverket 2013b) have proposed reference intake values of and , respectively.

Results

Distribution of the Individual Margin of Exposure

The results obtained from an application of the IPRA method applied in the present study to the risk assessment of PFAS are illustrated with cumulative distribution functions (cdf). In Figure 3, the individual exposure distribution (IEXP) is plotted as an inverse cdf, indicating the percentage that exceeds the concentration on the x-axis. Figure 3 left inverse cdf describes PFOS IEXP and Figure 3 right inverse cdf presents the combined exposure of PFOS and PFHxS, with respective 90% CIs. Also in Figure 3 right inverse cdf, the dashed cdf describes the distribution of the human IBMDs, with 90% CIs, indicating the percentage that reaches the critical effect size (10% decrease, on the y-axis) in relation to the necessary serum concentrations (on the x-axis). The IMoE distributions are plotted in Figure 4 (PFOS and PFHxS combined exposure scenario) and Figure 5 (PFOS only). The calculated median PoCE for combined exposure was 2.1% (90% CI: 0.4%–13.1%), meaning that for every 100,000 individuals, between 400 and 13,100 were likely to have a decrease in their total serum levels. For the exposure scenario that considered PFOS exposure only, the calculated median PoCE for the entire population was 1.0% (90% CI: 0.2%–6.2%), meaning that for every 100,000 individuals, between 200 and 6,200 were likely to be affected (Table 5).

Figure 3.

The left inverse cdf describes PFOS IEXP and the right cdf represents the IEXP for the combined exposure of PFOS and PFHxS, with respective 90% confidence intervals. In the same figure, the right dashed cdf describes the distribution of the human individual benchmark doses (IBMD, with 90% confidence intervals), indicating the percentage that reaches the critical effect size (10% , decrease, on the y-axis) in relation to the necessary serum concentrations (on the x-axis). The IMoE distributions are plotted in Figure 4 (PFOS and PFHxS combined exposure scenario).

Figure 4.

Probability of critical exposure [] in the combined PFOS and PFHxS exposure scenario, for the Ronneby, Sweden, residents living in the areas provided with PFAS-contaminated drinking water for at least 1 y (2013). Note: IMoE, Individual Margins of Exposure; PFHxS, perfluorohexane sulfonic acid; PFOS, perfluorooctane sulfonic acid; Prob, probability.

Figure 5.

Probability of critical exposure [] due to PFOS exposure, for the Ronneby, Sweden, residents living in the areas provided with PFAS-contaminated drinking water for at least 1 y (2013). Note: IMoE, Individual Margins of Exposure; PFOS, perfluorooctane sulfonic acid; Prob, probability.

Table 5

Median (5th and 95th percentile) probabilities of critical exposure ( ) for individuals living continuously in Ronneby, Sweden, for at least 1, 10, or 29 y (before 2013).

Category (%)	≥1y in Ronneby (n=1,845)		≥10y in Ronneby (n=1,176)		≥29y in Ronneby (n=506)
	PFOS and PFHxS	PFOS	PFOS and PFHxS	PFOS	PFOS and PFHxS	PFOS
All	2.1 (0.4, 13.1)	1.0 (0.2, 6.2)	2.7 (0.6, 16.5)	1.2 (0.3, 7.8)	3.5 (0.7, 21.8)	1.5 (0.3, 10.4)
Women only	3.9 (0.8, 21.6)	1.7 (0.4, 10.6)	4.8 (1.0, 27.2)	2.1 (0.5, 13.1)	6.2 (1.2, 34.7)	2.6 (0.6, 17.0)
Men only	0.08 (0.02, 2.9)	0.04 (0.01, 0.7)	0.1 (0.03, 3.6)	0.05 (0.02, 1.0)	0.13 (0.03, 5.8)	0.06 (0.02, 1.6)

Note: IMoE, Individual Margins of Exposure; PFHxS, perfluorohexane sulfonic acid; PFOS, perfluorooctane sulfonic acid; PoCE, probability of critical exposure; Prob, probability.

A separate gender analysis revealed that the estimated PoCE was mainly distributed among women (Table 5). For the combined exposure, the PoCE for women were 3.9% (90% CI: 0.8%– 21.6%) and for men 0.08% (90% CI: 0.02%–2.9%) (Table 5). The estimated PoCE for PFOS only was about half of that for the combined exposure, that is, 1.7% for women (90% CI: 0.4%–10.6%) and 0.04% for men (90% CI: 0.01%–0.7%) (Table 5). When the result was separated in relation to duration of residence, the highest PoCE [6.2% (90% CI: 1.2%–34.7%)] was observed among women with nearly three decades of residency, for the scenario combining co-exposure to PFOS and PFHxS (Table 5).

Sources of Uncertainty

The contribution of the different sources of uncertainty for the PoCE in the scenario combining PFOS and PFHxS exposure is described in Figure 6. The multiple linear regression model explained 86.2% of the variance. The most important contributors to uncertainty were extrapolation between subchronic and chronic exposure duration (60.8%), intraspecies toxicodynamic extrapolation (17.8%), interspecies toxicodynamic extrapolation (11.4%), and BMD estimates (10.0%). Monte Carlo and exposure uncertainty sources were estimated to contribute 0.0%.

Figure 6.

The contribution of the different sources of uncertainty to the estimation of the PoCE [ ()] for co-exposure to PFOS and PFHxS, for the residents of Ronneby, Sweden, living in the areas provided with PFAS-contaminated drinking water for at least 1 y (2013). Note: BMD, benchmark dose; PFAS, per- and polyfluoroalkyl substances; PFHxS, perfluorohexane sulfonic acid; PFOS, perfluorooctane sulfonic acid; PoCE, probability of critical exposure; Prob, probability; TD, toxicodynamics.

The contribution of the different uncertainty sources for the overall distribution of IMoEs in the scenario combining PFOS and PFHxS exposure showed that the duration (78.7%) and the interspecies toxicodynamics extrapolations (21.3%) were the two main contributors to uncertainty, whereas the remaining uncertainty sources did not contribute to the uncertainty (see Figure S2).

Results according to the Deterministic Approach

Assuming a daily total liquid intake of 2 L and a bodyweight of for adults [i.e., the default values recommended by the EFSA Scientific Committee (2012)], an average daily intake of and was estimated for the Ronneby population. In relation to the established TDI values, the PFOS estimated MoEs ranged between 0.018 and 1.3 (Table 1).

The exposure to PFHxS was compared with the reference values given by the ATSDR (2018) and the Swedish National Food Agency (Livsmedelsverket 2013b). The estimated MoE ratios were 0.59 and 147, respectively. The estimated exposure would be interpreted as safe if the MoE values are or unsafe if values are .

Discussion

Tiered risk assessments usually begin by using a deterministic approach, such as the MoE, or by comparing the estimated exposure with a reference dose, such as the TDI. Deterministic assessments are simple to carry out, often use readily available data, and produce results that are straightforward to interpret. This low tier is often established as a conservative approach for risk assessment, providing a quick hazard characterization (Meek et al. 2011; WHO/ILO/UNEP 2014). The MoE approach has, however, several shortcomings, such as the unknown degree of uncertainty associated with the estimate and limited resourcefulness in the evaluation of mixture exposures (Meek et al. 2011; WHO/ILO/UNEP 2014). As shown in the present study, the standard interpretation of a situation where the exposure is lower than the reference dose as safe may be an underestimation of the actual risk. Departing from the same data, the MoE indicated no concern, whereas the IPRA approach showed a more appreciable risk. Thus, in situations when the estimated exposure indicates that health risks exceed the acceptable range, high-tier approaches such as the probabilistic analysis might be required. These higher-tier approaches will take the variability in the population more precisely into consideration. In this case, the employed IPRA method starts by describing a dose–response relationship as a mathematical function (i.e., by the use of BMD instead of NOAEL), addressing differences in the exposure (by individual data instead of median or maximum values) and uncertainty (by the use of distributions instead of default uncertainty factors). Finally, it results in a CI that describes the population percentage at risk in relation to a predefined level of effect, given all uncertainties considered.

Group-Specific Risk Estimates with the IPRA

A great feature of the IPRA approach is that it allows estimating the PoCE for subgroups of the study population. In the present study, we estimated the PoCE in relation to gender and exposure duration (i.e., years of continuous residency in the contaminated area) (Table 5). These groups differ both in susceptibility (women are considered to be more sensitive in relation to thyroid effects) and exposure levels (individuals living in the area for many years are more likely to have higher serum levels). Other groups that might be of interest in relation to thyroid effects are different age groups (elderly women have a higher background incidence of hypothyroidism), pregnant women, neonates (thyroid hormones are critical during early development) or groups with specific dietary lifestyles (vegetarians often have lower iodine intake, which may influence the susceptibility in relation to thyroid hormones). Groupwise results may be important from a risk communication and risk mitigation perspective and constitute examples of how the IPRA can be useful for describing the risk for subsets of a population. In the present study, higher PFOS and PFHxS serum levels were confirmed in residents living for longer periods in the areas with contaminated drinking water (Table 5). Given the results, in combination with the higher rate of thyroid disease described among women (and increasing with age) (Andersson et al. 2019; Boelaert and Franklyn 2005) and the fact that the BMD analysis of the study by Seacat et al. (2002) showed a higher sensitivity among female monkeys, we consider women at older age to be at higher risk for chemically induced lowering of total levels.

IPRA Uncertainty Analysis

In the IPRA model, we included as sources of uncertainty the duration extrapolation, inter- and intraspecies toxicodynamic variabilities, exposure, BMD, and Monte Carlo uncertainties. The use of distributions for extrapolations purposes when estimating the risk of adverse effects can reduce the uncertainty of the result, avoiding the pitfall of the uncertainty factors’ arbitrariness.

The biggest contributor to uncertainty in our study was found to be the extrapolation from subchronic-to-chronic exposure ( at the PoCE and 78.7% for the overall IMoE distributions). The distribution for this factor was based on the WHO’s subchronic-to-chronic extrapolation factor that in turn was based on the paper by Bokkers and Slob (2005), which researched mice and rat studies (WHO/ILO/UNEP 2014). This duration extrapolation is likely to be suboptimal for a scenario with monkeys. Moreover, the extrapolation from a 6-month study in primates to a long, yet unknown, exposure in humans carries a high degree of uncertainty.

The toxicodynamic variability () between individuals was the second most important source of uncertainty at the PoCE (17.8%). This variability is both an individual-specific factor, also perceived as sensitivity and is also dependent on the toxicodynamic properties that vary among chemicals. Given that the specific interindividual toxicodynamic values for PFOS and PFHxS are unknown, the surrogate distribution recommended by the IPCS was used (WHO/ILO/UNEP 2014). The interindividual toxicokinetic variability source is absent in the present study, because the authors used internal concentrations, that is, the individually measured serum PFAS concentrations.

The extrapolation factor for interspecies differences in toxicodynamics () relates to how chemicals interact differently with the biological receptors and targets in different species (WHO/ILO/UNEP 2014). The contributed to a lesser extent to the uncertainty of the PoCE estimation (11.4%) and to the overall IMoE distributions (21.3%). The distribution is based on a mouse-to-rat interspecies extrapolation because data are lacking regarding monkey-to-human extrapolation factors (WHO/ILO/UNEP 2014).

BMD uncertainty contributed 10% to the PoCE estimation and 0% to the overall IMoE distribution uncertainty. This uncertainty source is related to the quality of the dose–response data. The actual BMD contribution to the overall IMoE distribution uncertainty is likely to be higher than 0%. The underestimation is partly related to the fact that the multiple linear regression method does not explain 100% of the variance.

The Monte Carlo uncertainty source contributed 0% to the uncertainty for PoCE estimation and to the overall IMoE distribution. This was mainly explained by the very high number of Monte Carlo runs.

Generally, uncertainty contributions can be reduced by employing distributions based on more precise data as, for example, duration extrapolation and inter- and intraspecies toxicodynamic distributions (van der Voet and Slob 2007). The WHO’s IPCS ranked the general contributions to uncertainty as “1. Use of a NOAEL as the POD; 2. Intraspecies variability; 3. Duration extrapolation; 4. Interspecies TK/TD differences; 5. Interspecies body size scaling” (WHO/ILO/UNEP 2014). The present IPRA analysis resulted in a slightly different ranking (). It is clear that the use of BMD modeling, considering all data in the dose–response curve, decreases the uncertainty. Even more precise results would have been achieved if a chronic toxicity study had been available and if we had a better understanding of the interindividual toxicodynamic variability. These areas are therefore in need of further research.

Comparisons between the Deterministic and Probabilistic Approaches

For PFOS, the deterministic ratio between estimated daily intake and the TDI ranged from to 60, depending on the selected guideline value. Similarly, the risk described for PFHxS is highly dependent on the selected reference dose, but generally showing values of less concern. The results of deterministic approaches are often used to guide prioritization and to perform screening assessment. The results are easy to communicate as they characterize the average exposure as safe or unsafe. Even though the most recent TDI values indicate that the exposure in Ronneby constitute an appreciable risk for human health, the vagueness of the estimate fails to describe the fraction of the population that might be at risk of having adverse effects (i.e., a 10% decrease in levels) due to the exposure.

As shown in Table 1, the results obtained by the deterministic approach are highly dependent on the selected POD and arbitrariness of the uncertainty factors applied to it. In the present study, we used the same key data and critical effect as EFSA (2008). Our distributions for extrapolation factors are also similar to the default uncertainty factors used by EFSA. For example, for duration extrapolation we used a distribution with the GM of 2 (with 5th–95th percentiles of 0.5 to 8, respectively), whereas the EFSA CONTAM panel used a fixed uncertainty factor of 2. Following the argumentation for a tiered risk assessment approach (Meek et al. 2011), the deterministic MoE should be considered as a more conservative approach than the higher-tier probabilistic measures, that is, only if the exposure exceeds the TDI is there a need for a higher-tier probabilistic approach. To compare the measured serum levels with an intake-based reference value such as TDI, we used the measured levels of in drinking water and estimated an intake of , which is below the TDI derived for by EFSA for PFOS in 2008 (EFSA 2008). This contrasts with the IPRA-based result indicating a median PoCE for women (living in Ronneby for at least 1 y) of 1.7% (90% CI: 0.4%–10.6%), that is, up to 1 of 10 are likely to have an effect. The limited ability of deterministic approaches to characterize uncertainty and variability may in fact lead to an underestimation of the risk given that “without risk” (MoE approach) and “up to 1 of 10” (IPRA approach) represent two very different narratives.

The contrasting results obtained by the lower-tier deterministic MoE approach and the higher-tier IPRA indicate a more pronounced risk by the latter. This might be surprising as a probabilistic analysis is considered to be less conservative (WHO/ILO/UNEP 2014). The main reason for the underestimated risk by the deterministic approach seems to be related to the misconception of the NOAEL as a zero effect level. The NOAEL is highly dependent on study design and statistical power. In the study by Seacat et al. (2002), the number of animals per group was rather low () and thus is the power of the study. The pairwise comparison (t-test) reported showed significance at the high- and middle-dose groups, leading the EFSA CONTAM Panel to interpret this as a likely absence of effect at lower doses, establishing a TDI of based on the study’s low dose of (EFSA 2008). Because the BMD analysis takes all data and variation into account and expresses the uncertainty through a CI, we suggest that the BMD should be considered as the preferred approach in risk assessment, as also proposed by an increasing number of expert opinions (EFSA Scientific Committee 2017; NRC 2001).

In summary, the traditional deterministic risk assessment procedures address the problem of uncertainty by making use of conservative or worst-case scenarios. Using conservative scenarios to deal with uncertainty is not desirable because they are purported to be conservative, but are not necessarily so. A conservative scenario is considered to be unrealistic and may result in an overly conservative risk assessment leading to unnecessarily stringent risk mitigation. In addition, it is impossible to explicitly quantify how conservative the risk assessment is. As shown in the present study, there is a substantial possibility that the deterministic approach MoE underestimates risks, specifically when a NOAEL-value is used as a POD.

The IPRA approach is especially useful when there is a concern, such as when the estimated exposure exceeds the TDI. Our results indicate that longer residence times, especially among women, lead to a greater risk of reduced levels of total . This effect might not be observable at an individual level, but a probabilistic analysis certainly results in more precise estimates of the population at risk of adverse effects.

In the present study, we assumed that the internal dose PFOS and PFHxS shared the same potency on a molar basis. Given the structural similarities, it is likely that PFHxS shares the toxicological profile of PFOS, especially when exposure is assessed in serum rather than in intake. The assumptions of similarity are supported by the similar order of magnitude of the estimated relative potency factors (RPF) for PFOS and PFHxS estimated by the Dutch National Institute for Public Health and the Environment (Zeilmaker et al. 2018) and Massachusetts Department of Environmental Protection (MassDEP 2019). By using internal exposure measures (serum levels), toxicokinetics differences have been circumvented, remaining only the toxicodynamic differences that have been addressed by the usage of extrapolation factor distributions and . In addition, other PFAS were also present in Brantafors waters, but at considerably lower concentrations than PFOS and PFHxS (Table 2). We excluded these chemicals from the present analysis as the serum levels were also significantly lower than those of PFOS and PFHxS (Li et al. 2018).

The impact of other potential sources of uncertainty, such as the measurement error in the PFAS quantification and the uncertainty regarding the exposure metric, cannot be disregarded. The human PFAS-exposure was based on a single serum determination, which was performed at the earliest 6 months after exposure termination. Given the long average half-lives for these compounds, 3.4 y for PFOS and 5.3 y for PFHxS (Li et al. 2018), the chronic exposure described and the great number of combinations performed, it was assumed that the variation in the exposed population was reflected in IEXP distributions and on the analysis performed. Moreover, exposure did not contribute significantly to the uncertainty of this analysis (Figure 6; see also Figure S2).

Risk assessment of PFAS mixtures is an emerging issue. As exposure to these chemicals mainly occurs in complex mixtures, a risk assessment approach that can combine different chemicals in an integrated manner should be preferred. Both the MoE and the IPRA have been developed for single substances but can be applied for mixtures by, for example, using relative potency factors. How to approach risk assessment of complex mixtures has been discussed elsewhere (Kienzler et al. 2016; Meek et al. 2011). In this probabilistic analysis, we showed how the PoCE estimate differs when considering exposure to PFOS alone or to PFOS and PFHxS combined. We consider the latter to be more relevant and, thus, combined exposure of chemicals assumed to act on a similar mode of action should be assessed groupwise using dose addition.

Final Remarks

Taken together, our findings suggest that the IPRA can be used to evaluate the risk for decreased thyroid hormone levels related to exposure of PFOS and PFHxS. We find that the highest risk is likely to be found among elderly women with long-term residency in the contaminated area. In addition, we show that deterministic methods, although they are considered to be conservative, are likely to underestimate the risk and our study exemplifies a shift of narrative following the use of a probabilistic approach compared with the traditional deterministic approach. We conclude that probabilistic risk characterization represents an important step forward with the ability to adequately analyze group-specific health risks such as in infants, pregnant women, and individuals with long term residency history in the contaminated area. Moreover, quantifying the uncertainty and its underlying sources is desirable as it will complete the hazard characterization by improving the awareness of uncertainties among stakeholders, leading ultimately to more accurate chemical risk assessment and risk management.

Click here for additional data file.

Click here for additional data file.