Exploring Route-Specific Pharmacokinetics of PFAS in Mice by Coupling in Vivo Tests and Physiologically Based Toxicokinetic Models

Background: Oral ingestion, inhalation, and skin contact are important exposure routes for humans to uptake per- and polyfluoroalkyl substances (PFAS). However, nasal and dermal exposure to PFAS remains unclear, and accurately predicting internal body burden of PFAS in humans via multiple exposure pathways is urgently required. Objectives: We aimed to develop multiple physiologically based toxicokinetic (PBTK) models to unveil the route-specific pharmacokinetics and bioavailability of PFAS via respective oral, nasal, and dermal exposure pathways using a mouse model and sought to predict the internal concentrations in various tissues through multiple exposure routes and extrapolate it to humans. Methods: Mice were administered the mixed solution of perfluorohexane sulfonate, perfluorooctane sulfonate, and perfluorooctanoic acid through oral, nasal, and dermal exposure separately or jointly. The time-dependent concentrations of PFAS in plasma and tissues were determined to calibrate and validate the individual and combined PBTK models, which were applied in single- and repeated-dose scenarios. Results: The developed route-specific PBTK models successfully simulated the tissue concentrations of PFAS in mice following single or joint exposure routes as well as long-term repeated dose scenarios. The time to peak concentration of PFAS in plasma via dermal exposure was much longer (34.1–83.0 h) than that via nasal exposure (0.960 h). The bioavailability of PFAS via oral exposure was the highest (73.2%–98.0%), followed by nasal (33.9%–66.8%) and dermal exposure (4.59%–7.80%). This model was extrapolated to predict internal levels in human under real environment. Discussion: Based on these data, we predict the following: PFAS were absorbed quickly via nasal exposure, whereas a distinct hysteresis effect was observed for dermal exposure. Almost all the PFAS to which mice were exposed via gastrointestinal route were absorbed into plasma, which exhibited the highest bioavailability. Exhalation clearance greatly depressed the bioavailability of PFAS via nasal exposure, whereas the lowest bioavailability in dermal exposure was because of the interception of PFAS within the skin layers. https://doi.org/10.1289/EHP11969


Introduction
Since the 1950s, per-and polyfluoroalkyl substances (PFAS) have been widely used in a variety of industrial and household fields owing to their excellent surface activity, thermal and chemical stability, and hydro-and oleophobic nature. 1 These qualities inevitably led to their ubiquitous occurrence and distribution in global ecosystems. 2 In particular, a growing body of evidence of toxic effects of some PFAS, such as perfluorooctanoic acid (PFOA) and perfluorooctane sulfonate (PFOS), prompt global regulations and limitations on them and on related compounds. 3A couple of PFAS with chain length shorter than PFOA and PFOS, such as perfluorohexane sulfonate (PFHxS), have been mass-produced assuming their lower accumulation potential and toxicities. 46][7] Thus, PFHxS has been proposed for restriction under the Stockholm Convention. 8umans are exposed to persistent organic pollutants mainly through dietary consumption of food and water, inhalation of air, ingestion of particulate matter, and skin contact with dust, clothing, or cosmetics. 9Epidemiological studies and exposure investigations reported that oral ingestion was the main exposure route of PFAS. 10 However, an increasing number of studies have shown that nasal and dermal exposure also made a distinct contribution to PFAS load in the human body, particularly in some special populations. 9,11It was reported that exposure via respiration accounted for 82% and 77% of the total intake of PFOS and PFOA in children, respectively, in a dust-rich environment. 11The detection of PFAS in masks gives another route for inhalation exposure. 125][16] For example, frequent seafood consumption was associated with elevated PFAS serum concentrations in Inuit people in Greenland, 14 whalers in the Faroe Islands, 15 and fishery employees in China. 16 significant correlation was also found between the levels of PFOA in the outdoor activity hand-wipes and human serum.17 High-level PFAS concentrations were observed in the serum of workers who often inhaled aerosols of PFAS.18 However, up to now, there is sparse information on the quantitative relationship between route-specific external exposure levels and internal concentrations of PFAS.Compared with traditional compartment models, physiologically based toxicokinetic (PBTK) model can simulate the absorption, distribution, metabolism and elimination (ADME) of pollutants by properly including description of the physiological parameters of absorption sites and kinetic parameters of the absorption processes, providing a powerful tool to understand the pharmacokinetics via different exposure routes.1 To our knowledge, no dermal or nasal PBTK model for PFAS has been developed, let alone been used in the comparison between different exposure routes.In addition, only PFOS or PFOA was selected as the target compound in most of the previous PBTK models without considering other PFAS such as PFHxS.[22][23][24] Only the proportion of a chemical that reaches the systemic circulation is considered to be bioavailable, which is generally <100% for organic pollutants based on in vivo mouse model.25,26 The bioavailability of a compound is dependent on the exposure routes because the compound undergoes different absorption processes such as digestion in gastrointestinal tract, skin interception, and pulmonary gas exchange.][28][29] However, there is sparse information about the dermal or nasal exposure to PFAS.In a study in which rats received 6 h of nose-only exposure to PFOA-spiked aerosols, the time to reach the peak concentration (T max ) of PFOA in male rat plasma was ∼ 3 h after exposure, 30 which was even quicker than with a single oral exposure (10 h).27 To our knowledge, the bioavailability of PFAS following nasal exposure is unknown.For dermal exposure, one of our recent studies using an in vivo rat model observed that dermal penetration of PFAS in rats with C F > 5 [including six perfluoroalkyl carboxylic acids (PFCA), two perfluorosulfonate acids (PFSA), and two polyfluoroalkyl phosphate diesters (diPAP)] lasted a long period with a T max of 72 h. 31 It suggests that 4.1%-18.0% of the applied PFAS were absorbed into the rat blood 144 h after dermal dosing, which was dependent on their physicochemical properties.31 Due to the few data and inconsistent exposure doses and durations applied in these studies, it is difficult to determine the influence of exposure routes on absorption kinetics, and in-depth studies are warranted.Most current risk assessment studies on PFAS via different exposure pathways were based on the external levels of PFAS in the exposure media, which might overestimate the exposure risks because the pharmacokinetics and bioavailability of PFAS were not considered.17,32 It is imperative to develop models to accurately evaluate the internal exposure levels of PFAS by comprehensively considering the multiple exposure routes and the specific bioavailability.
In this study, PFOA, PFOS, and PFHxS were selected as target compounds, considering that they are highly detected in various environmental matrices, with various functional groups or chain lengths, recalcitrance for biotransformation, 27,33 which is conducive to tracking and comparing their in vivo pharmacokinetics.Single-dose exposure experiments were conducted following separate and joint oral, nasal, dermal, and intravenous administration.The uptake and elimination kinetics of PFAS via different routes and the route-specific bioavailability of PFAS were investigated systematically, which allowed us to calibrate and validate the route-specific PBTK models for PFAS.Then the established models were evaluated under repeated exposure scenarios derived from previous studies, and then we further extrapolated to humans.This study established a reliable platform to quantitatively describe the correlation between external exposure via different routes to internal body loads of PFAS.

Mice Treatment
Male CD-1 mice of 33 ± 3 g were purchased from Beijing Vital River Laboratory Animal Technology Co., and the animal experiments were approved by the Animal Ethics Committee of the Institute of Radiation Medicine, Chinese Academy of Medical Sciences (Approval Number: IRM-DWLL-2021189).The mice were housed in polypropylene mouse cages with a 12 h light/dark cycle and were given free access to regular mouse chow and water.
The mice were randomly assigned into nine groups after 1 wk of acclimatization (Table S2).They were dosed through separate exposure routes, including single gavage (SG), nasal drip (ND), and intradermal injection (ID; outer thigh skin) to simulate oral, nasal, and dermal exposures, respectively.ND introduces liquid into the nostrils, allowing the compound to pass through the respiratory tract autonomously into the lung.In brief, the mouse was held up with their nostrils pointing straight up, and then 50 lL of the mixed PFAS solution was dropped into one nostril through a syringe.To control the applied dosage accurately, ID, which bypasses the process of epidermal penetration, was used.The skin on the outer thigh was carefully shaved before injection.The needle was gently inserted into the skin for a few millimeters while being held virtually parallel to the skin's surface.A bleb will develop following an effective injection.To calculate the bioavailability of PFAS via different exposure routes, a tail intravenous injection (IV) group was set up to calibrate related parameters of the PBTK models (discussed later).The mice were fasted overnight before administration.Mixed PFAS solutions of PFHxS, PFOS, and PFOA at equal concentrations were prepared.Low (1 mg=kg body weight for each PFAS; Table S2, Groups 1-4) and high dosage levels (5 mg=kg body weight for each PFAS; Table S2, Groups 5-8) and a blank control consisting of 0.5% of Tween-20 in pure Milli-Q water were set up for each exposure route, and all treatments were conducted with four replicates (n = 4).Following the guidance for administered volume, [34][35][36] the mice in the SG, ND, ID, and IV groups were dosed with 330 lL of 100 g=L (or 500 g=L), 50 lL of 660 g=L (or 1,320 g=L), 50 lL of 660 g=L (or 1,320 g=L), and 100 lL of 330 g=L (or 1,650 g=L) the mixed solutions, respectively, for the low-level (or high-level) dosed groups so that the dosages in the mice were similar among the different route exposure groups.On 4 h, 16 h, and days 1, 2, 4, 8, 16, and 32 post dosage, the mice were euthanized through cervical dislocation to collect plasma (from the orbital vein), liver, kidney, lung, and skin of the outer thighs (only sampled in orally, nasally, and intravenously dosed groups).The collected tissue was immediately stored at −20 C before analysis.
To simulate the realistic scenario where organisms and humans are usually exposed to PFAS via multiple exposure pathways simultaneously, another group of mice was exposed to PFAS through the three routes concomitantly, each at a dosage of 1 mg=kg body weight (Table S2, Group 9) (n = 4).Plasma, liver, and kidney were sampled on days 2, 16, and 32 and processed following the procedure described above.

Extraction and Instrumental Analysis
The mice livers, kidneys, and lungs were homogenized with phosphate-buffered saline (PBS) by a Freezer Mixer (1:3, 1:1, 1:2 wt/vol, respectively).The mouse skin was cut into pieces with anatomical scissors.Before sample extraction, 1 ng of isotopically labeled standards were added to a 15 mL polypropylene tube that contained the homogenized tissue samples (10 lL of plasma, 10 lL of liver homogenate, 20 lL of kidney homogenate, 50 lL of lung homogenate, 100 mg of skin, respectively). 37ext, 1 mL of 0.5 M TBAH (pH = 10), 2 mL of 0:25 M sodium Environmental Health Perspectives 127012-2 131(12) December 2023 carbonate buffer, and 5 mL of MTBE were added and then shaken at 250 rpm for 20 min.The organic phase was collected in a new PP tube after centrifuging at 6,000 rpm for 10 min.The extraction process was repeated, and the organic extracts were combined and then almost dried out under nitrogen.Following that, 5 mL of methanol was added to reconstitute the mixture.A Cleanert PEP cartridge (Bonna-Agela Technologies) was used to purify the concentrated extract after it had been preconditioned with 10 mL of methanol and 10 mL of water.This process was followed by washing with 5 mL of 20% methanol in water.The target chemicals were pumped dry and then eluted with 10 mL of methanol.Under nitrogen flow, the eluate was evaporated to dryness before being reconstituted in 0:5 mL of methanol.After that, samples were filtered for ultra-high performance liquid chromatography-tandem mass spectrometry (UPLC-MS/ MS) analysis using a 0:22-lm filter.The three target PFAS were analyzed on a Vanquish UPLC system combined with a TSQ Altis quadrupole mass spectrometer (ThermoFisher Scientific) with a BEH C18 column (Table S3).

Quality Assurance and Quality Control
PFAS were quantified using isotopically labeled standards (2 ng=mL).The average recovery of PFAS in the spiked samples was 103 ± 6:15% (Table S4).Thus, the concentrations were reported directly without correction by recoveries.Interday and intraday precision and accuracy of the assay were determined by analyzing three quality control samples on the same day (n = 6) and three consecutive days (n = 6).The recovery of intraday and interday injections ranged from 94.8% to 104.6% and 92.6% to 102%, respectively.The coefficients of variation of intraday and interday recovery were 3.60% and 3.78%.
Duplicate analysis was included for every 10 samples, and the average relative standard deviation of the duplicate analysis was 3.70%.A procedure blank and one solvent blank (HPLC grade methanol) were carried out after 10 samples to monitor for any background contamination and instrument carryover.No target compounds were observed in the procedure blanks or solvent blanks.Blind duplicate samples and standards were analyzed for each batch of 10 samples to examine the stability of the instrument.The method detection limit is shown in Table S5.

PBTK Modeling
Description of the route-specific PBTK models.Based on the PBTK models for prediction of PFAS in mice, rats, monkeys, and humans reported in literatures, 24,38 the PBTK models via the individual exposure routes in this study consisted of compartments representing plasma, liver, kidney, lung, dermis, and the rest of the body (Figure 1).Only stomach, intestines, and chamber were included in the oral and dermal exposure groups, respectively.Lung and skin were merely considered as the storage compartments in the nonnasal and nondermal routes, respectively, where no absorption, metabolism or elimination occurred.All the tissues Structures of the PBTK models for PFAS in the mouse following oral (A), nasal (B), dermal (C), and IV exposure.They were modified based on Chou and Lin 38 and Loccisano et al. 24 Note: IV, intravenous injection; K 0C , rate of uptake from stomach into the liver, 1=ðh × BW −0:25 Þ; K absC , rate of absorption of chemical from intestine to liver, 1=ðh × BW −0:25 Þ; K bileC , bile elimination rate constant, 1=ðh × BW −0:25 Þ; K p , the permeability coefficients for dermis, cm/h; K t , affinity constant, ng/mL; K unabsC , rate of unabsorbed dose to appear in feces, 1=ðh × BW −0:25 Þ; K urineC , urine elimination rate constant, 1=ðh × BW −0:25 Þ; K vol , loss rate from chamber, 1/h; P DerCha , dermis-to-chamber partition coefficient, unitless; P LungA , lung-to-end exhaled breath partition coefficient, unitless; Q FilC , the fractional blood flow to filtrate, unitless; Q KC , the fractional blood flow to kidney, unitless; Q LC , the fractional blood flow to liver, unitless; Q LungC , the fractional blood flow to lung, unitless; Q pC , the alveolar ventilation rate, mL/h/g; Q RestC , the fractional blood flow to rest of body, unitless; Q SkinC , the fractional blood flow to skin, unitless; SA, exposure surface area, cm 2 ; T mc , maximal reabsorption rate, ng=h=BW 0:75 .

Environmental Health Perspectives
127012-3 131( 12) December 2023 are perfusion limited and well mixed. 39Biotransformation was not included because of the strong stability of PFHxS, PFOS, and PFOA in organisms. 27Lung and dermis served as a storage compartment for nonnasal and nondermal exposure, respectively.It is acknowledged that PFOA and PFOS are strongly bound to proteins in plasma such as albumin in mammals. 40Han et al. reported that >90% of PFOA was bound to serum albumin in rat and human blood. 40The thermodynamic model assumes that only the free fraction of a chemical is pharmacologically active, 41 which is supported by the fact that the PBTK models based on the free fraction of PFAS usually generated the best fit of the plasma and tissue concentrations. 42,43In this study, the initial value (0.09) of the free fraction was taken from Chou et al. 38 The vascular and extravascular areas are thought to be well stirred with quick equilibrium in the flowlimited model.The venous concentration of a chemical leaving a tissue is equal to its free concentration in the tissue, which is calculated as the total concentration in the tissue multiplied by the quotient of the free fraction and tissue-specific partition coefficient.
The oral absorption of PFAS was described in a twocompartment gastrointestinal (GI) model. 38After absorption, PFAS are transported directly from GI tract to the liver via the portal vein.A first-order constant K 0 (1/h) describes the uptake rate from the stomach, whereas the first-order rate constant K abs (1/h) refers to the uptake in the intestine.The mass balance differential equations (MBDE) describing oral uptake are provided as follows: where R ST and R SI are the change rates of the amount of PFAS (ng/h) in the stomach and small intestine, respectively; A ST and A SI are the amount of PFAS in the stomach and small intestine (ng), respectively; GE is the gastric emptying rate (mL/h); K unabs is the rate constant of unabsorbed PFAS that appear in the feces (/h).Nasal exposure is considered as a single dose to the lung directly. 44Output via exhalation is described by the pulmonary ventilation rate Q p (mL/h) and lung to end-exhaled breath (alveolar air) partition coefficient P LungA (unitless).The MBDE describing nasal uptake is provided and explained below 45 : where R Lung is the change rate of the amount of PFAS (ng/h) in lung; Q Lung is the blood flow rate to lung (mL/h); C A is the PFAS concentration in plasma (ng/mL); C Lung is the concentration of PFAS in lung (ng/mL); P Lung is lung to plasma partition coefficient (unitless); and Free is the free fraction of PFAS in plasma (unitless).Continuous inhalation from air was neglected because the concentration of PFAS in the mice of the control group was extremely low.ID injection creates a chamber on the dermis surface in which the dose is administered, and then PFAS diffuse across the dermis.A two-compartment model separating chamber from the dermis was used to predict the dermal absorption of PFAS. 46Elimination of PFAS via secreted sebum was considered.By incorporating Fick's Law for diffusion, the MBDE for chamber is: where R Cha is the change rate of the amount of PFAS in the chamber (ng/h); K p is the permeability coefficient for the dermis (cm/h); SA is the exposed surface area (cm 2 ); C Der is the concentration in the dermis (ng/mL); P DerCha is the dermis-tochamber partition coefficient (unitless); C Cha is the concentration in the chamber (ng/mL); K vol is the loss rate from the chamber (1/h).
A Cha is the amount of PFAS in the chamber (ng).The MBDE for the dermis is: where R Der is the change rate of the amount of PFAS in the dermis (ng/h); Q Skin is the blood flow rate to the dermis (mL/h); P Der is the dermis-to-plasma partition coefficient (unitless).In one of our recent studies, PFAS were applied on the skin surface of the rats, and the absorption percentages of PFAS were determined (7.74%, 7.21%, and 8.30% for PFHxS, PFOS, and PFOA, respectively). 31When the dosage was multiplied by the absorption percentage obtained in that study, the model was successful to simulate the concentration of the target PFAS in the rat blood, dermis, and liver (Figure S1).To avoid overestimating the internal level due to intradermal injection, the applied dosage was corrected by the absorption percentage when simulating the pharmacokinetics of dermal exposure.The first-order rate constant (K urine , 1/h) was used to describe the urinary excretion rate of PFAS.It has been speculated that renal resorption through transporters in the proximal tubule of the kidney is responsible for the prolonged plasma halflives of PFAS. 47Thus, the filtrate compartment described the saturable renal resorption with the maximal reabsorption rate (T m , ng/h) and transporter affinity constant (K t , ng/mL). 22The MBDEs for the liver, kidney, and filtrate are as follows: where R K and R Fil are the change rates of the amount of PFAS in the kidney and filtrate, respectively (ng/h); Q K and Q Fil are the blood flow rates to the kidney and filtrate (mL/h), respectively; C K and C Fil are the concentration in the kidney and filtrate, respectively (ng/g); A Fil is the amount of PFAS in the filtrate (ng); and P K is the kidney-to-plasma partition coefficients.K abs and K 0 are only incorporated in the oral model.The MBDE for the liver is as follows: where R L is the change rate of the amount of PFAS in the liver (ng/h); Q L is the blood flow rate to the liver (mL/h); C L is the concentration in the liver (ng/g); A L is the amount of PFAS in the liver (ng); and P L is the liver-to-plasma partition coefficient (unitless).It was estimated that ∼ 97% of the secreted PFOS was reabsorbed from bile, 48 but this reabsorption was not included in the current model due to a lack of information regarding the relevant parameters.To decrease uncertainty of the model, the "storage" compartment was set to represent the simplified biliary pathway with the first-order biliary elimination rate (K bile , 1/h).
The MBDE for noneliminating tissues, such as lung, skin, and the rest of the body in the oral model, is as following: where R Nt is the change rate in the amount of PFAS in these tissues (ng/h); Q Nt is the blood flow rate to them (mL/h); C Nt is the Environmental Health Perspectives 127012-4 131( 12) December 2023 concentration in them (ng/g); and P Nt is the no-eliminating tissue-to-plasma partition coefficient.Parameterization.][51][52] The fitted chemical parameters are summarized in Table S7.The parameters of blood flows and maximal reabsorption rate were scaled as body weight ðBWÞ 0:75 ; the parameters of the first-order rate constants were scaled as BW −0:25 ; and the partition coefficients, affinity constant, and tissue volumes were scaled linearly to body weight. 53Low-and high-level dosed groups were used to calibrate and validate the PBTK models, respectively.The models were calibrated by fitting the measured concentrations of the target compounds in the plasma, liver, kidney, lung (except nasal route), and skin (except dermal route) of the dosed mice.It was hypothesized that exposure routes did not affect the rates of metabolism and elimination of PFAS. 54,55Therefore, parameters related to renal resorption, the body weight normalized maximal reabsorption rate (T mc and K t ) and bile elimination rate constant (K bileC ) were calibrated by the results obtained from IV (1 mg=kg).Additional calibration of the urine elimination rate constant (K urineC ) was conducted to fit the amount of PFAS eliminated through urine as shown in previous studies. 28,29,56Then they were applied in the calibration of the PBTK models via other routes. 57 preliminary sensitivity analysis was conducted to identify sensitive parameters followed by calibrations using the Levenberg-Marquardt algorithm. 42The purpose was to reduce the computing load and improve the quality of the performance.The routespecific bioavailability of PFAS was calculated by dividing the predicted area under the concentration curve (AUC) via oral, nasal, or dermal route by the AUC via IV. 58Based on the predicted concentration-time data, the elimination half-life (T 1=2 , d) was calculated as follows 59 : where k e is the first-order elimination rate (/d), C end (ng/mL), and C max (ng/mL) are the predicted concentration of PFAS at the end (day 40) and the maximum predicted concentration of PFAS in plasma, respectively.T end (d) and T max (d) is the time at the end point of prediction (40 d in this study) and time to reach the peak concentration.
Model evaluation.By contrasting the predicted concentrationtime kinetic profiles with the experimental results, the performance of the PBPK model was assessed. 60If the predicted values were within a factor of two of the observed values, it was regarded as satisfactory. 60The degree of fit between the log 10 -transformed observed and predicted values was assessed using linear regression.
A local sensitivity analysis was carried out by adjusting each parameter by 1% stepwise from the original value.The normalized sensitivity coefficient (NSC) was determined to examine the impact on the area AUC of plasma (AUC plasma ) because it reflects the relative amount of PFAS absorbed into the systemic circulation. 38The equation to calculate NSC is 38 : where r is the response variable; Δr is the change in the response variable caused by 1% increase in the parameter value; p is the parameter's original value, and Δp is 1% of the parameter's original value.
Extrapolating from Single-to Repeated-Dose Exposure Scenario Because organisms in the environment are usually exposed to pollutants continuously, the PBTK model established from the single-dose scenario was extrapolated to a repeated-dose exposure scenario.We collected the data sets obtained from rat or mice models reported in literature (Table S8) and compared them with the results predicted by the model.The exposure scenarios, body weight, and alveolar ventilation (if necessary) of rats or mice were consistent with those in the literature.2][63] For nasal and dermal exposure routes, only two studies were available for model evaluation of repeated nasal (1 mg=m 3 PFOA, 6 h/d, 5 d/wk for 3 wk) and dermal exposure (20 mg=kg=d PFOA for 10 d). 30,64The original exposure scenario of nasal exposure was adjusted to 953:7 lg=kg=wk for 30 wk (when steady state was reached) derived from the body weight (216 g) and alveolar ventilation rate (31:74 mL=h=g) of the rat to fit our model. 30,49trapolating from Mouse to Human Because oral intake is the main exposure route for general people and there are abundant data on human exposure to PFOS, the established oral PBTK model of PFOS was extrapolated to human assuming that all the estimated daily intake (EDI) enters the human body orally.The physiological parameters of mice were replaced with those of human (Table S9).S10). 650][71] After that, the nasal and dermal models were extrapolated to human by using the Free, K bileC , and rest of body-to-plasma partition coefficient (P Rest ) value calibrated in the oral human model.Poothong et al. determined the EDI values of PFOS through diet and dust digestion, air inhalation, and skin absorption for general people, which were used to estimate the steady-state concentration of PFOS in human blood based on the established human PBTK models. 32The sum of the simulated steady-state concentration of PFOS via each route was considered the total concentration of PFOS in human blood.

Data Analyses
All model simulations were conducted using R (version 4.1.2;R Development Core Team).The model codes were adapted from the code in a previous study. 42The PBTK model was coded in R package "mrgsolve."Linear regression analysis and plotting were performed with GraphPad Prism 8 (GraphPad Software, Inc.).The model code is provided in the supplemental file, "Model Codes.txt."

Physiological Status of Mice
No mortality was observed during the experiment.There was no significant difference between the exposed and control groups in the body weight and liver:body weight ratio, which suggests that no distinct significant toxic effect was observed throughout the experiment (Tables S11 and S12).The concentrations of PFHxS, PFOS, and PFOA in the plasma of the control group were about 2-4 orders of magnitude lower than the exposure groups (Table S13; average 4.63, 2.67, and 6:28 ng=mL, respectively), and therefore contamination from background was considered to be negligible.They were all detected in all the target tissues of the test groups, suggesting that they were absorbed into the tissues through the individual specific exposure routes.

Model Evaluation
Figure 2 and Figure S2 compare the predicted accumulation kinetics of PFHxS, PFOS, and PFOA using the individual PBTK models with the observed results for the purpose of calibration (dose at 1 mg=kg) and validation (dose at 5 mg=kg).The routespecific PBTK models calibrated at a dose of 1 mg=kg well simulated the concentration-time curves at a dose of 5 mg=kg (Figures S3-S14).Overall, 98.8% and 84.2% of the predicted concentrations of PFAS in the tissues dosed at 1 mg=kg and 5 mg=kg, respectively, fell within a predicted-to-observed ratio of 2:1.The predicted concentrations of the three PFAS generated from the individual models correlated with the measured ones very well (R 2 = 0:95-0:98 for the dosage of 1 mg=kg and R 2 = 0:85-0:92 for the dosage of 5 mg=kg).The slopes of the linear fitting equation were very close to 1 (1:02 ± 0:14 for the dosage of 1 mg=kg and 1:11 ± 0:17 for the dosage of 5 mg=kg, p < 0:0001), suggesting the established models are reliable.The predicted amount of these three PFAS in urine also fit the results in previous mice studies very well (Figure S15).At the end of the simulation (40 d), the cumulative excretion of PFHxS in urine was the highest (11.7%-24.0% of the dosage), followed by PFOA and PFOS (2.28%-11.15%,and 0.694%-2.04%,respectively) (Figure S16).The sensitive parameters with NSC >10% for AUC plasma in the low-dosage groups are shown in Figure S17.Free-, body weight-, and liver-related parameters [P L , K bileC , and fractional volume of liver (V LC )] were sensitive parameters common to all exposure routes and all three compounds (Figure S17A and S17B).Specially, the AUC plasma of PFHxS was highly sensitive to the kidney-related parameters [K urineC , T mc , fractional blood flow to filtrate (Q FilC ), and the fractional volume of filtrate (V FilC )].The sensitivity of these parameters and Cardiac output (Q CC ), Hematocrit (Htc), K 0C , and P Rest were highly dose-dependent (Figures S17C and S18C).For nasal exposure, P Lung , P LungA , fractional blood flow to lung (Q LungC ), and the alveolar ventilation rate (Q pC ) were highly sensitive parameters to AUC plasma for all PFAS.For dermal exposure, K vol , K p , and SA were unique sensitive parameters, which only had a high influence on the AUC plasma of PFOS.Gastric and intestinal parameters (K 0C and K absC ) had an obvious effect on AUC plasma of PFHxS for oral exposure.

Pharmacokinetics and Bioavailability of PFAS in Mice
Figure 3A-C illustrates the simulated concentration-time curves of PFHxS, PFOS, and PFOA in the plasma, respectively.The T max of PFAS in plasma (Figure 3D) was the longest in the dermal exposure groups (1.42-3.46d), followed by the oral (0.708-1.71 d) and nasal exposure groups (0.0417 d).For dermal and oral exposure, the order of T max was PFOS (3.46 d)   (2.67 for PFHxS; 2.02 for PFOA; 0.534 for PFOS), respectively.The maximum concentration (C max ) of each PFAS in the plasma via different exposure routes was the highest for oral exposure, followed by nasal, and then dermal exposure (Figure 3E).Among all exposure routes, the concentration of PFHxS in the plasma was the highest, followed by PFOA and PFOS.The concentration of PFOS in the liver was the highest, followed by PFOA and PFHxS.As predicted, PFOS had the highest  S1.Note: IV, intravenous injection.
In the mice receiving joint exposure via oral, nasal, and dermal routes concomitantly, the concentrations of PFAS in the plasma, liver, and kidney on days 2, 16, and 32 were measured (triangles in Figure 4).Figure 4 illustrates that the measured concentrations were very close to the summed concentrations of the three single-exposure routes (squares), suggesting that there was no mutual interference between the exposure routes.Thus, a combination of the three individual models was applied to predict the joint exposure, and the prediction was in agreement with the measured results (triangles dots in Figure 4).

Prediction of the Repeated Exposure in Mouse and Extrapolation to Human
For repeated oral exposure, especially for the 104-wk exposure of PFOS, all the predictions were in good agreement with the observed PFOS concentrations in both the plasma and liver of the mice dosed at 0.024, 0.098, 0.242, and 0:984 mg=kg=d (Figure 5).The model also fit well with the PFHxS, PFOS, and PFOA concentrations in the plasma and liver of the mice dosed at 0.023 and 0:23 mg=kg=d for 7 d orally (Figures S19 and S20) and PFHxS concentration in the rat plasma dosed at 0:3 mg=kg=d for 42 d orally (Figure S21).The modeled concentration of PFOA after repeated nasal and dermal exposure in the rat plasma agreed very well with the measured result (Figures S22 and S23).Overall, 95.0% of the mean ± SD of the measured data under the repeated oral, nasal, and dermal dose scenarios were within the range of 2 times the predicted values.
The results indicate that the human model performed very well in prediction of the PFOS concentration in human blood (Figure 6).The Free, K bileC , and P Rest in the human PBTK model were much lower than that in the mouse model (Table S10).The simulated PFOS concentration in human plasma via separate oral, nasal, and dermal exposure is shown in Figure S24.Oral intake was estimated to make the most contribution to the internal level of  PFOS in the human model.It is worth noting that although the daily intake of PFOS through the nasal route (3:5 pg=kg=d) was 18.4 times that of dermal route (0:19 pg=kg=d), the contribution of the nasal route to the internal level was comparable to that of dermal route (0.0235% and 0.0355%, respectively) (Table S14).

Discussion
In this study, mice were exposed to PFAS via oral, nasal, dermal exposure routes singly or jointly.The route-specific PBTK models were established to investigate the route-specific effects on the pharmacokinetics of PFHxS, PFOS, and PFOA.For the first time to our knowledge, nasal and dermal routes were included in the PBTK models of PFAS.These PBTK models successfully simulated the pharmacokinetics of PFAS in mice after low-and high-level exposures.The results suggest that exposure routes significantly affect the uptake of PFAS in mice, resulting in the highest bioavailability of PFAS after oral exposure, followed by nasal and dermal exposure.The combined model of the three exposure routes performed well to describe the accumulation of PFAS via concurrent multiple-route exposure, which happens in realistic environments.The models were also powerful to predict the concentrations of PFAS in the plasma and various tissues of rodents and humans under repeated exposure in long-term situations.

Route-Specific PBTK Model for PFAS
As shown in Figures 2 and Figure S2, the simulated data obtained from the four individual PBTK models correlated with the observed data very well, which meet the World Health Organization (WHO) model precision criteria, 60 confirming that the established routespecific PBTK models were capable of simulating the singledose exposure of PFHxS, PFOS, and PFOA in mice.Deviations between the predicted and measured values at a few time points might be attributed to the dose-dependent pharmacokinetic process of PFAS in vivo and changes in the physiological parameters of the mice over time. 42The models were used to predict the PFAS concentrations in the serum of mice dosed at 1 mg=kg by Lou et al., 72 Chang et al., 28 and Sundstrom et al. 29 The model  S3.Note: Con, concentration; PFHxS, perfluorohexane sulfonate; PFOA, perfluorooctanoic acid; PFOS, perfluorooctane sulfonate; SD, standard deviation.
Environmental Health Perspectives 127012-8 131( 12) December 2023 well predicted the T max , though it slightly underestimated the observed concentrations in serum (Figure S25).This slight underestimation might be attributed to the mixed PFAS exposure used in this study, leading to competition among the PFAS with different carbon chain lengths and structures when binding to the limited protein cavity sites. 73Similarly, the weaker binding affinity of PFOA to liver fatty acid binding protein (K d-LP = 197 ± 13) than PFOS (K d-LP = 81 ± 7) may account for its lower PL under mixed PFAS exposure in this study than that exposed to PFOA alone in previous studies. 72,74,75UC plasma is a dose-dependent parameter.BW and Free were the most sensitive parameters for AUC plasma of target PFAS and determined the dose and the amount absorbed into the systemic circulation, respectively.The estimated Free (<10%) was comparable to the value obtained in previous studies. 28,38,72In the present study, all the PFAS were sensitive to the distribution and excretion process through the liver.Previous PBTK models were mainly applied to PFOS, which also found that Free, K bileC , and P L were key parameters affecting the PFOS concentrations in plasma and tissues in rodents. 42Renal excretion and reabsorption played a decisive role for PFHxS, which might be because PFHxS prefers to partition in plasma. 29This finding suggests that the models were effective in describing the partition behaviors of PFAS with different physicochemical properties.

Differences in Pharmacokinetics and Bioavailability of PFAS
7][78][79] For example, it was modeled to take 3.79 d for PFOS in the chamber to decrease below 1% of the given dose after dermal exposure, whereas it took only 0.75 and 0.04 d in the stomach and lung after oral and nasal exposure, respectively (Figure S26).This finding suggests that the dermis acts as a barrier to depress the penetration of PFAS; thus, it took a prolonged period for them to be absorbed into the bloodstream.Continuously slow release of contaminants at the sites of dermal absorption could potentially lead to drastically different biological responses than would be expected following other Application of the oral PBTK model to predict PFOS concentrations in plasma and liver in rat that were exposed every day for 104 wk, followed by 1 wk of clearance.Comparisons of the modeled vs. observed PFOS concentrations (mean ± SD) in rat plasma and liver dosed at 0.024 (A1 and A2), 0.098 (B1 and B2), 0.242 (C1 and C2) and 0:984 mg=kg=d (D1 and D2), respectively (n = 5).The experimental data (dots) are from Butenhoff et al. 63 The shaded area represents a predicted-to-observed ratio between 0.5:1 and 2:1.The corresponding numeric data are shown in Excel Table S4.Note: Con, concentration.

Environmental Health Perspectives
127012-9 131(12) December 2023 routes.For example, subcutaneous administration of CPF in the rats greatly prolonged inhibition of cholinesterase activity in the brain relative to oral exposure. 80It was reported that the rate at which pollutants were absorbed through the skin into the blood was limited by the penetration rate. 81After percutaneous exposure, the T max of PFOS in plasma was greater than that of PFOA, which may be attributed to the weaker permeability and lower volatility of PFOS (K p = 1:33 × 10 −5 cm=h, K vol = 3:23 × 10 −7 =h) than PFOA (K p = 1:47 × 10 −3 cm=h, K vol = 6:99 × 10 −7 =h).PFOS is more hydrophobic and has a larger molecular volume; thus, it was more difficult to pass through the water-rich dermis containing living cells. 82he T max of PFAS in blood was the highest after nasal exposure (0.96 h), which was influenced by the physiological structure of the nasal cavity.The olfactory epithelium comprises almost 50% of the total nasal surface area, which would promote the delivery of PFAS to the olfactory bulb, 83 and the abundant capillaries of alveolar epithelial cells expedite the rapid transfer of compounds to blood. 84In addition, anion surfactants were confirmed to be able to promote the nasal absorption of a compound by changing the membrane potential of the mucosa and damaging the close connection between epithelial cells. 85PFAS as anionic surfactants may also be absorbed by the nasal cavity rapidly through these mechanisms.It was modeled that nasal absorption was rapid for many hydrophobic organic pollutants, such as naphthalene, which was independent of storage or metabolism in lung. 86However, this model does not represent absorption of PFAS through the nasal mucosa and upper airways, which may underestimate the absorption rate and proportion of PFAS.
The elimination rate of PFOA was higher than that of PFOS and PFHxS (Figure 3F).The predicted T 1=2 was close to that reported in previous studies using mouse models (30.5-42.8d for PFOS, 24.9-26.8d for PFHxS, and 14.7-24.1 for PFOA) and followed the same order in male rats (7.46-82 d for PFOS, 15.9-28.7 d for PFHxS, 1.83-13.4for PFOA). 27,72It was predicted that renal reabsorption played an important role in their elimination because the T 1=2 was associated with T mc =K t for these three PFAS.It was reported that PFAS were excreted by organic anion transporters (Oat1 and Oat3) and reabsorbed by an organic anion transporting polypeptide (Oatp1a1) in the rats. 87Therefore, the substance-specific reabsorption of PFAS may be related to the different binding abilities with the renal transport proteins.
The modeled results showed that among the three exposure routes, the bioavailability of PFAS decreased following the order of oral >nasal >dermal route.A similar order was observed for 2,2 0 ,4,4 0 -tetrabromodiphenyl ether. 76The route-specific C max is a key factor in determining the AUC and therefore bioavailability (Figure 3H).The highest bioavailability via oral exposure indicated that GI was most efficient of the three routes in absorption of PFAS accompanied by the increased proportions in the blood and tissues (Figures S26A, S27A, and S28A).In comparison with oral exposure, PFAS absorbed via nasal exposure might be subjected to exhalation clearance, thus reducing the final bioavailability of PFAS via this route.Despite the low volatility of PFAS due to their ionic nature, the model predicted that 52.6%, 68.9%, and 38.6% of single-dose PFHxS, PFOS, and PFOA, respectively, were excreted by exhalation after 40 d, suggesting that perfluorosulfonate acids (PFSA) had higher exhalation rate than their carboxylic acids (PFCA) counterparts (Figures S26B,  S27B, and S28B).During exhalation, a large number of bioaerosol particles with large specific surface areas are produced. 88FSA were shown to have strong binding affinities with bioaerosol particles, promoting their exhalation and reduced the overall bioavailability. 89A transport carrier such as this is absent in the dermal exposure 90 ; thus, the predicted loss of PFAS through secreted sebum was relatively low at the end of dermal exposure, which was 12.6% and 2.70% for PFOS and PFHxS, respectively, and PFOA was not eliminated in this way.Combining these results, PFAS due to dermal exposure displayed the lowest bioavailability among the three routes because of their interception by the epidermis.It was reported that hydrophobic compounds are directly delivered into the lipid or protein pool of tissues after oral exposure in a way similar to IV injection, where binding to proteins depresses their removal. 84However, for nasal and dermal exposure, more PFAS exist in free form, and their absorption is mainly via passive diffusion driven by the concentration gradient, favoring interception within the skin layers or exhalation. 84s for the joint exposure via oral, nasal, and dermal pathways simultaneously, for the first time to our knowledge, we report that there was no synergy or antagonism on the accumulation of PFAS in the tissues among different exposure routes.This finding allowed us to evaluate the relative contribution of each route to the internal doses of PFAS.The combined model was capable of predicting the joint exposure results successfully, thus providing a powerful platform to predict the internal body burdens of PFAS due to concomitant multiroute exposure.

Evaluation of the Repeated Exposure of PFAS
In general, the models successfully simulated a majority of the observed data derived from short-, medium-, and long-term repeated dose scenarios via different exposure routes.For the repeated oral and nasal exposures, there was a relatively large error between the predicted and observed values at a few time points, which might result from dynamic changes in physiological parameters and concentration dependence of pharmacokinetics.For the case of 104-wk daily oral exposure followed by 1 wk of clearance, the model underestimated the PFOS concentrations in the rat plasma or liver at week 14, whereas it overestimated at week 105 (Figure 5).At the early stage (week 14) following PFOS exposure, Butenhoff et al. 63 demonstrated that rats were in a healthy physiological state, and thus the gastrointestinal absorption capacity might be stronger than predicted.At 105 wk after exposure, Butenhoff et al. 63  ) and verification (Australia, [69][70][71] Canada, [66][67][68] and United States 27 ) results for human model of PFOS.It is assumed that all the estimated daily intake entered into the human body through oral route.The shaded area represents a predicted-to-observed ratio between 0.5:1 and 2:1.The corresponding numeric data are shown in Excel Table S5.Note: PFOS, perfluorooctane sulfonate; US, United States.
Environmental Health Perspectives 127012-10 131(12) December 2023 nephritis.In addition, Butenhoff et al. 63 observed that the decrease in food intake for the rat at later exposure stage also contributed to the overestimation.The model successfully predicted the concentration of PFOA in the rat plasma after repeated nasal exposure, whereas it overestimated the time to reach equilibrium (Figure S22).This overestimation may be due to the uncertainty in the values of physiological parameters such as alveolar ventilation and the body weight of the rats.
In comparison with the mice PBTK models, much lower Free, K bileC , and P Rest values calibrated in the human models well explains the much longer half-life of PFOS in humans reported previously (3.3-27 y). 27The difference between the route-specific external and internal exposures suggest that the relative contribution to the internal exposure levels is relevant to the exposure routes.For people who are exposed to PFAS through the dermal route frequently, considering the absorption efficiency of PFAS would generate a more accurate assessment of exposure risk.
Although the PBTK model here can be applied to humans, more experiment data are necessary to optimize the human PBTK model and identify the differences in the route-specific absorption processes between human and mice.For example, studies have shown that human skin is considered to be less permeable than mouse skin. 91The airways of mice are directly connected to the alveoli, whereas additional respiratory bronchioles are present in humans, 92 which may prolong the transition time in the respiratory airways.In addition, the characteristics of binding to proteins, such as the number of binding sites, dissociation constants, and a time-dependent description of the Free, need to be characterized in detail. 22,24n conclusion, this study explored the oral, nasal, and dermal exposure of PFAS simultaneously and clarified the routespecific pharmacokinetics of PFAS in mice.Based on the experimental results, route-specific PBTK models were established for, to our knowledge, the first time to predict the concentrations of PFAS in plasma and tissues following single or joint exposure as well as long-term repeated-dose scenarios.Oral exposure exhibited the highest bioavailability, likely due to a quick and effective absorption.However, dermal and nasal exposures exhibited lower PFAS bioavailability, perhaps due to a distinct hysteresis effect in dermal exposure and clearance via exhalation following nasal exposure.The model was capable of predicting long-term repeated exposure in mice, which was further extrapolated to humans.For people with high nasal and dermal exposure to PFAS, such as firefighters and infants, the contribution of nasal and dermal exposure cannot be excluded from the perspective of long-term exposure.The results improve our understanding of the impacts of exposure routes on ADME kinetics of PFAS in rodents and humans.More studies on nasal and dermal exposure to PFAS are urgently needed to reduce parameter uncertainty.The model would be improved if a full description of the variance in the parameters over time is incorporated.Finally, it is necessary to consider the influence mechanisms of these exposure routes, such as the parameterization related to transmembrane transport, enterohepatic circulation, and protein binding.

Figure 2 .
Figure 2.An overall evaluation of goodness of model fit for PFAS between the model-predicted and the observed concentration following 1 mg=kg oral (A), nasal (B), dermal (C), and IV (D) exposure.The solid diagonal line represents equal observed and predicted values.The dotted line indicates a predicted-toobserved ratio of >2:1 or <0:5:1.Linear fitting equation and R 2 were estimated through the linear regression model.The corresponding numeric data was shown in Excel TableS1.Note: IV, intravenous injection.

Figure 3 .
Figure 3. Simulated concentration-time curves of PFHxS (A), PFOS (B), and PFOA (C) in the plasma of mice following 1 mg=kg single oral, nasal, dermal, and IV exposure, and the time to peak (T max ) (D), predicted maximum concentration (C max ) (E), elimination half-life (T 1=2 ) (F), and bioavailability (G) of PFAS.Linear fitting between relative C max (the ratio of oral or dermal or nasal C max to IV C max ) and bioavailability (H); the dotted line represents the 95% confidence interval.T 1=2 are presented as the mean value of these three routes.Error bar represents the standard deviation.Linear fitting R 2 and p-value were estimated through the linear regression model.The corresponding numeric data are shown in Excel TableS2.Note: IV, intravenous injection; PFHxS, perfluorohexane sulfonate; PFOA, perfluorooctanoic acid; PFOS, perfluorooctane sulfonate.

Figure 4 .
Figure 4.The predicted results of combined model via multiple routes (lines).The triangles represent the concentrations (mean ± SD) of PFHxS, PFOS, and PFOA (letters A, B, and C, respectively) in plasma, liver, and kidney (numbers 1, 2, and 3, respectively) on day 2, 16, and 32 following the joint exposure (1 mg=kg for each route) (n = 4).The squares represent the sum of the concentrations in mice following a single exposure of 1 mg=kg via oral, nasal, or dermal exposure.The corresponding numeric data are shown in Excel TableS3.Note: Con, concentration; PFHxS, perfluorohexane sulfonate; PFOA, perfluorooctanoic acid; PFOS, perfluorooctane sulfonate; SD, standard deviation.

Figure 5 .
Figure 5. Application of the oral PBTK model to predict PFOS concentrations in plasma and liver in rat that were exposed every day for 104 wk, followed by 1 wk of clearance.Comparisons of the modeled vs. observed PFOS concentrations (mean ± SD) in rat plasma and liver dosed at 0.024 (A1 and A2), 0.098 (B1 and B2), 0.242 (C1 and C2) and 0:984 mg=kg=d (D1 and D2), respectively (n = 5).The experimental data (dots) are from Butenhoff et al.63The shaded area represents a predicted-to-observed ratio between 0.5:1 and 2:1.The corresponding numeric data are shown in Excel TableS4.Note: Con, concentration.

Figure 6 .
Figure 6.Model calibration (Norway65 ) and verification (Australia,[69][70][71] Canada,[66][67][68] and United States27 ) results for human model of PFOS.It is assumed that all the estimated daily intake entered into the human body through oral route.The shaded area represents a predicted-to-observed ratio between 0.5:1 and 2:1.The corresponding numeric data are shown in Excel TableS5.Note: PFOS, perfluorooctane sulfonate; US, United States.