Physiologically based pharmacokinetic modeling of perfluoroalkyl substances in the human body

Currently, there are limited data on the levels of perfluoroalkyl substances other than perfluorooctane sulfonic acid and perfluorooctanoic acid in the human body. Most of this information has been extracted from biological monitoring of plasma while the occurrence of perfluoroalkyl substances in other human tissues is rarely studied. The objective of the present study was to develop a physiologically based pharmacokinetic model to assess the concentration of perfluoroalkyl substances in human tissues, based on an existing model previously validated for perfluorooctane sulfonic acid and perfluorooctanoic acid. Experimental data on concentrations of perfluoroalkyl substances in human tissues from individuals in Tarragona County (Catalonia, Spain) were used to estimate the values of some distribution and elimination parameters needed for the simulation. No significant correlations were found between these parameters and the chain lengths. The model was finally validated for five perfluoroalkyl substances.


Introduction
In recent years, perfluoroalkyl substances (PFASs), a group of organic compounds, have attracted social and scientific attention. PFASs are highly persistent in the environment, bioaccumulative in living organisms, potentially toxic, as well as being subject to longrange atmospheric transport and categorized as persistent organic pollutants (POPs) (Wang et al. 2009). Many PFASs have been detected worldwide in a wide range of environmental and biological samples (Butt et al. 2010;Domingo 2012;Fromme et al. 2009;Weiss et al. 2013). Due to the high strength of carbon fluorine bonds, PFASs are not metabolized and poorly eliminated from the human body (Lau 2012;Lau, Butenhoff, and Rogers 2004). PFASs are mainly distributed to well-perfused tissues such as plasma, liver, and kidney (Squadrone et al. 2014). Half-lives in human blood have been estimated to be 5.4 years for perfluorooctane sulfonic acid (PFOS), between 2.3 and 3.8 years for pefluorooctanoic acid (PFOA), and 8.5 years for perfluorohexane sulfonic acid (PFHxS) (Lau 2012). Shorter half-lives have been estimated for perfluorohexanoic acid (PFHxA), perfluorobutane sulfonic acid (PFBS), and perfluorobutanoic acid (PFBA) (32, 30, and 3 days, respectively) (Lau 2012;Olsen et al. 2009;Russell, Nilsson, and Buck 2013). Toxicological studies of PFASs have been conducted in respect to hepatotoxicity, reproductive toxicity, and immunotoxicity, as well as their effects on body weight, development, and on the levels of cholesterol and thyroid hormones (Lau et al. 2007;Stahl, Mattern, and Brunn 2011).
In 2002, the US company 3M, the primary global manufacturer of PFASs, phased out the production of PFOS (Wang et al. 2009). In 2009, this compound was included in the annex B of the Stockholm Convention restricting its production and use, while PFOA became a serious candidate to enter the annex. In 2006, the U.S. EPA created the PFOA 2010/2015 stewardship program aimed at reducing 95% of the PFOA releases and product content by 2010, as well as completely eliminating its production in Western countries by 2015 (U.S. EPA 2014). Due to these regulatory restrictions, the producers of PFASs started to look for alternatives to PFOS and PFOA. The company 3M has been developing a new line of products based on PFBS (Renner 2006). In turn, other companies became partner in the FluoroCouncil (Bowman 2015), aimed at developing a new line of products based on PFHxA (Buck et al. 2011). Due to the higher elimination rates from humans and the lower toxicity of PFBS and PFHxA, PFOA has been progressively replaced in industrial applications (Russell, Nilsson, and Buck 2013;Olsen et al. 2009). PFOS and PFOA have been extensively studied in regard to their pharmacokinetics, but there is a lack of such studies for other PFASs (Borg et al. 2013).
Similarly to other substances, the use of experimental animals has been recurrent to identify the main toxic endpoints of PFASs and their distribution in the body (Inoue et al. 2012;Kowalczyk et al. 2012). Human data are particularly scarce because of the limitation to obtain tissue samples (Olsen et al. 2003;Maestri et al. 2006). Alternatives to animal testing are promoted to avoid or minimize the use of living animals, with in silico tools being promising alternatives (Benfenati et al. 2010). PBPK models may be used to elucidate the pharmacokinetics (the study of the time-course absorption, distribution, metabolism, and elimination (ADME)) and pharmacodynamics (the study of the site of action and the resulting effects) (PK/PD) of chemicals in the human body (Medinsky 1995;Wang et al. 2013). Recently, various PBPK models have been developed to estimate the distribution of PFOS and PFOA within the human body (F abrega et al. 2014;Loccisano et al. 2011;Loccisano et al. 2013;Andersen et al. 2006;Tan, Clewell III, and Andersen 2008). It has been suggested that the most important mechanism of removal of PFOS and PFOA in the human body is urinary elimination, with a renal resorption mechanism in the filtrate compartment (Andersen et al. 2006;Tan, Clewell III, and Andersen 2008). The free fraction (the fraction not bound to plasma albumin) of PFOS and PFOA is also an important parameter to determine the distribution of PFOS and PFOA. However, there still exist some gaps and notable uncertainties around the PK/PD of PFASs different from PFOS and PFOA.
Based on the previously obtained data for PFOS and PFOA (F abrega et al. 2014), a PBPK model to estimate the concentrations of nine additional PFASs in different human tissues was developed. Overall, up to 11 PFASs with different chain lengths were considered for the PBPK modeling: PFBS, PFHxS, PFOS, perfluorodecane sulfonic acid (PFDS), PFHxA, perfluoroheptanoic acid (PFHpA), PFOA, perfluorononanoic acid (PFNA), perfluorodecanoic acid (PFDA), perfluoroundecanoic acid (PFUnDA), and perfluorotetradecanoic acid (PFTeDA). The tissues studied in the present model were selected according to their pharmacokinetic relevance in the distribution and elimination of PFASs. Plasma, gut, liver, fat, kidney, filtrate, bone marrow, brain, and lungs were selected as compartments of interest, while the rest of the body was integrated into an additional single compartment. Experimental data on PFAS concentrations in human tissues from individuals in Tarragona County (northeast of Spain) were used to estimate the values of some distribution and elimination parameters needed for the simulation. The correlations with the chain length of each respective PFAS were also studied. Afterwards, the model was validated by using data of intake and plasma concentration from Andøya Island (Northern Norway) (Rylander et al. 2009;Haug et al. 2010).

PBPK model for PFASs
The conceptual structure of the PBPK model is depicted in Figure 1. The intake of chemicals was considered to occur only by food and water consumption through the gut. Kidney was selected because it is the organ where elimination takes place (Han et al. 2012). Liver, bone marrow, kidney, lungs, and brain were chosen because quantifiable concentrations of PFASs have been found in autopsy tissues (P erez et al. 2013). When the percentage of detection was 0, the value used for validation purposes was half of the detection limit. Finally, plasma was used as a descriptor of the systemic circulation. The PBPK model was based on flow-limited equations (Thompson and Beard 2011).
The concentrations of PFASs in non-eliminating tissues (bone marrow, fat, brain, lungs, and the rest of the body) were estimated as follows: where C i is the concentration in the tissue i (ng L ¡1 ), Q i is the blood flow to the tissue i (L h ¡1 ), free is the unbound fraction in the plasma, C a is the arterial concentration (ng L ¡1 ), K p is the partition coefficient, and V i is the volume of the tissue i (L). The concentrations of PFASs in the gut were calculated by means of the following expression: where C g is the concentration in the gut (ng L ¡1 ), Q g is the blood flow to the gut (L h ¡1 ), free is the unbound fraction in the plasma, C a is the arterial concentration (ng L ¡1 ), K p is the partition coefficient, Intake is the hourly ingestion of PFASs (ng h ¡1 ), and V g is the volume of the gut (L). For the liver compartment, the concentrations of PFASs were estimated as follows: where C l is the concentration in the liver (ng L ¡1 ), free is the unbound fraction in the plasma, Q l is the blood flow to the liver (L h ¡1 ), C a is the arterial concentration (ng L ¡1 ), Q g is the blood flow to the gut (L h ¡1 ), C g is the concentration in the gut (ng L ¡1 ), K p,g is the partition coefficient in the gut, K p,l is the partition coefficient in the liver, and V l is the volume of the liver (L). For the kidney compartment, the equation used was as follows: where C k is the concentration in the kidney (ng L ¡1 ), Q k is the blood flow to the kidney (L h ¡1 ), free is the unbound fraction in the plasma, C a is the arterial concentration (ng L ¡1 ), K p is the partition coefficient in the kidney, T m is the resorption maximum (ng h ¡1 ), C fil is the concentration in the filtrate compartment (ng L ¡1 ), K t is the affinity constant (ng L ¡1 ), and V k is the volume of the kidney (L). Finally, the levels of PFASs in the filtrate compartment were estimated by means of the following mathematical expression: where C fil is the concentration in the filtrate compartment (ng L ¡1 ), Q fil is the blood flow to the filtrate (L h ¡1 ), free is the unbound fraction in the plasma, C a is the arterial concentration (ng L ¡1 ), T m is the resorption maximum (ng h ¡1 ), K t is the affinity constant (ng L ¡1 ), and V fil is the volume of the filtrate compartment. Calculations were performed for the adult population of Tarragona County. Water intake was assumed to be 1.0 L d ¡1 (Nissensohn, Castro-Quezada, and Serra-Majem 2013), the mean body weight was set at 70 kg. The daily intake of PFASs through food and water for the population of Tarragona County (Table 1) was obtained from Ericson et al. (2009) and Domingo et al. (2012). A constant intake was considered for the lifespan of the individuals under study. The mean lifetime of the individuals was set to 80 years.
The distribution of PFASs in the human body was described by the partition coefficient (P k ), which is defined as the concentration of a chemical in a specific tissue in relation to its concentration in blood. To calculate P k , the concentration of each PFAS found in autopsy tissues from subjects of Tarragona County was divided by the mean blood level of the same compound found in samples of donors from the same area of study (P erez et al. 2013;Ericson et al. 2007). P k s were assessed for liver, bone marrow, brain, lung, and kidney. For the remaining tissues (adipose, gut, filtrate, and the rest of the body), a constant value of the P k s, coming from the data on PFOA in rats, was used (Kudo et al. 2007). The concentration in these tissues is not very high and, as a consequence, the added uncertainty is not likely to be incremented.
The elimination parameters of the PBPK model, resorption maximum (T m ) and affinity constant (K t ), and the free fraction of PFASs in plasma were also assessed for the case study of Tarragona County. The data on the body burdens of PFASs in plasma, as well as those in a number of tissues (liver, bone, kidney, lungs, and brain), were available from previous studies (P erez et al. 2013). PFASs were considered to be bound to serum albumin in a fraction of 97%, while the free fraction (3%) was available for accumulation in other tissues and elimination through urine (Chen and Guo 2009;Bischel, Macmanus-Spencer, and Luthy 2010;Han et al. 2003). In order to assess the elimination parameters (T m , K t ) and free fraction, the model was coded and parameterized, being the simulated plasma concentration fitted with experimental values of PFASs (Ericson et al. 2007). Values of T m , K t , and free fraction constants were obtained here. The model was coded by using Berkeley Madonna TM v 8.3.18, using Stiff method to solve the differential equations (Macey, Oster, and Zahnley 2009).

Validation of the PBPK model
The validation step consisted in executing the PBPK model for another population. Adults living in Andøya Island (Norway), for which data on intake and plasma concentrations of PFHxS, PFOS, PFHpA, PFOA, and PFNA were available (Rylander et al. 2009;Haug (Huizer et al. 2012) and confirmed by sensitivity analysis of our PBPK model (unpublished results), T m , K t , and free fraction were the variables with the highest sensitivity in the PBPK model. Minimum and maximum values for T m and K t were set as 0.3 times the coefficient of variation (CV) (Allen, Covington, and Clewell 1996). In turn, 0.1% and 3% were considered as the minimum and maximum percentages of the free fraction, respectively (Huizer et al. 2012). The daily intake of food and water consumption by the adult Norwegian population was set at 1.2, 18, 8.2, 31, and 9.5 ng day ¡1 for PFHxS, PFOS, PFHpA, PFOA, and PFNA, respectively (Haug et al. 2010). Levels of PFASs in plasma were taken from a biological monitoring study of 44 women and 16 men living in the coastal population of Northern Norway, which was performed in September 2005 (Rylander et al. 2009). Mean plasma concentrations of PFHxS, PFOS, PFOA, and PFNA were 1.8, 33, 4.4, and 0.95 mg L ¡1 , respectively. Although PFHpA was not detected, the limit of detection (0.26 mg L ¡1 ) was used to validate the model. Finally, model validation was based on the comparison of experimental data on PFASs in plasma and the range of values predicted by the PBPK model.

Results and discussion PBPK model for PFASs
A number of assumptions were taken into account before the PBPK modeling development and application. All the variables, including volumes, cardiac output, partition coefficient (P k ), intake, and elimination parameters, were considered to be constant along time. Data on time volumes and cardiac output were taken from Brown et al. (1997). Due to their high absorption rates observed in rats, PFASs were assumed to be completely absorbed (Hundley, Sarrif, and Kennedy 2006;EFSA 2014). In turn, since PFASs under study are poorly metabolized in animals, metabolism was considered negligible (Ophaug and Singer 1980;Hundley, Sarrif, and Kennedy 2006). Although other exposure pathways, such as air inhalation and dermal absorption, might have a relevant contribution to the total intake of PFASs (Ericson Jogsten et al. 2012; D'Eon and Mabury 2011), water intake and food consumption were considered as the only routes of PFAS entrance to the human body. For confirmation purposes, values of dust ingestion of PFASs from the population living in the same area of study were compared with the dietary and water intake of PFASs, being the former considerably lower (Ericson Jogsten et al. 2012). In addition, since PFASs were not detected in air samples, air inhalation was considered negligible (Ericson Jogsten et al. 2012). The partition coefficients of 11 PFASs are summarized in Table 2. The P k values ranged from 0.001 for PFDS, PFDA, and PFTeA in the liver, as well as PFTeA in the bone marrow, to 201.6 for PFHxA in the brain. The correlation between the carbon chain length and those physicochemical and biochemical parameters of each respective PFAS was also studied. No significant correlations were noted between the PFAS chain length and P k (Figure 2). However, the available experimental information is relatively scarce, since PFAS concentrations from only 20 individuals had been obtained.
Based on the results of animal studies (Andersen et al. 2006;Tan, Clewell III, and Andersen 2008), the elimination of PFOS and PFOA is mainly urinary, with a resorption mechanism in the filtrate compartment. Once in the urine, the chemicals are resorbed back to the plasma in a saturable mechanism, driven by resorption maximum (T m ) and affinity constant (K t ). T m and K t are analogous constants to the kinetics of the MichaelisÀ Menten reactions for the enzymatic reactions. This mechanism was expected to be responsible of the high half-lives of PFOS and PFOA in human blood. However, it has not been proven yet to be involved in the elimination of shorter chain PFASs. In the present study, a resorption mechanism was assumed for all the compounds, including C4 and C6 shorter chain PFASs (Han et al. 2012). The elimination constants (T m and K t ), and the value of free fraction of PFASs in plasma, were calibrated by running multiple simulations of the PBPK model until fitting the simulation results with experimentally  obtained PFAS concentrations in plasma. The elimination constants resulting from the PBPK modeling are summarized in Table 2. The final results for T m ranged from 6.1 to 2456.3 mg h ¡1 , K t ranged from 0.018 to 30 mg L ¡1 , and free fraction ratios ranged from 0.001 to 0.03. Significant correlations between the carbon chain length and any of these constants were not found (p > 0.05) (Figure 2). The model was able to simulate the pharmacokinetic behavior of PFASs with a shorter chain than PFOS and PFOA, such as PFBA or PFBS, whose half-lives in blood are considerably lower. To the best of our knowledge, this is the first attempt to develop in humans a PBPK model for PFASs other than PFOS and PFOA. According to the present PBPK model, the urinary elimination, the resorption mechanism, and the strong protein binding are the hallmark of the ADME of PFASs in the human body. Although other mechanisms of elimination may play an important role (Harada et al. 2007), data are currently too scarce to accurately evaluate their influence. Moreover, this is a model that should be improved by including more refined individual pharmacokinetics data of each individual chemical compound. Anyhow, it is undoubtedly of great interest as information on the body burdens of PFASs, other than plasma, is not available. The estimated steady-state concentration of PFASs in target tissues is shown in Table 3. The simulations followed a trend with a nearly lineal part at the beginning, reaching a plateau after 20À30 years. The result of the simulations found in tissues depends on the P k , as well as on the elimination constants and the daily intake. In the liver, the highest concentrations corresponded to PFOS, PFHxA, PFHpA, PFOA, and PFNA, with values ranging from 1.32 to 127.6 ng g ¡1 (PFNA and PFOS, respectively). In contrast, the minimum values of other long-chain PFASs (PFDS, PFDA, PFTeDA, and PFUnDA) were also found in the liver (0.0002, 0.0005, 0.0005, and 0.0016 ng g ¡1 , respectively). This fact is surprising, as liver is usually considered the main target organ where PFASs are accumulated (Yeung et al. 2013). Furthermore, there are no evidences of low concentrations of long-chain PFASs in human liver (Karrman et al. 2010). PFOA was the predominant compound in bone marrow and lung (205.30 and 99.70 ng g ¡1 , respectively). Unexpectedly, the highest concentration of PFHxA (39.19 ng g ¡1 ) was estimated in the brain. Since PFHxA is rapidly cleared from biota, with elimination half-lives of 0.5À1.5 months in humans (Russell et al. 2013), high levels of this compound should For model parameterization, data on PFAS concentrations in blood and other human tissues from two independent studies, performed in Tarragona County in 2007 and 2008, were used (P erez et al. 2013Ericson et al. 2007). The number of analyzed autopsy samples was limited to 20, therefore, contributing to a higher uncertainty of the PBPK model results. There is a lack of reliable experimental data for their use in model parameterization. Thus, limited available data may contribute to higher uncertainty in the PBPK model due to temporal variability (different years of sampling) and physiological variability (age group differences) associated with the experimental data. Moreover, sources of variability, like changes of body weight and cardiac output along time, which were not taken into account, should be explored in further investigations.

Validation of the PBPK model
For validation purposes, the current PBPK model was applied in a case study in Norway, where information about dietary intake and plasma concentration of a number of PFASs was available (Haug et al. 2010;Rylander et al. 2009). We assume that the data of ingestion found by Rylander et al. (2009) is representative of the dietary intake of the average population of Norway. The concentrations of PFHxS, PFOS, PFHpA, PFOA, and PFNA were assessed in 21 samples of selected foods and beverages, including meat, fish, vegetables, fruits, eggs, milk, cereals, bread, and drinking water (Haug et al. 2010). The samples were purchased between 2008 and 2009 in Norway, and the daily consumption of foodstuffs for the adult population was assessed by using statistical data. The total intake of PFASs for the adult Norwegian population was estimated in 1. 2, 18, 8.2, 31, and 9.5 ng day ¡1 for PFHxS, PFOS, PFHpA, PFOA, and PFNA, respectively. In turn, plasma concentrations of PFASs for the Norwegian population were also obtained from the scientific literature (Rylander et al. 2009). Samples from 44 men and 16 women (aged 26À60) living in Andøya Island (Northern Norway) were collected, and the PFAS content was determined. Mean concentrations of PFHxS, PFOS, PFOA, and PFNA in plasma in Norway were 1.8, 33, 4.4, and 0.95 pg L ¡1 plasma, respectively. In turn, PFHpA was not detected in any sample. The experimental concentrations were compared with the results of the PBPK simulation. The ranges of experimental concentrations were reported for all the compounds, with the exception of PFHpA, which was not detected (Rylander et al. 2009). Since the limit of detection (LOD) for PFHpA was provided (0.26 ng g ¡1 ), this value was used to validate the model. The minimum and maximum values for the simulation were assessed by running the model with the minimum and maximum values of T m , K t , and free fraction. These have been identified as the most uncertain parameters of the PBPK model (Huizer et al. 2012;Allen, Covington, and Clewell 1996), being also identified as the most sensitive (unpublished results). The comparison of the simulated and experimental concentrations of five different PFASs in plasma from people living in Norway is depicted in Figure 3. Excepting for PFHpA, which was undetected in plasma, the experimental concentrations of PFASs were slightly higher than the modeled results, therefore, underestimating the internal dose of PFASs. This fact could be linked to a potential interaction among PFASs, as well as between these and other chemicals present in blood. In the past, PBPK modeling has been increasingly used to assess the pharmacokinetics of chemical mixtures (El-Masri, Mumtaz, and Yushak 2004;Dennison et al. 2004). These approaches are based on the comparison of the levels of simulated mixture responses with those anticipated from the individual response addition. The main objective is to determine the interaction threshold, which is defined as the combined total dose of chemicals at which interactions become significant in terms of joint toxicity of a mixture (Mumtaz et al. 2012). More specifically, complex interactive effects of PFOS and PFOA in zebra fish embryos have been observed, being additive, synergistic, or antagonistic according to the mixture ratios of individual chemicals (Ding et al. 2013). Anyhow, since the range of experimental data fell within the range of simulation results, the PBPK model was confirmed to be validated and applicable to PFHxS, PFOS, PFHpA, PFOA, and PFNA, although validity of the model should not be judged on the closeness of the data alone. The model robustness, plausibility, and extrapolation uncertainty should be part of the process (Rescigno, Beck, and Thakur 1987). Consequently, further investigations are still needed in order to incorporate more refined PK/PD data in the model, as well as to reduce its uncertainty.

Conclusions
A PBPK model was here developed and validated to simulate the body burdens of nine PFASs, additionally to PFOS and PFOA. The model structure was characterized by owning a mechanism of urinary elimination, a mechanism of resorption of the chemicals from urine to plasma, and a strong protein binding for all the compounds. This PBPK model was able to simulate compounds with important pharmacokinetic and pharmacodynamic differences, including half-lives in plasma (Lau 2012). Our findings support that the regulation of both short-chain and long-chain PFASs is done by the above-mentioned pharmacokinetic mechanism. Although further evidence is still necessary, this is the first and successful attempt to describe and simulate the pharmacokinetics of PFASs by using these mechanisms. Furthermore, the PBPK model was validated for five of the compounds by means of a case study in Norway. The relatively small differences between the experimental and modeled results are a good indicator of the reliability of the model. However, some biochemical and mathematical aspects, such as the involvement of other elimination mechanisms and the uncertainty of parameterization data, deserve further research. Moreover, because of the scarcity of pharmacokinetics data for PFASs, other than PFOS and PFOA, differences among the compounds could not be studied in depth by PBPK modeling. Since the generation of quantitative time-course data-sets is essential for the validation of PBPK models, an increase of well-conducted human biomonitoring investigations should be also enhanced to assure the validity of these models and their suitability for estimating PFAS body burdens.