In the environment of the modern world, chemicals of anthropogenic and natural origin are ubiquitous. The diversity of chemicals that risk assessors must appraise is large. For example, in the food sector this includes anthropogenic contaminants such as pesticides, biocides, food and feed additives, pharmaceuticals, air pollutants, persistent organic pollutants, heavy metals, perfluoroalkyl substances, brominated flame retardants, dioxins, and those of natural origin (marine biotoxins, mycotoxins, etc.) to name but a few. In this context, human risk assessment of chemicals aims to quantify exposures in human subpopulations from relevant sources and exposure routes (exposure assessment), to identify and characterize adverse effects and determine safe levels (hazard identification and characterization) as well as quantify risks associated with such exposures (risk characterization) (Ingenbleek et al., 2021).

Hazard identification and hazard characterization requires an understanding of what the human body does to the chemical, known as “toxicokinetics (TK)” and what the chemical does to the body, known as “toxicodynamics (TD).” Over the last century, the TD dimension and the derivation of acute and chronic safe levels of exposure for human health has been mostly addressed using animal toxicological studies in test species (rat, mouse, rabbit, and dog). Such studies allow the identification of an apical endpoint which is an observable outcome in the whole animal including clinical signs or pathologies indicative of a disease state that can result from exposure to a toxicant. For a given chemical, the key apical endpoint is identified based on the role of the endpoint in causing adversity and the dose-response relationship observed in the critical toxicity study in a test species (rat, mouse, rabbit, and dog). This principle is based on the requirement that all relevant adverse effects should be considered in order to achieve a sufficient level of protection. As a consequence, all reliable toxicological studies and relevant apical effects are considered to identify the highest dose that does not produce any statistically significant increase in the incidence of an adverse effect. Such a dose or concentration, known as the reference point (RP) or point of departure (PoD), is defined as the point on a toxicological dose–response curve established from experimental data that corresponds to an estimated no-observed-adverse-effect level (NOAEL) or low effect level. PoDs are used as the basis for the derivation of safe levels of human exposure known as health-based guidance values (HBGVs) (Ingenbleeks et al., 2021).

The most common RPs or PoDs are the NOAEL and the benchmark dose (BMD). The NOAEL approach uses statistical methods to identify the tested dose with no significant effect compared to the control group. The BMD approach fits a dose–response model(s) to a complete dose–response dataset to identify the benchmark dose lower confidence limit (BMDL) for a selected observed level of effect, the benchmark response (BMR) (e.g., a 10% response). The BMD is increasingly preferred by regulatory agencies, but its use is often limited by test design (Bokkers and Slob, 2005; Bokkers and Slob, 2007; European Food Safety, 2009; EFSA Scientific Committee et al., 2017). From the RP or PoD, HBGVs are derived most often by applying a default uncertainty factor (UF) of 100 allowing for interspecies (10-fold) and inter-individual differences (10-fold) when no TK and TD information is available. Alternatively, when chemical-specific or pathway-specific TK or TD information exists, more informed UF values can be used. For chemicals with a threshold for toxicity (non-genotoxic) chemicals, HBGVs in the food and feed safety area include the acceptable daily intake (ADI) for food additives, feed additives and pesticides, tolerable daily intake (TDI) for contaminants, upper limits (UL) for vitamins and minerals, and, for acute effects, the acute reference dose (ARfD) (Sheffer, 2009; EFSA Scientific Committee, 2012).

Since PoDs and HBGVs mostly rely on toxicological studies in test species, the international scientific community has been involved in considerable research and validation efforts to reduce animal testing and provide alternative-to-animal testing methods (i.e., in vitro and in silico) known as new approach methods (NAMs). The development of an alternative to animal human safety testing strategy for chemicals has been described as being akin to seeking the Holy Grail (Louisse et al., 2016). Efforts toward achieving this goal have been ongoing since publication of the United States National Research Council (NRC) report “Toxicity Testing in the 21st Century: A Vision and a Strategy” (NRC, 2007). In Europe, the development of such a strategy received impetus following the full marketing ban enforced under the EU Cosmetics Regulation (EC 1223/2009) in 2013 for cosmetic products and ingredients tested on animals anywhere in the world (Coecke et al., 2013).

In practice, alternative to animal methods invariably refer to an in vitro bioassay-based strategy that ideally uses human cell lines primarily for the determination of a RP. Importantly, in vitro concentration–response data must be converted to in vivo dose responses in order to be used in human safety testing of chemicals. This activity is known as quantitative in vitro to in vivo extrapolation (QIVIVE) (Bale et al., 2014; Hartung, 2017). Examples of QIVIVE increasingly involve the application of physiologically based kinetic (PBK) modeling–based reverse dosimetry for the translation of in vitro to in vivo responses and the derivation of in vivo BMDs (Adam et al., 2019; Boonpawa et al., 2017; Li et al., 2017; Louisse et al., 2016; Louisse et al., 2010; Louisse et al., 2012; Punt et al., 2017; Shi et al., 2020; Strikwold et al., 2017a; Strikwold et al., 2013; Strikwold et al., 2017b; Zhang et al., 2020; Zhao et al., 2019). Within this approach, all parameters, other than input dose or exposure, are held fixed at central values. An optimization routine is implemented in order to minimize the discrepancy between a target in vivo concentration, predicted by the PBK model and a given in vitro concentration. The dose concentration which corresponds to the target in vitro concentration is considered to be a surrogate for the in vivo concentration. However, these studies did not account for PBK model structure or parameter uncertainty. Therefore, an algorithm, described in detail previously (McNally et al., 2018), was developed to extrapolate in vitro concentration–response to in vivo dose–response relationships while applying a rigorous statistical framework for accommodating uncertainty in both the parameters of the PBK model, the quality of fit of the model to measure biological monitoring data, and a consideration of how this affects an in vivo dose–response relationship in the context of QIVIVE (Judson et al., 2011). This is important since the level of detail (fidelity) captured in the model could have a bearing on model output (Rowland et al., 2017). Understanding and quantifying the level of uncertainty in each step of a chemical safety assessment with NAMs is important for the development of confidence in this approach (Berggren et al., 2017). The workflow uses global sensitivity analysis (GSA), PBK modeling, approximate Bayesian computation (ABC), and Markov Chain Monte Carlo simulation to convert in vitro concentration–response data to in vivo dose–response data (McNally et al., 2018).

There are a number of advantages with regard to exposure or dose reconstruction provided by this probabilistic approach. First, defining informative prior distributions around parameters converts a deterministic model 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. 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).

In this report, we tested the workflow using HTS in vitro concentration–response data for perfluorooctanoic acid (PFOA). The data were obtained from the ToxCast/Tox21 database on the US EPA Chemistry Dashboard (Williams et al., 2017) and translated to in vivo dose responses with a PBK model for PFOA (Worley et al., 2017). An in vivo BMD for PFOA intake (ng/kg BW/day) was calculated from the in vivo dose responses and compared to the intake also derived from a BMD used 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 (EFSA Contam Panel, 2020) and previously for increases in serum cholesterol (EFSA Contam Panel et al., 2018).

PFOA is a synthetic chemical comprised of a fully fluorinated eight carbon chain with a carboxylic acid functional group. It was used in the manufacture of many consumer products including fast food wrappers, non-stick cookware, a stain-resistant coating used on carpets and other fabrics (Bartell et al., 2010). PFOA is both hydrophobic and lipophobic, does not break down in the environment, contaminates drinking water sources, and accumulates in food chains (Bartell et al., 2010; Rayne and Forest, 2009; Shin et al., 2011). It is not metabolized in the human body, and the half-life is estimated to be between four and five years (Emmett et al., 2006; Steenland et al., 2010). Epidemiological studies support a possible association with liver, pancreatic, and testicular cancer (Lau et al., 2007; Barry et al., 2013; Vieira et al., 2013) and breast cancer (Bonefeld-Jorgensen et al., 2011).

The purpose of this study, however, was not to propose an animal-free risk assessment for PFOA, since it is recognized that much work is still needed to demonstrate in vitro to in vivo concordance for systemic, chronic exposures to environmental xenobiotics. Instead, this study serves to further demonstrate the utility of the algorithm in anticipation of the acceptance of in vitro concentration–response data in chemical risk assessment.

Nevertheless, a workflow was followed which would be similar to most risk assessment approaches. For example, epidemiological studies reported that PFOA is a suspected endocrine disruptor (Pierozan et al., 2018) associated with a risk of breast cancer (Bonefeld-Jorgensen et al., 2011). Other studies observed an increase in cellular triglyceride levels and in vivo expression of genes involved in cholesterol metabolism (Fletcher et al., 2013) and gene sets related to “PPAR signaling,” “lipid metabolism,” “fatty acid beta oxidation,” and “tRNA amino-acylation” which are related to “cholesterol biosynthesis” which may be associated with hepatic steatosis (Louisse et al., 2020). Functional thyroid disease was observed in a large cohort of people exposed to PFOA in drinking water contaminated from a mid-Ohio River Valley chemical plant (Winquist and Steenland, 2014) where reduced expression of parathyroid hormone 2 receptor which may increase risk for conditions related to parathyroid hormone signaling was also observed (Fletcher et al., 2013; Galloway et al., 2015). Finally, recent epidemiological studies found associations with effects on the immune system (reduced antibody titers) (EFSA Contam Panel, 2020) which were hypothesized to result from a dysregulated cytokine/chemokine response and impaired neutralizing antibody response (Lee et al., 2017). Therefore, in vitro datasets were selected as measured responses that could be related to the mechanistic understanding related to PFOA toxicity and observations in human exposure studies, such as, estrogen receptor binding activation associated with breast cancer, pregnane X receptor binding associated with hepatic steatosis, thyroid hormone receptor α binding leading to thyroid disease, and immunotoxicity (pro-inflammation from cytokines).

In this work, we applied the workflow of McNally et al. (2018) in order to account for parameter and model uncertainty within a computationally efficient framework. This the first time the workflow has been applied to a human PBK model and human in vitro cell line data. A technical discussion of the framework and in particular a justification for the inclusion of model uncertainty are covered in detail in McNally et al. (2018) and therefore is not repeated here. However, a thorough comparison between the “reverse dosimetry” approach adopted in the workflow and the iterative forward dosimetry “dose matching” approach to QIVIVE that is frequently adopted in the literature is of merit. In our approach, the uncertainties and variabilities in model parameters have been specified through probability distributions, the sensitivity of model output to uncertainties in model inputs has been tested through sensitivity analysis, uncertainty has been refined (although not completely eliminated) through calibration, and finally parameter uncertainties and model uncertainty have been accounted for during QIVIVE. In contrast, the “dose matching” approach does not account for uncertainty at any stage of the modeling process. Results are inevitably sensitive to the default parameters assumed by researchers, but there is currently no mechanism for quantifying this sensitivity.

The PBK model for PFOA was written in MCSim syntax, compiled and subsequently run using the generated executable. MCSim was an appropriate platform given the long simulation period (up to 350,000 h which took around 10 h on the Microsoft Azure IaaS) and the requirement for hundreds of thousands of model evaluations (covering the processes of development and testing, GSA, calibration, and QIVIVE). However, while MCSim is particularly well suited to computationally intensive computations such as MCMC which can be executed at the command line, the command line interface was not well suited to the overall workflow adopted for this work. Interaction with the model executable was undertaken in three distinct ways. RVis was used for development and testing of the model and for GSA. RVis acts as a user-friendly front end to an MCSim executable and supplies input, runs, processes and reports model outputs from the executable. In addition to providing a user-friendly front end for interacting with the model, RVis parallelizes the runs in batch-run operations over the available cores on a personal computer, such that the run time for a computationally expensive technique such as eFAST (GSA), which requires tens of thousands of runs, may be substantially reduced. Calibration was performed using MCMC (a technique for which MCSim is particularly well suited and run through the command line). Model output from MCMC was subsequently processed using R scripts. QIVIVE was also performed using R scripts which supplied inputs to, ran, and processed outputs from the executable. The computation for this process was exported to Azure as described in Methods. At present the full workflow requires significant expertise in PBK modeling, statistics, and programing and is not widely accessible. A focus of current work is to develop a user-friendly module within RVis for QIVIVE, which would substantially reduce the entry barrier to the technique.

This study has demonstrated that the availability of freely accessible in vitro concentration–response data for environmental pollutants from the Tox21 and ToxCast high-throughput in vitro screening programs could be valuable and effective in chemical risk assessment. The credible intervals for the in vivo BMDL values for Intake for serum cholesterol (ATG_PXRE_CIS_up) with and without the application of a CSAF or default UF encompass the BMDL for Intake derived by the EFSA. However, only the in vivo BMD values for Intake for decreased antibody titer (BSK_3C_uPAR_down) with application of the CSAF encompassed the value derived and used by the EFSA. These results suggest that the in vitro concentration–response data could be used to predict in vivo dose response successfully.

At present we address each in vitro concentration in turn and estimate corresponding in vivo concentrations, accounting for model and parameter value uncertainty with a distribution of in vivo concentrations resulting from each in vitro concentration. From this data, we derive a dose–response curve from central estimates and calculate a BMD. This BMD can be used in the risk assessment process, with standard uncertainty factors applied as necessary. With a minor change to methodology, our approach could be used to derive CSAFs. The principle change would be to derive the full sequence of in vivo concentrations, corresponding to a given set of uncertain parameters and the full sequence of in vitro target concentrations, within a single step. The dose–response data could be interpreted as corresponding to a given individual in a population, and a BMD estimated from the data. Through repeating and generating a sample of BMDs, uncertainty factors may be derived. In application, there are considerable technical challenges in achieving a reasonable acceptance rate for proposals and computational efficiency. This is a priority area of research going forward.

However, there are a number of caveats that must be considered. For example, we applied a simple approach to mitigate the use of nominal, applied concentrations. Estimation of the available free concentration of the test compound is a function of serum protein and lipid composition of the media, and it is foreseen that more accurate estimates of bio-actively available in vitro concentrations will lead to different and more accurate extrapolated in vivo concentrations. In addition, the majority of assay protocols appear to be standardized where the in vitro concentrations and spacing are identical for all datasets, spanning three orders of magnitude, ranging from 41.41 to 41,407 μg/L. This may not be ideal for the derivation of a BMDL for which the concentration spacing might need to be altered to more faithfully capture the changes in response. More generally, the assays used in this study involved binding of unmetabolized chemical to a receptor which would initiate a cascade of interactions assumed to cause the perturbation of various biochemical pathways that precede overt toxicity. Therefore, these assays are limited to reactions that do not require metabolism of chemical to an active moiety. Indeed, the challenge of predicting a priori the rate at which a potential toxicant is metabolized in vivo still remains unsolved (Wetmore et al., 2012).

It is also important to document known limitations of the assays. For instance, assays using fluorescent readouts can give unreliable results for compounds that are themselves fluorescent (e.g., azo dyes). With cell-based assays, simultaneous cytotoxicity measurements are usually needed because cytotoxicity can confound the target-specific readout. Despite the known limitations of in vitro assays, there are plenty of examples of increased utility when combined with PBK models. For example, if in vivo clearance is included, approximately two orders of magnitude variation in biological potency could be captured when predicted directly using HTS in vitro measurements (Knudsen et al., 2015).

While our ABC algorithm is more appropriate for prioritized chemicals for which a thorough risk assessment is justified, the ability to quantify model parameter, structure, and data uncertainty is of paramount importance for the development of confidence in model use. The ability to do this within an intuitive, freely available software tool would make this approach accessible to a wider user base. This is an important step forward to implement the use of NAMs in chemical risk assessment since QIVIVE and reverse dosimetry have been described as the critical “endgame” in the workflow of predictive toxicology (Knudsen et al., 2015). QIVIVE is essential in order to transition away from animal model–based toxicology to entirely in vitro/in silico-based toxicological science (Knudsen et al., 2015). It is recommended to further test the current probabilistic workflow allowing the derivation of BMDL from the integration of data from in vitro assays, PBK modeling, and QIVIVE approaches, through relevant case studies for chemicals of relevance to the food safety area and chemical risk assessment in general.