Perfluorooctanoic acid (PFOA) is a fully fluorinated substance largely used in the manufacture of consumer products including food wrappers, non-stick cookware, and coatings (Bartell et al., 2010). PFOA is also known as a “forever chemical” due to its persistence both in the environment as it breaks down very slowly and potential to bioaccumulate in the human body (Rayne and Forest, 2009; Bartell et al., 2010; Shin et al., 2011). At a mechanistic level, PFOA has been implicated in increasing gene expression involved in cholesterol biosynthesis and metabolism (Fletcher et al., 2013; Louisse et al., 2016). Additionally, recent epidemiological studies identified critical effects on the immune system (Bampidis et al., 2020), which may result from a dysregulated cytokine/chemokine response and impaired neutralising antibody response (Lee et al., 2017). Current EU regulation (Article 25, EC 2020/2184) provides that member states shall take the necessary measures to ensure safe levels of a sum of per- and polyfluoroalkyl substances in water intended for human consumption by 2026. Given these effects, a thorough risk assessment of PFOA, including the derivation of an accurate health-based guidance value remains an important need.

Momentum for a paradigm shift in chemical risk assessment is growing (Hartung and Tsatsakis, 2021). Reduction in animal testing and its replacement with non-animal testing procedures received considerable impetus following the enforcement of the European Union (EU) Cosmetics Regulation (EC 1223/2009) in 2013. This directive imposed a full marketing ban in Europe for cosmetic products and ingredients tested on animals anywhere in the world. Other directives such as EC 2010/63/EU (amended by Article 6 of Regulation 2019/2010) and REACH (EC 1907/2006) also strengthen political pressure for progress towards replacing animal testing in science and the use of alternative test methods for hazard assessment of chemicals.

Alternative test methods that seek to minimise animal testing, known as New Approach Methodologies (NAMs), include in silico methods and their integration with in vitro bioassays that use human cells for determining a molecular point of departure (PoD). In the context of New Generation Risk Assessment (NGRA) (Johnson et al., 2022; Barutcu et al., 2023), the application of “omics” technologies to derive PoDs has been growing rapidly (Farmahin et al., 2017; Shah et al., 2021; Alcaraz et al., 2022; Beal et al., 2022; Jin et al., 2022; Kuo et al., 2022; Gou et al., 2023; Malinowska et al., 2023; Min et al., 2023; Ramirez-Hincapie et al., 2023; Reardon et al., 2023; Song et al., 2023). PoDs are used as the basis for the derivation of safe levels of human exposure known as Health Based Guidance Values (HBGV) (Ingenbleek L et al., 2021).

The benchmark dose (BMD) approach fits dose–response models to a complete dose–response dataset to identify the dose that corresponds to the lowest biologically relevant change in the response (i.e., benchmark response or BMR). The BMD is reported with a credible interval, e.g., the lower and upper confidence limits of the benchmark dose (BMDL and BMDU) for a selected observed level of effect (More et al., 2022).

The BMD is increasingly preferred by regulatory agencies, but its use is often limited by test design, which in turn are limited by current OECD guidelines for in vivo assays (Bokkers and Slob, 2005; Bokkers and Slob, 2007; European Food Safety, 2009; Hardy et al., 2017).

In vitro concentration-response data must be converted to in vivo dose-responses for use in human safety testing of chemicals. This process is known as quantitative in vitro to in vivo extrapolation (QIVIVE) (Yoon et al., 2012; Bale et al., 2014; Yoon et al., 2015). QIVIVE increasingly involves the application of physiologically based kinetic (PBK) modelling. In most previous QIVIVE studies, all PBK parameters other than input dose or exposure were fixed at central values. The discrepancy between a target in vivo dose, predicted by the PBK model, and a given in vitro concentration was minimised by an optimisation routine. The concentration corresponding to the target in vitro concentration was considered a surrogate for the in vivo dose. These studies, however, did not account for structural uncertainty in the PBK model, nor parameter value uncertainty (Louisse et al., 2010; Louisse et al., 2012; Strikwold et al., 2013; Louisse et al., 2016; Boonpawa et al., 2017a; Boonpawa et al., 2017b; Li et al., 2017; Punt et al., 2017; Strikwold et al., 2017; Adam et al., 2019; Zhao et al., 2019; Shi et al., 2020; Zhang et al., 2020). It is also well known that the amount of biological mechanistic detail described in a PBK model can have a bearing on model output (Rowland et al., 2017; McNally et al., 2018; Loizou et al., 2021).

Another limitation of earlier QIVIVE studies is the use of applied or nominal in vitro concentrations. This means that no consideration is made of the fate and distribution of the chemical in the in vitro test system. This could be a significant omission because the concentrations of chemical that are bioavailable to interact with the cells could be substantially lower than the nominal concentration (Tanneberger et al., 2010; Armitage et al., 2014; Groothuis et al., 2015; Kramer et al., 2015; Armitage et al., 2021; Proença et al., 2021). Chemicals that are highly lipophilic are particularly prone to interacting with constituents of the reaction medium and reaction vessel geometry and composition (Armitage et al., 2014; Proença et al., 2019; Armitage et al., 2021; Proença et al., 2021). Currently, there are several in silico tools that calculate the proportion of the nominal concentration that is free and presumably active (e.g., the in vitro mass balance Armitage model), and therefore available in the reaction medium to be taken up by cells (Tanneberger et al., 2010; Armitage et al., 2014; Comenges et al., 2017; Fischer et al., 2017; Fisher et al., 2019; Armitage et al., 2021). The Armitage model can estimate: (i) chemical fate and transport; (ii) cell partitioning; (iii) sorption to the reaction vessel; and (iv) considers the experimental set up (size of well-plate). Ultimately, it can simulate the active, free concentration of a chemical in the cell medium, which can be assumed to be equal to the organ or tissue efferent concentration, that is, the concentration exiting the organ or tissue and entering the venous blood compartment described in the PBK model. The Armitage model can also be used on its own to design and interpret in vitro experiments, and in combination with PBK models to perform in vitro to in vivo extrapolation (Comenges et al., 2017; Bell et al., 2018; Armitage et al., 2021). It is important to mention that the mass balance model has its own limitations and although it accounts for cellular uptake and excretion, it does not consider metabolic activity or the activity of specific cell transporters.

Here, we calculate a HBGV from the in vitro benchmark concentration modelling of previously published gene expression changes in human primary liver cell spheroids mediated by PFOA (Rowan-Carroll et al., 2021). The in vitro transcriptomics data were translated to in vivo dose-responses with a PBK model for PFOA (Worley et al., 2017; Loizou et al., 2021). A TDI (ng/kg BW/day) was calculated from the in vivo dose-responses and compared to the TDI values recommended by the European Food Safety Authority (EFSA) in the Scientific Opinion on the Risk to human health related to the presence of perfluoroalkyl substances in food for effects on the immune system (Bampidis et al., 2020) and previously for increases in serum cholesterol (Knutsen et al., 2018). Our work was achieved using a QIVIVE workflow described previously for ethylene glycol monoethyl ether (EGME) (McNally et al., 2018), PFOA (Loizou et al., 2021) and bisphenol A (Loizou et al., 2022) that accounts for PBK model structure uncertainty, parameter value uncertainty, and the calculation of free concentration of chemical in vitro by applying a rigorous statistical (probabilistic) framework to accommodate uncertainties. Overall, the purpose of this study, however, was not to propose an updated animal-free risk assessment for PFOA, since we recognise that much work is still needed to demonstrate in vitro to in vivo concordance for systemic, chronic exposures to environmental xenobiotics. Such a demonstration is beyond the scope of this study. Instead, this study serves to further demonstrate and build confidence in the utility of “omics”-derived molecular PoD measurements in combination with the multi-step in silico QIVIVE workflow in anticipation of the acceptance of in vitro concentration-response data in chemical risk assessment.

The step-by-step workflow used in this study is depicted in Figure 1. Each of the following subsections are mapped according to the five steps in our workflow.

In the foreseeable future, risk assessments are increasingly likely to incorporate animal-free testing methods and gain acceptance. There are extensive and ongoing efforts worldwide to increase the uptake of NAMs, including to gain their acceptance by regulatory agencies. Next-Generation Risk Assessment (NGRA)-aligned frameworks embrace the notion of “protecting not predicting”, which has the goal of protecting humans from harmful exposure levels to a plethora of substances found in the environment, food, drinking water, and other sources. This is where this piece of work is positioned—aligned with NGRA principles, considering the uncertainties of each method integrated into the proposed in silico workflow, and demonstrating the utility of well-documented NAMs within a case study.

Our work builds on the gene expression data generated by (Rowan-Carroll et al., 2021) upon exposure of liver spheroids to different concentrations of four PFAS, including PFOA, over a time course of 14 days. Our study also relied on a probabilistic approach that combines a calibrated PBK model (Loizou et al., 2021) and a population based QIVIVE algorithm for deriving exposure doses. There are several advantages regarding exposure or dose reconstruction provided by this probabilistic approach. First, by defining informative prior distributions around parameters, a deterministic model converts to a population model, which can account for inter-individual variability. Second, the probabilistic approach is appropriate for systems where tissue dose is not necessarily linearly related to external exposure. For instance, there is evidence demonstrating that interindividual variability of blood protein concentrations can affect the relationship between external exposure and internal doses of PFAS due to their particular toxicokinetic behaviour (Fischer et al., 2024). Because per- and polyfluoroalkyl substances tend to bind to serum proteins and nuclear receptors, their half-lives also tend to be increased; on the other hand, if a disease lowers the quantity of serum proteins, there will be more unbound PFAS circulating, hence a higher susceptibility to their adverse effects is also expected (Zaccherini and Bernardi, 2019; Nielsen et al., 2024; Rosato et al., 2024). Finally, this combination can extract population variability and multiple routes of exposure information integrated within pharmacokinetic data (McNally et al., 2012; McNally et al., 2018).

We propose a solution to the oftentimes arbitrary choice of a biologically relevant molecular PoD made from the observed perturbations of single genes by first determining the most sensitive pathways to the chemical exposure via a gene set enrichment analysis of functionally linked gene responses. To fulfil this approach, we derived molecular PoDs from the publicly available transcriptomics data using BMDExpress3, which in turn allowed for ranking genes according to their BMDL values. The question we sought to answer was straightforward; our aim was to demonstrate the utility of an in silico workflow designed to extrapolate in vitro concentration-response data related to the responses of the most sensitive and statistically relevant pathways activated by exposure to low, environmentally-relevant concentrations of PFOA to in vivo doses. The pathway enrichment analysis highlighted “fatty acid metabolism”, “MYC targets V1 and V2”, and “unfolded protein response” pathways as the most sensitive to PFOA exposure. ECH1 was the most sensitive gene among the 60 genes in the fatty acid metabolism pathway. The BMDL value for ECH1 (0.25 µM) was 70-fold lower than the BMDL for PPARα (17.33 µM), which is consistent with studies demonstrating the relation between PFOA exposure and activation of PPAR-related pathways (Kaiser et al., 2022).

PPAR-related pathways play an important role in lipid metabolism in humans; PPAR receptors family is composed by three subtypes—α, δ, and γ—and differ in their spectrum of activities. PPAR α is mainly expressed in the liver, but can also be found in the heart, bones, and muscles (Kamata et al., 2023). PPARα activation has been described by Yang et al. (2023) as an important molecular initiating event in PFOA/PFOS induced dysregulation of hepatic lipid metabolism both in humans and rodents; besides, the authors also demonstrated that hepatic ACOX1 in the peroxisomal β-oxidation pathway is also a critical molecular target associated with PFOA/PFOS dysregulated lipid metabolism and resultant hepatotoxicity. In this work, we identified ACOX1 as the second most sensitive gene, with a BMDL of 0.37 µM, and the gene set enrichment analysis recognised it as part of the fatty acid metabolism pathway.

BMDL concentrations lower than 1 µM activated responses in at least 9 genes in the fatty acid metabolism pathway, suggesting the onset of molecular events that could eventually lead to downstream phenotypic outcomes. Such low concentrations derived from gene expression data are anticipated to translate to conservative TDIs using our workflow to integrate transcriptomics concentration-response data with a calibrated PBK model. These results are also useful for the study of molecular key events of relevant adverse outcome pathways and for the question of how many genes need to be initially perturbed for the test system to be considered perturbed.

The TDI derived using our computational workflow was in the same range and in general slightly more conservative than the values derived by EFSA. Loizou and colleagues (2021) applied the same workflow, using in vitro concentration-response data from the ToxCast/Tox21 database. That study focused on concentration-response data related to estrogen, pregnane X and thyroid receptor binding, associated with breast cancer, altered serum cholesterol levels, thyroid disease/decreased antibody titre. As shown in Table 2, the TDI derived using transcriptomics data, which interrogates a broader range of molecular responses, was more conservative than those values derived using cell-based high-throughput screening in vitro data.

An aim of this study was to demonstrate the utility of a workflow designed to quantify uncertainty within the different stages of data analysis and processing and to account for human inter-individual variability. The health-based guidance value determined in this study is not proposed for the future regulation of PFOA. Further work is required to provide the confidence in the use of in vitro to in vivo NAMs required by regulatory agencies. Furthermore, this study has a number of limitations. Firstly, the absence of parameters for liver spheroids in the current version of the mass balance model used in this work. Considering that the liver spheroids were composed of primary human hepatocytes and Kupfer cells, we have addressed this limitation by using the parameters for HepG2 cells according to standard parameters of an in vitro ToxCast assay for PFOA. This limitation is linked to the uncertainty around the freely dissolved concentrations of PFOA, as our approach considered a common cell culture and not a spheroid-based biological test system.

Another important limitation to be discussed is how different mathematical (in silico) models address the estimation of the freely dissolved fraction of a chemical within an in vitro test system. Different methods with different parameters have been proposed and to determine the most adequate method to estimate the mass balance of the different chemical classes is anything but trivial (Proença et al., 2021). Moreover, determining the applicability domain of such models is also challenging. In our study, we have tested both the updated version of the in vitro mass balance model (Armitage et al., 2021) and the virtual cell-based assay (VCBA) model (Proença et al., 2019). We decided to proceed with the updated version of Armitage’s model (Armitage et al., 2021) as it had the capability to simulate a ToxCast assay with larger microplate settings and handle ionizable organic compounds such as PFOA. It is noteworthy that the correction of the nominal dose to derive an available fraction is an important step of the workflow. We have reviewed some of the available models, their data requirements, and limitations in the introduction, but a broader review can be found elsewhere (Proença et al., 2021). As with all models there is uncertainty associated with predictions, and different models may be better for certain classes of chemicals. A particular model may have both uncertainty and a bias in predictions of the freely dissolved fraction for certain chemical classes (Heringa et al., 2004; Teeguarden and Barton, 2004; Treijtel et al., 2004; Treijtel et al., 2005; Madureira et al., 2014; Kwon et al., 2020) or may be more general (Gulden and Seibert, 2003; Armitage et al., 2014; Stadnicka-Michalak et al., 2014; Comenges et al., 2017; Fischer et al., 2017; Fisher et al., 2019). In this work there has been no explicit consideration of uncertainty and a potential bias in the calculations for the free fraction. However, we note that model uncertainty, arising from a multitude of (unspecified) sources, is explicitly considered in the QIVIVE through the ABC algorithm, therefore we consider this uncertainty to be already addressed within the workflow. Any uncertainty arising from the correction of the nominal dose is in any case trivially small relative to the case of an uncorrected nominal concentration used within calculations.

In the present work, despite the lack of specific parameters for liver spheroids in the current version of Armitage’s model (Armitage et al., 2021), the freely dissolved concentrations of PFOA was estimated to be 2.9% of the nominal concentrations. Previously (Loizou et al., 2021), tried to predict the freely dissolved concentrations of PFOA as 2.3% of the nominal concentrations based on its LogP. It was not based on any specific model, but solely relying on the estimations for benzo [a]pyrene (which has a similar LogP to PFOA) calculated and published by (Proença et al., 2019).

Another important limitation is that our in vivo benchmark calculation accounted for the data of just one exposure time point of the original study (day 10), with a visual comparison with the other time points. The transcriptomics data is functional; however, it is a snapshot of the exposure to a chemical at that time point. As depicted in Figure 2, the number of perturbed genes in the fatty acid metabolism pathway increases and reaches its maximum at day 10, and decreases on day 14, as the gene expression pattern returns towards baseline.

Kaiser et al. (2022) mapped PFOA to AOPs using the AOP-helpFinder tool (http://aop-helpfinder.u-paris-sciences.fr/index.php) (Kaiser et al., 2022). The concept of an AOP is to relate the knowledge about a molecular initiating event (MIE) with an adverse outcome (AO) through several key events (KE), connected by key events relationships, occurring at different biological levels (Ankley et al., 2010). Their study highlighted the relation between PFOA and AOPs within AOP-wiki using epidemiological, in vitro, and in vivo published studies. Their results also reinforced the link between PFOA and PPAR-related pathways, which is consistent with the results of the pathway enrichment analysis we performed.

Recently, Chen et al. (2022) investigated the links between perfluorooctane sulfonate (PFOS) and molecular pathways using a similar approach and highlighted the response that this perfluorinated chemical induces to cholesterol- and fatty acid metabolism-related pathways (Chen et al., 2022). However, despite using a PBK model, the authors did not estimate the free concentrations of PFOS and only accounted for the extra quantity of albumin naturally produced and secreted to the media by the liver spheroids. Also, to derive their human equivalent doses they applied an uncertainty factor of 3 for toxicodynamic and 10 for toxicokinetic differences, which resulted in an overall UF of 30. Our use of a probabilistic PBK model accounted for the toxicokinetic UF. Therefore, the application of an UF of 3.16 to account for the toxicodynamic component was required in our study. The derivation of a similar HBGV to that of (Chen et al., 2022), where those authors did not estimate the free in vitro concentration of PFOS, may be explained by their use of an overall UF of 30. However, in both studies, the use of transcriptomics data coupled with PBK and QIVIVE allows for the identification of gene-pathway-level effects triggered by exposure to a given stressor, which can provide vital information to support chemical risk assessment. Another important point in common to both studies is the potential to identify different, and still comparable, sensitive pathways and diseases that would not be detected without the use of “omics data.”

In conclusion, we have developed an in silico workflow that integrates transcriptomics data with a calibrated PBK model and quantitative in vitro to in vivo extrapolation algorithm to derive a tolerable daily intake for perfluorooctanoic acid. Our approach incorporated an estimation of the in vitro free concentrations of PFOA and a pathway enrichment analysis that allowed us to identify the most sensitive molecular pathway related to the exposure to this “forever chemical”.