Advancing understanding of human variability through toxicokinetic modeling, in vitro-in vivo extrapolation, and new approach methodologies

Abstract

The merging of physiology and toxicokinetics, or pharmacokinetics, with computational modeling to characterize dosimetry has led to major advances for both the chemical and pharmaceutical research arenas. Driven by the mutual need to estimate internal exposures where in vivo data generation was simply not possible, the application of toxicokinetic modeling has grown exponentially in the past 30 years. In toxicology the need has been the derivation of quantitative estimates of toxicokinetic and toxicodynamic variability to evaluate the suitability of the tenfold uncertainty factor employed in risk assessment decision-making. Consideration of a host of physiologic, ontogenetic, genetic, and exposure factors are all required for comprehensive characterization. Fortunately, the underlying framework of physiologically based toxicokinetic models can accommodate these inputs, in addition to being amenable to capturing time-varying dynamics. Meanwhile, international interest in advancing new approach methodologies has fueled the generation of in vitro toxicity and toxicokinetic data that can be applied in in vitro-in vivo extrapolation approaches to provide human-specific risk-based information for historically data-poor chemicals. This review will provide a brief introduction to the structure and evolution of toxicokinetic and physiologically based toxicokinetic models as they advanced to incorporate variability and a wide range of complex exposure scenarios. This will be followed by a state of the science update describing current and emerging experimental and modeling strategies for population and life-stage variability, including the increasing application of in vitro-in vivo extrapolation with physiologically based toxicokinetic models in pharmaceutical and chemical safety research. The review will conclude with case study examples demonstrating novel applications of physiologically based toxicokinetic modeling and an update on its applications for regulatory decision-making. Physiologically based toxicokinetic modeling provides a sound framework for variability evaluation in chemical risk assessment.

Introduction

Toxicokinetics (TK) is the study of xenobiotic fate in an organism, which requires the consideration of absorption, distribution, metabolism, and excretion (ADME) [1]. Alternatively, toxicodynamics (TD) describes the study of the effects a chemical may have on an organism, which in turn will require distinct considerations from TK [2]. Simply put, TK reflects what the body does to a chemical; TD reflects what the chemical does to the body.

A significant challenge worldwide in chemical or drug safety evaluations is understanding the extent and consequences of interindividual and population variability [3, 4]. Such variability can stem from differences in TK that lead to differing internal concentrations despite equivalent external exposure levels; or differences in TD, represented by differing susceptibilities to an adverse effect [5]. These research areas (i.e., TK and TD) are analogous to pharmacokinetics (PK) and pharmacodynamics (PD) as applied to pharmaceuticals, but with broader consideration of a wider range of exposures, dose–response, and effects or adverse outcomes that are important in toxicology during chemical risk assessment [6]. As such, the majority of the text, except when capturing historical PK concepts, will use the phrases TK and TD to describe TK/PK and TD/PD concepts. Regardless of which space is being considered, variability assessment requires distinct strategies for adequate consideration in human health risk assessment and has been an active area of research and discussion for well over 30 years [3,4,5, 7,8,9].

Of particular note is the vast amount of knowledge gained regarding differences that contribute to TK variability, which include physiology, genetics, ontogenetics, and exposomics as depicted in Fig. 1. Data from genomics technologies have shown that numerous polymorphisms exist across enzymes and biological targets, and how these vary by ethnicity and life-stage [10, 11]. A concerted effort fueled by the recognition of developmental changes in enzymes has provided a rich source of Phase I and II metabolic enzyme abundance data across life-stages [12,13,14,15]. The subsequent emergence of targeted proteomics as a tool to monitor enzyme and transporter abundance has further expanded the ability to evaluate differences in distribution and metabolism [16, 17]. Recent progress in understanding xenobiotic-microbiome interactions is poised to inform lines of inquiry regarding impacts on chemical biotransformation, fate, and potential interactions that will contribute to variability [18, 19]. Population variability characterization remains a highly active research area.

Fig. 1

Drivers of toxicokinetic (TK) variability. Physiological drivers include factors such as blood flows. Genetic factors include polymorphisms in enzymes such as CYP2D6, which results in varying rates of compound metabolism by poor metabolizers (P M), intermediate metabolizers (IM), extensive metabolizers (EM), and ultrarapid metabolizers (UM). Ontogenetic factors include changes in CYP expression over the course of pre and postnatal development. Adapted from “Fig. 1: Key life stages” from “Guidance on the risk assessment of substances present in food intended for infants below 16 weeks of age,” by EFSA Scientific Committee, used under CC BY-ND 4.0 [20]. Exposomic drivers include diet, lifestyle, environmental justice, and exposure/co-exposures

In parallel, the toxicology community has made significant strides in developing new approach methodologies (NAMs) to evaluate chemical safety. Increased consideration of alternative testing strategies that incorporate experimental human in vitro methods and in silico modeling approaches has sparked substantial investment in industry, academic, and regulatory sectors with the aim of incorporating such models and systems into hazard identification [21,22,23]. NAMs allow for the interrogation of mechanistic pathways or endpoint-specific effects that may provide insight into chemical modes of action. Moreover, NAMs can be adapted for high-throughput screening and can circumvent the need for interspecies extrapolation by using human-specific in vitro models [24, 25]. While not designed to provide a full replacement for in vivo studies, these NAMs are beginning to be applied to inform certain decision contexts (e.g., prioritization) and even waive in vivo studies given sufficient supporting evidence, as they continue to mature [26, 27].

To harness the advantages of NAMs for human health risk assessment, assay readouts must be translated into in vivo metrics. This can be achieved through in vitro-in vivo extrapolation (IVIVE). Physiologically based toxicokinetic (PBTK) frameworks are particularly well suited for IVIVE by bridging the gap between in vitro and in vivo TK both qualitatively and quantitatively [28]. These PBTK models allow for incorporation of in vitro data and in silico predictions to derive human equivalent doses [28]. Consideration of variability of these data allows for extrapolation of a range of in vivo equivalent doses across individuals rather than a single value, which shall more accurately reflect the distribution of population responses that occur from exposure to a chemical [29].

This review will describe major drivers that contribute to TK/TD variability and provide background on the structure and evolution of TK and PBTK models, designed to address a varied range of exposure scenarios across a wide range of data inputs. This will be followed by a state of the science update describing current and emerging modeling strategies and emerging experimental and modeling strategies that incorporate population and life-stage variability, including the increasing application of IVIVE with PBTK models in pharmaceutical and chemical safety research. The review will conclude with case study examples demonstrating novel applications of PBTK modeling that demonstrate the power and potential of these models to evaluate variability.

Drivers of human TK variability

In order to provide context for incorporation of variability into PBTK modeling, it is necessary to understand the key drivers of human TK variability. Factors that contribute to population TK variability can be divided across four different groups: physiological, ontogenetic, genetic, and exposomic. Table 1 provides a summary of these factors and considerations regarding their impacts on variability.

Table 1 Contributors to TK variability

Physiological factors

Differences in the underlying physiology involved in ADME can contribute to variations in TK. Although variability does indeed exist within and between individuals within the same life-stage, notable physiological differences are more pronounced when comparing adults (e.g., 20–50 years old) to children or elderly (e.g., > 60 years old), or groups with differing health conditions/status (e.g., heathy subject vs. one with liver disease, non-pregnancy vs. pregnancy) [30,31,32,33]. The impact of physiological differences is perhaps best conveyed through comparisons between young and elderly adult populations (e.g., 25-year-old vs. 65-years-old). The elderly generally show reduced hepatic metabolism and renal clearance due to a variety of age-related physiological changes, such as reduced liver blood flow, a reduction in hepatic cytochrome P450 (CYP) half-life, and a reduction in muscle mass and water content [30]. Both the fat-to-muscle ratio and water content can change the volume of distribution (V_d), the theoretical body volume in which a drug is distributed to reach the same concentration observed in the plasma, e.g., an increase in the former, or loss of the latter would increase the V_d for lipophilic compounds and decrease the V_d for hydrophilic compounds. Moreover, there is a reduced ability in the elderly to counter environmental exposures and stressors, as evidenced by an accumulation of lipofuscin, which has been linked to a drop in the ability to offset oxidative stress, e.g., a decline in brain antioxidant levels [34]. In line with this, brains in the elderly show an increased sensitivity to neurotoxicants, which may in part be caused by reduced activity of acetylcholinesterase, the enzyme that breaks down the neurotransmitter acetylcholine, and reduced functional reserves of the aging brain [30]. In addition, liver and kidney disease are more common in the elderly, further altering metabolism and clearance as discussed below [30].

Alternatively, infants and children are characterized by differing physiologic composition as their bodies evolve and develop into adults [35]. Infants up to three months of age have an elevated water: lipid ratio, generally resulting in an increase in the half-life and V_d of hydrophilic chemicals [31, 36]. Children up to six years of age have a larger brain to body weight ratio and greater blood flow to the central nervous system, thereby potentially increasing the distribution and half-life of a chemical in the brain [31, 37,38,39]. Other physiological differences include tissue volumes and blood flow, and immaturity of respiratory, renal, and hepatic systems in children [40]. Although more of an exposure factor, children also generally have increased inhalation and food intake in proportion to body weight. This is of particular concern due to children’s increased exposure to soil and house dust, which can be contaminated with toxicants [41]. There is also evidence for differences in the absorption of chemicals such as an increase in gastrointestinal absorption of lead, inorganic mercury, and other metals in children [31].

Apart from children and the elderly, TK differences are observed in numerous disease populations. Perhaps the most notable of these is liver cirrhosis. Numerous liver diseases exist, with cirrhosis being the most common, and itself having numerous etiologies, which may differentially influence PK [42]. In general, cirrhosis results in impaired hepatocyte function, changes in blood flow, and a reduction in plasma proteins, which subsequently result in changes in drug clearance due to altered liver size, CYP expression, plasma protein binding, and hepatic blood flow [43, 44]. Changes in glomerular filtration rate are also frequently seen [43]. Numerous studies have shown a reduction in enzyme and transporter expression in cirrhosis that correlates with disease severity; these generally result in a reduction in clearance and increase in the area under the curve (AUC) [43, 44]. Non-alcoholic fatty liver disease is the primary contributor to cirrhosis, even becoming more common in pediatric populations, and has its own set of symptoms, including changes in gastrointestinal pH [42]. Many symptoms of non-alcoholic fatty liver disease can be linked to comorbidities, such as obesity, which is the likely source of the reduction in cardiac output and CYP2E1 activity and expression seen in non-alcoholic fatty liver disease [42].

Indeed, obesity, which now affects 42% of adults in the U.S., results in a number of physiological and biological changes that influence TK [45,46,47]. The elevated body mass index of obesity increases cardiac output and liver and kidney weight and blood flow, while adipose blood flow is decreased [45, 47]. These physiological changes result in an increase in the V_d of lipophilic compounds and subsequently their half-life [46, 47]. Apart from increases in clearance by CYP2E1, there is also a reduction in CYP3A4 expression [42, 47]. Changes in plasma protein levels are also seen [47].

Another disease that results in significant alterations in PK, including in clearance, transporter function, and V_d, is chronic kidney disease [48, 49]. Kidney disease results both in a reduction in glomerular filtration rate as well as in active tubular secretions [50]. Plasma protein binding is also altered, causing an increase in uremic solutes and their subsequent accumulation in blood and tissues, which has been associated with a reduction in hepatic metabolism and both hepatic and renal transport [50]. PK is affected not only for renally eliminated drugs, but also for non-renally cleared drugs, particularly those cleared through the CYP2D6 pathway [49]. Each of these factors tend to increase with disease severity [48].

Recent work shows that one critical driver behind the changes in PK that result from disease, particularly chronic inflammatory diseases, seems to be inflammation [44]. For instance, the reduction in enzyme and transporter expression seen in liver cirrhosis has been linked to the upregulation of inflammatory cytokines [44]. An elevation in inflammatory markers and adipokines is also seen in obesity [45]. For additional perspective on the impact of inflammation on TK, we refer the reader to a recent review [51].

Genetic factors

Genetics can play a major role in interindividual variability throughout the lifespan, depending on the functional consequences of the genetic variant and the compound [52]. Polymorphisms have been shown for nearly all the CYPs, resulting in functional variability, particularly for CYP2D6 and CYP2C19 activities, leading to differences in metabolism and subsequently in chemical half-lives [53]. CYP2D6 is commonly cited as an example of the role of genetic variability on interindividual variation due to its major role in metabolism of pharmaceutical compounds and wide distribution in patient responses. There are over 133 CYP2D6 variants, leading to 60–100-fold variation in CYP2D6 activity [10]. Varying phenotypes resulting from common alleles can be grouped as poor metabolizers—those that show an increased and delayed AUC, or extensive metabolizers and ultrarapid metabolizers—those that show a lower, shorter, AUC (Fig. 1) [54]. Poor metabolizers make up 1–8% of the population, intermediate metabolizers 0.4–11%, normal 67–90%, and ultrarapid 1–21%, depending on ethnicity [55, 56]. Other CYPs that have shown variability with clinical relevance include CYP2C9, the primary enzyme in warfarin metabolism, and CYP2C19, the major enzyme in omeprazole metabolism. Dosing for warfarin is commonly determined based on CYP2C9 genotype of the patients [57]. Phase II enzymes can also exhibit phenotypic variation. Of the Phase II uridine 5'-diphospho-glucuronosyltransferases (UGTs), UGTs 1A1, 1A7, 2B15, and 1B17 are the enzymes most likely to show functional consequences [52]. Sulfotransferase 1A1 is also suggested to have phenotypic variations. Further details on six such polymorphic xenobiotic metabolizing enzymes (i.e., CYP2D6, CYP2E1, aldehyde dehydrogenase 2 (ALDH2), paraoxonase 1 (PON1), glutathione transferases (GSTs), and N-acetyltransferases (NATs)) can be found elsewhere [12]. Polymorphisms in transporters may also contribute to TK variability, one example of which has been reported for organic cation transporter 2 (OCT2) [58].

The frequency of certain polymorphisms varies by ethnicity, which contributes to demographic differences in metabolism. Allelic frequencies for many metabolizing enzymes have been estimated, particularly for the groups most commonly studied, such as Caucasians and Japanese [59, 60]. Variation is not solely due to polymorphisms in the coding region of the enzyme. The phenotypic variability may also arise from variants in noncoding regions such as the promoter, microRNAs that regulate expression, and other transcriptional regulators [10]. While many polymorphisms have been identified, not all of these have been characterized, limiting the ability to draw conclusions on functional consequences of these polymorphisms. It should also be noted that polymorphisms may lead to heightened or poor metabolism, which may or may not lead to an adverse outcome; so although TK variability may be noted, whether it results in TD differences needs to be determined [12].

Ontogenetic factors

Ontogenetics, also known as developmental pharmacology, refers to the study of how compounds are differentially processed and metabolized in the body throughout the course of development. It plays a major role in exposure response variability in children, which particularly involves enzyme and transporter ontogenies [40, 61]. Most metabolic enzymes are present at very low levels, if at all, during fetal development [31, 62]. CYP3A7 is an exception, with its highest expression during the first trimester, staying stable throughout gestation, until its expression levels fall after birth, when CYP3A4 begins to predominate among the CYP3As. Some enzymes, such as CYPs 2C9, 2D6, 2E1, and 3A4, and UGTs show a rapid increase in expression after birth, though this is highly variable [62, 63]. Other enzymes are slower to express, such as CYP1A2, which does not reach adult levels until approximately one year of age [31, 64]. Despite a lack of certain enzymes, there may be metabolic compensation for some compounds due to substrate overlap. For instance, CYP3A7 shares many substrates with CYP3A4, and sulfotransferases share several substrates with UGTs, providing alternate routes for these to be conjugated and eliminated [52]. Intra- and inter-study comparisons show that there are windows of biological hypervariability in enzyme development, adding yet another layer of complexity to neonatal metabolic function [65]. Although life-stage-specific abundances of CYPs and UGTs readily involved in drug metabolism are well-characterized, several enzymes involved in non-drug chemical metabolism have only limited data available, leading to greater uncertainties in such evaluations. As the field matures more data are being collected [66,67,68]; indeed, more experimental data are needed to provide greater confidence in their application, in modeling or in development of quantitative structure–property relationship prediction tools.

Evaluation of transporter abundance and developmental changes will also be important in characterizing TK variability. As membrane-bound proteins present in numerous tissues (e.g., kidney, liver, brain, placenta), transporters actively uptake or efflux endogenous and exogenous compounds, thus impacting chemical distribution in the body [69, 70]. With the emergence of targeted proteomics as a robust experimental tool for membrane protein abundance, tissue-specific, ontogenetic evaluations have been a highly active research area in recent years [16, 71]. It is worth noting that levels of P-glycoprotein (P-gp), an efflux transporter present in several metabolic tissues, such as brain, placenta and mammary tissue, have consistently shown lower levels in preterm and newborn samples compared to older children and adults [72]. An investment in understanding age-dependent abundances of placental and mammary transporters will yield critical information that could be directly applied in maternal–fetal PBPK models to evaluate fetal drug or chemical exposure [73]. Recent efforts that compile documented interactions between environmental pollutants and transporters have noted a range of activities, spanning from substrates for efflux transporters (endosulfan, methoxychlor: P-gp); uptake transporters (2,4-D, Aflatoxin B₁: organic anion transporter 3). PBTK models will be critical in broadening our understanding of transporter involvement in TK variability.

Exposomic factors

Exposomics aims to comprehensively assess all types of exposures and/or stressors that an individual may encounter over their lifetime [74]. Exposomic factors cover a huge range, from lifestyle-related to environmental and geographic factors, and have the potential to modulate both TK and TD variability. Alcohol and tobacco consumption, diet-related factors (e.g., nutritional status, cruciferous vegetable intake), obesity, and contraceptive use, all impact expression of a range of Phase I and II metabolic enzymes [52]. Xenobiotic co-exposure through drug use, particularly polypharmacy, can often impact metabolism, with drug-drug interaction a major research area for the pharmaceutical industry [75]. This is particularly relevant for the elderly, who are often taking many medications at once. Apart from drugs, other xenobiotics may impact metabolic rates, such as cadmium and lead exposure, or exposure to environmental carcinogens more generally [76]. It is of particular note that xenobiotic exposures may have critical windows of susceptibility, most of which occur during embryonic development or the first few years of life [77]. In utero chemical exposures have been shown to cause gene expression changes in the fetus, which could lead to changes in enzyme expression. For children, it is important to also consider that they have more hand-to-mouth activities and floor contact.

Consideration of the influence of gastrointestinal microbiota on metabolism is an emerging exposomic factor of relevance for both TK and TD variability [19]. Although unlikely to be important for xenobiotics rapidly absorbed in the upper small intestine, it may be a factor for low-solubility, low-permeability compounds and those subject to biliary excretion [18, 78]. Interindividual host microbial genome differences have been attributed to differences in gut microbiome composition and drug metabolism [79, 80] and as such are likely to elicit variability in microbial-compound metabolism. Moreover, many environmental pollutants have been shown to be subject to the reductive and hydrolytic microbial enzymes such as azoreductases, nitroreductases, and β-glucuronidases, including nitrotoluenes, pesticides, polychlorinated biphenyls, and azo dyes [81]. Additionally, many of these chemicals can affect or inhibit microbiome composition as demonstrated in animal studies, including chlorpyrifos, bisphenol A, and traffic-related air pollution [82,83,84]. The interplay between these interactions and their implications for human health assessment are quite complex, still requiring significant elaboration across multiple experimental models. More analysis and insights can be found in two excellent reviews on this topic [18, 19].

Until recently, an often overlooked exposomic factor is societal, which encompasses a myriad of factors including dietary, social, and environmental stressors. Consideration of geographic factors that can impact xenobiotic exposure, such as the toxic substance profile of the area, housing quality (e.g., asbestos and paint), the types of heating and cooking, pest control methods, urbanization (e.g. traffic-related air pollution), and climate (e.g., exposure to diseases and the types of outdoor play areas) is important [85]. Moreover, these factors are often related to environmental justice, with certain populations experiencing a greater burden than others, e.g., non-white people living in areas with increased pollution and toxicant exposure propagated from financial service discrimination (i.e., red-lining), with simultaneously elevated stress due to housing and economic pressures [85,86,87]. Indeed, a consensus recommendation following a 2020 National Academies of Science Engineering and Medicine workshop focused on strategies to integrate aging and environmental health research called for the implementation of a compound exposome approach to ensure adequate consideration of unique social and neighborhood factors [88].

Additional considerations

Apart from the discussed contributors to variability, it is likely apparent that delineations between such factors are not always clear-cut. Some population characteristics, such as sex, disease, and ethnicity, contribute to TK variability in overlapping yet distinct ways; consideration of population inputs regarding physiology and genetics for instance may require a nuanced approach. Another potential confounder in this field is the bias of available data, as the majority of studies are based on pharmaceutical compounds and healthy Caucasian individuals. Relevance to environmental pollutants possessing disparate physicochemical properties and other subpopulations is unclear without further evaluation. On the other hand, not all factors that may contribute to variability will have an impact on TK or TD leading to an adverse effect. For instance, if blood flow is the limiting factor to metabolism, variability in metabolic rate or enzyme levels may not impact steady-state levels [12]. Variability in metabolism and whether it is ultimately linked to an adverse outcome depends on whether such metabolism leads to clearance from the system, or formation of a reactive, toxic metabolite [9]. As such, problem formulation and case-by-case considerations are important to ensure TK and TD variability are appropriately considered and evaluated for decision-making.

PBTK modeling and incorporation of variability

PBTK modeling

PBTK methods can range from noncompartmental—where xenobiotic blood concentrations per unit time measurements are used to inform estimates—to compartmental methods that can represent varying degrees of complexity. The concept of multi-compartmental PK models for the simulation of PK data was introduced as early as 1937 (Teorell), but widespread use has only occurred recently largely due to initial limitations of mathematical complexity and presumed data requirements [89]. The advent of reliable tools to predict required inputs (e.g., V_d) using physicochemical properties has facilitated increased usage [89,90,91]. Initial problem formulation is required to ensure all relevant scenarios for all relevant ADME characteristics are considered: for instance, which exposure routes are relevant (e.g., dermal, inhalation, ingestion) will in turn dictate which routes of absorption need to be addressed and how bioavailability (i.e., amount of chemical reaching the bloodstream) is to be calculated [2].

Compartmental models can range from more simple forms that consider little more than a central plasma compartment connected to peripheral compartments, with the application of rate constants to estimate PK; to more complex models that incorporate multiple tissue types and accompanying biology and physiology (Fig. 2) [92]. In more complex PBTK models designed to evaluate population variability, input parameters are described as either system-level (i.e., population-specific) or chemical-specific, wherein either TK inputs for a particular chemical and/or the accompanying physicochemical properties are considered to derive the needed input parameters (Fig. 3). If in vitro or in silico estimates are employed, in vitro to in vivo scaling may be required. These approaches of using in vitro or in silico measurement to estimate parameter values for populating PBTK models are referred as a bottom-up approach [93]. To incorporate variability in such approaches, information on the variability (e.g., probability distributions of a parameter for the population being modeled) can be incorporated during model simulation, as captured in Fig. 2.

Fig. 2

Evolution of PK Models. A. General structure of Classical TK, Generic PBTK, and Complex PBTK models. B. Monte Carlo (MC) simulations can be used to incorporate interindividual variability. k_a: rate of absorption; k_e: rate of elimination; CL: clearance; Q: blood flow

Fig. 3

PBTK model input streams and considerations

PBTK modeling and variability simulation

Incorporation of variability in PBTK models can be performed to capture both physiological and biochemical variability to simulate TK in a population of individuals rather than an average subject. Equations describing distributions of system parameters for a model are derived from distributions of data based on real populations or patients. PBTK models can then be used to estimate internal concentrations from external exposures for these virtual individuals, which can arise from any population or life-stage of interest for which such data are available. The proprietary software Simcyp, a modeling and simulation software supported by a consortium of pharmaceutical companies with a mutual need for modeling and simulation [55], has a varied range of preparameterized population libraries that encompass human gestational, pediatric, and elderly life-stages and several ethnic and disease populations, as well as a few preclinical species [43, 94, 95]. Of the open-source tools designed specifically for chemical PBTK modeling, PopGen possesses a virtual human population generator that can be used to predict anatomical, physiological, and Phase I metabolic variation in healthy humans using human biomonitoring data [96]. Another platform designed specifically for NAM chemical risk evaluations developed at the United States Environmental Protection Agency (EPA) is the httk R package, which can make predictions of TK profiles for multiple chemicals in various species [97]. The httk package includes different types of TK models with varied complexity that allow the incorporation of in vitro TK data, physicochemical properties, and species-specific physiological data. With each type of TK model, the plasma or tissue concentration–time profile can be generated using forward dosimetry for individuals or populations. In the simplest model of the httk package, only intrinsic clearance rates (Cl_int) and fraction unbound in plasma (f_up) values are needed to run the model.

Virtual population simulations to capture physiologic variability

Both Monte Carlo (MC) simulation and Bayesian approaches can be used to incorporate population variability in PBTK modeling by generating and simulating virtual populations. MC, an approach that can be correlated to ensure variability is constrained to biologically feasible limits, is the preferred method for incorporating population variability [94, 96,97,98]. One example of MC application is in the httk package HTTK-POP, which incorporates interindividual variability for ten U.S. demographic groups based on National Health and Nutrition Examination Survey (NHANES) biomonitoring data [99]. To develop the representative populations, MC sampling of physiological parameters for the ten demographic groups from the NHANES cohort was performed, and then used to simulate interindividual physiological variability for each group, as well as experimental variability in f_up and Cl_int. As isozyme-specific information was not available, five percent of the population was assumed to be poor metabolizers. Exposure information was also considered for each of the demographic groups, showing the elderly population to be at the greatest risk, due to reduced clearance.

With incorporation of measures of human TK, such as Cl_int and f_up, the distributions of doses that could elicit bioactivity can be predicted, providing a measure of interindividual variability. Interindividual TK and physiological variability have been incorporated in an IVIVE-PBTK approach using Simcyp [29]. The Simcyp PBTK model offers advantages over other PBTK models for assessing interindividual variability in TK and TD, as it provides a range of subpopulation and life-stage libraries that include information on variation in isozyme expression as well as physiological characteristics [29, 100]. Information on isozyme expression allows for consideration of both genetic and ontogenetic contributors to TK and TD variability, providing more accurate estimates of Cl_int, particularly for children under 5 years of age [94, 100, 101]. Wetmore et al. generated isozyme-specific clearance rates on nine chemicals across 13 different Phase I and II enzymes and incorporated this data into Simcyp, performing PBTK modeling for a range of sensitive groups to identify the degree of variability in steady-state concentration (C_ss) concentrations [29]. This bottom-up approach utilizing Cl_int values at the level of individual isozymes allows great flexibility in revising system parameters, such as for new demographic, genetic and physiological data to make assessments of interindividual variability in TK parameters. Several other studies have incorporated TK variability using similar approaches [102, 103]. A study on methyleugenol used both individual human liver fractions and specific isozymes to estimate clearance rates and combined this in vitro data with literature-derived interindividual variability in a PBTK model with a MC simulation approach [102]. A study by the same group incorporated a developmental endpoint for calculating a phenol chemical-specific adjustment factor (CSAF), making for a semi PBTK-PD approach. Interindividual variability in oral absorption and UGT expression were incorporated into the PBTK model using MC simulation [103].

PBTK modeling can also help to elucidate the key drivers underlying the manifestations of population variability in PK seen in disease, as well as optimize dosing in disease populations [43, 44, 46, 49]. Indeed, the European Medicines Agency supports the use of PBTK models for clinical study design and a number of commercial and open-source PBTK modeling platforms include disease populations in their platforms [43, 44, 104, 105]. These populations can be continuously refined based on additional data.

Moreover, sensitivity analysis can be performed to identify the major drivers of interindividual variability. Al-Subeihi et al. [102] found bioactivation by CYP1A2, epoxidation by CYP2B6, and the apparent kinetic constants for oxidation and sulfation to be the greatest contributors to metabolic variation between individuals. Grzegorzewski et al. [106] explored the contribution of physiological parameters on CYP2D6 metabolic activity, showing these to have little impact and contribute similarly, independent of CYP2D6 activity. In a PBTK model of bromodichloromethane in pediatric populations, CYP2E1 was found to be the greatest contributor to variability in pharmacokinetics [107]. Such sensitivity analyses are useful in focusing efforts to better understand population variability.

IVIVE

IVIVE, a process using in vitro data to predict in vivo phenomena, is a critical translational step in the conversion of NAMs into metrics that are useful in risk assessment. IVIVE is used both to scale in vitro measures of TK parameters to the corresponding in vivo parameter, or to convert an in vitro concentration, typically related to assay bioactivity, to an in vivo effect dose, which can be compared against in vivo data or exposure [28].

IVIVE to support TK modeling requires assay- and parameter-specific considerations are met to ensure quality data are generated and appropriate scaling is achieved. Many of the assays designed to evaluate key parameters such as plasma protein binding and hepatocyte clearance have undergone rigorous evaluation as they were developed for use in the pharmaceutical industry, with more recent evaluations conducted to support guidance provided by the Organisation for Economic Cooperation and Development (OECD) [108,109,110]. Good in vitro method practices ensure that approved methodologies are followed, in conjunction with the associated negative and positive controls and any reference compounds to ensure assay reproducibility is maintained. Additional consideration is given to the use of extrapolation factors that may be required in instances where variances across in vitro systems (e.g., hepatocytes vs. microsomes) or assay conditions (differing but acceptable pH conditions) require certain adjustments for normalization. Finally, scaling factors required to scale up to represent organ-level clearance may vary, depending on the population or species of interest.

Incorporating TK and high-throughput screening data for equivalent dose estimation

IVIVE can also translate a bioactivity concentration from an in vitro assay (which can include TK/TD variability) into human-relevant doses. In this type of IVIVE application, reverse dosimetry is performed by incorporating chemical TK, which is typically realized by using a PBTK model, and to calculate an in vivo equivalent dose that could lead to a plasma or tissue concentration equal to the in vitro bioactivity concentration. The typical method to calculate the in vivo equivalent dose is described below:

$${\text{Administered}}\; {\text{equivalent}}\; {\text{dose}}\; \left( {\frac{{\frac{{{\text{mg}}}}{{{\text{kg}}}}}}{{\text{d}}}} \right) = \left( {\frac{{{\text{1mg/kg/d}}}}{{{\text{C}_{\text{ss}}}\left( {\mu {\text{M}}} \right)}}} \right)*{\text{in vitro}} \;{\text{bioactive}} \;{\text{concentration}}\; \left( {\mu {\text{M}}} \right)$$

Alternate methods to calculate points of departure for risk assessment are available as well, such as based on dose–response curves [111]. A key assumption to operationalize this IVIVE approach is that human plasma concentrations are equivalent to in vitro assay media concentrations. While nominal assay chemical concentrations are typically assumed to be the concentration available to elicit an effect, differential partitioning of chemicals between plastic, cells and media may affect the accuracy of such nominal concentration-based estimates, leading to a potential over/underestimation of bioactivity [112, 113]. Current practice uses nominal concentrations due to limited empirical data to evaluate available models; however, ongoing efforts to refine existing models, curate assay-specific details, and generate evaluation data may yield a path forward to incorporate distribution characteristics in future potency estimations [114,115,116].

The general PBTK-IVIVE workflow is outlined in Fig. 4. IVIVE for TK variability is perhaps best demonstrated in Wetmore et al. in which measured chemical TK parameters—hepatic clearance and f_up—were incorporated into an IVIVE-PBTK approach [117]. In addition, in vitro bioactivity data were incorporated to derive the oral equivalent dose, or administered equivalent dose (AED), that would produce a C_ss that might elicit bioactivity.

Fig. 4

General PBTK-IVIVE workflow. Various exposure routes can be simulated. In vitro TK data is generated and used to inform PBTK model through IVIVE of ADME parameters. In vivo TK data (e.g., plasma protein binding (PPB), transport, hepatic clearance) can also be used to inform the PBTK model. Estimated internal concentrations (e.g. C_plasma) are combined with in vitro toxicodynamic points of departure to derive an administered equivalent dose (AED), representing an external exposure that could elicit a similar bioactivity in vivo. These AEDs are compared against human exposure estimates to derive bioactivity exposure ratios (BERs), an ad hoc margin of exposure

Exposure risk evaluation using bioactivity data: bioactivity exposure ratios (BERs)

Using IVIVE, in vitro effective concentrations can be used to generate AEDs which can then be compared against exposure estimates to derive a BER, a margin of exposure metric, to inform chemical prioritizations for risk assessment [118]. The lower the ratio, the higher the potential for the exposure to elicit bioactivity.

In vitro bioactivity data for IVIVE can be obtained from large scale screening studies performed across numerous chemicals and covering a range of endpoints. Two such major data resources from the U.S. are the Tox21 and ToxCast programs. Tox21/ToxCast is a collaborative effort between the National Institutes of Health (NIH), FDA, and EPA to develop and implement in vitro high-throughput assays to screen thousands of chemicals across hundreds of bioactivity or toxicity endpoints [11, 118]. Chemical potency from these assays is commonly summarized by an activity concentration exerting 50% of maximum response (AC50) or lowest effective concentration value [119], although other metrics can be used to define in vitro points of departure. It is worth noting that while the Tox21 and ToxCast assays provide a valuable resource for obtaining in vitro bioactivity data, due to their high-throughput nature, potency values derived from these assays should not be taken at face value without a thorough review of data quality underlying concentration-response curves. Bioactive concentrations should always be confirmed, and potential flags such as borderline activity, curve overvitting, and background noise should be carefully assessed. The subsequent chapter will highlight some of the NAMs employed in these IVIVE and risk assessment approaches.

NAM applications of IVIVE and modeling to inform chemical risk evaluations

NAMs incorporating metabolism and physiology to quantitate life-stage TK variability

Some measures of human TK variability are well-characterized as they are readily measured or obtainable—these are primarily physiological parameters, such as blood flows and tissue volumes [120]. In contrast, fewer data are available for the effect of genetic and ontogenetic factors, such as clearance rates, and far less is known regarding exposomic factors. NAMs here hold a particular advantage in that they can be performed using human-specific biology and in a high-throughput manner, allowing for an increased sample size, which is necessary when evaluating interindividual variability. Such NAMs include in vitro assays for both TK and TD variability, ranging from subcellular measures of metabolism and transport to cellular and tissue-level endpoints, and in silico prediction tools [29, 121].

In vitro assays that can measure human TK variability exist for most tissue types, including liver, lung, and skin [24, 122]. Wetmore et al. developed a screening approach to characterize TK variability across different populations and life-stages by measuring isozyme-specific clearance rates [29]. In this initial approach, a set of nine chemicals was screened across 13 isozymes. Another major contributor to TK variability is transport of chemicals. Several in vitro approaches are available to measure transporter activity. One common technique is to overexpress certain transporters in cell lines, such as HEK293 or Caco-2 cells, and then measure cellular uptake of compounds [58]. In addition, inhibitors can be used to measure the contribution of individual transporters [123]. With the incorporation of known variation in transporter expression between populations and life-stages, these measures can help to estimate variability in transport process and subsequent C_ss using PBTK modeling.

NAMs for exposure prediction

While emphasis One major resource for deriving human exposure estimates in the U.S. is the NHANES study from the Centers for Disease Control, which measures exposure biomarkers—primarily metabolites—from a diverse array of individuals across the U.S. in both blood and urine [124]. To expand use of human exposure data from NHANES and incorporate the National Research Council (NRC) request for predictive tools to be used for safety assessment, ExpoCast was developed by the EPA as a high-throughput exposure prediction model [125]. In addition to chemicals measured in NHANES, ExpoCast provides exposure predictions for thousands of additional chemicals, covering nine life-stages and demographic groups.

TD variability

While emphasis is often placed on TK variability, TD variability is another driver of population variability [52]. One of the first NAMs developed to assess human TD variability used a set of 146 lymphoblastoid cell lines that were characterized through the 1000 genomes project [11]. This population genomics approach has particularly helped to characterize genetic contributors to human population variability. It has also been used to characterize TD variability under complex scenarios, such as in response to chemical mixtures and identifying a chemical’s mode of action. TD variability can also be assessed in in vitro studies using donated human samples obtained through clinical trials. These NAMs, such as those from Genoskin and the Institute for In Vitro Sciences, use patient samples from a diverse pool of donors covering a range of ethnicities and ages, which can be used to assess interindividual variability for endpoints such as efficacy and immune responses [122, 126].

An exciting novel avenue to address interindividual variability has appeared with the development of methods to culture human induced pluripotent stem cells (hiPSCs). hiPSCs can be generated from somatic cells of individuals with different genetic backgrounds, which can be used to represent the genetic diversity of the population. Burnett et al. developed a model to measure TD population variability in cardiomyocytes derived from hiPSCs that shows strong reproducibility, sensitivity, and specificity [121]. Using hiPSCs from 43 individuals, a size similar to that used in human clinical trials, Burnett et al. screened 134 chemicals for functional and viability endpoints. They found interindividual variability to be the primary contributor to total variability, with little contribution from sex and ancestry. Using such an approach, they were able to characterize interindividual variability in cardiomyocyte responses and look at specific genetic contributors to human TD variability. Due to their human-relevance, usability in biochemical, genetic, and genomic approaches, and their capacity for assessing endpoints and mechanisms, NAMs using hiPSCs are being developed for a range of other tissues and endpoints [127, 128].

TD can be linked to PBTK models to create a PBTK-TD model [129]. This model allows for tying the concentration at the site of action to the downstream effect of the compound and can help to provide mechanistic information underlying this effect. Further details on such PBTK-TD models can be found elsewhere [130, 131].

Consideration of population variability in regulatory decision-making

Incorporating IVIVE and modeling approaches to inform population variability has a role in several arenas, including to support pediatric drug development and dose selection in the pharmaceutical industry [132], and informing setting testing priorities and evaluations of population variability in the toxicology arena [27, 118]. In the pharmaceutical industry, PBTK modeling and simulation (i.e., model-informed drug development) can help determine clinical study needs, study design including dosimetry, and the need for population-specific product label warnings [133]. In the toxicology arena, the release of two NRC reports—“Toxicity Testing in the 21st Century”, advocating for the use of what were to become NAMs to inform toxicity assessment [134], and “Science Decisions: Advancing Risk Assessment,” advocating for the need to include population variability in toxicity testing [135], spawned numerous efforts utilizing IVIVE and PBTK models to evaluate population variability [29, 117, 118, 136,137,138].

UFs and CSAFs

Population variability has traditionally been accounted for with the use of uncertainty factors (UFs) in risk assessment, with the default UF for human safety assessment using animal data as 100X, consisting of a 10X factor for interspecies variation and a 10X factor for interindividual variability [139]. The 10X for interspecies can further be broken down into 4X for TK and 2.5X for TD, and the 10X for interindividual variability into 3.2X for both TK and TD, allowing for adjustment of these different types of uncertainties when additional information is available [140, 141]. These factors and subfactors have been adopted by organizations worldwide, including the World Health Organization (WHO)/International Programme for Chemical Safety (IPCS), European Food Safety Authority, Health Canada and the EPA [20, 27, 142].

In those instances where quantitative data on TK and TD are available, most commonly in a higher tier risk assessment, there is an opportunity to depart from the default UFs in favor of CSAFs. Guidance advising on CSAF derivation was published in 2005 by WHO/IPCS, an outcome of an international harmonization project [142]. In essence, availability of data-derived factors that allow grouping due to 1) non-mode of action information (e.g. allometric scaling, clearance-related) or 2) mode of action-related that incorporates TK and/or TD data allows for biologically-based adjustment factors that reduce uncertainty estimations, allowing for departure from the default UF.

Some investigative (i.e., non-regulatory) evaluations using modeling have shown that a UF of 10X may be insufficient for certain populations and life-stages, including children and the elderly, and for chemicals metabolized by polymorphic enzymes [143, 144]. Strikwold et al. [103] provided an example of a CSAF approach for phenol in which a CSAF of only two was calculated based on the degree of interindividual variability in the maximal plasma concentrations. Wetmore et al. found the degree of variability, indicated by the human TK adjustment factor (HK_AF), to vary from 1.3 to 13.1, with the pediatric life-stage typically showing the highest internal dose (C_ss) [118]. In an effort that estimated TK and TD variability for complex mixtures, Abdo et al. found that the default UF aligned well with the TD variability of the two mixtures tested, although whether this result holds for other mixtures is unknown [145]. Uncertainty estimation relevant for the specific problem under study is also an important consideration that requires further elaboration [146].

A recent state of the science update [63] reviewed regulatory and investigative CSAF development and consideration by regulatory bodies to provide next steps to foster future adoption of CSAFs. Preferably, CSAF development based on PBTK models should, for increased confidence, capture (1) the biological basis of the model structure and parameters; (2) performance of the model comparison and model simulations with experimental data; and (3) reliability of model predictions of dose metrics relevant to risk assessment (i.e., model testing, uncertainty, and sensitivity analyses).

Children

Children provide perhaps the most frequently assessed case of TK variability and its incorporation into human health risk evaluation. As a more consistent approach than simply deriving CSAFs “where such data exists,” Ginsberg et al. outlined a framework to guide risk assessments for childhood chemical exposures. In the paper outlining this approach, the ontogeny of metabolizing systems with chemical-specific exposures is used to compare child versus adult internal dose for five case study chemicals [9]. This approach is outlined in Fig. 5.

Fig. 5

Approach for incorporating TK variability for children into risk assessment, adapted from Ginsberg et al. [7]

Tiered testing strategies have been proposed to identify chemicals for further prioritization and making decisions on which populations to target for risk assessment [23]. The approach outlined by Ginsberg et al. could be incorporated into a tiered testing strategy, and a similar strategy might be developed for other subpopulations and chemical groups.

Drug development considerations

Approaches used in the pharmaceutical realm can provide additional examples for assessment and integration of TK variability that could be harnessed. In the drug safety arena, the U.S. FDA has been utilizing PK information and modeling to support drug safety for a number of years [11, 147]. The FDA accepts PBTK models to inform drug and biologics applications for potential waiving of clinical PK data [148]. Much of the current trend at the FDA is focused on predictive models, using quantitative structure activity relationships (QSARs) and PBTK models, such as for predicting drug-related toxicities, including drug-drug interactions, falling in line with the trend in the pharmaceutical industry. In 1991, roughly 40% of drug candidate attrition was due to ADME concerns. With increased use of ADME predictions and PBTK modeling, in 2000 that number dropped to 9% of all failures [149]. PK modeling can be used to predict therapeutic doses for certain populations, such as pediatrics, or other understudied patient populations. PBTK modeling on pediatric populations has received particular attention from stakeholders and is becoming commonplace in Phase I trials since companies are now required to submit data on pediatrics beforehand [92]. PBTK modeling has been used to determine dosing for children for numerous drugs, such as dolutegravir in neonates [150, 151]. In 2018 the FDA issued draft guidance for inclusion of pregnant women in clinical trials [152]. Similar PK modeling approaches have been used for evaluating variability brought on by diseases, such as renal impairment [42, 94, 153].

PBTK modeling has also allowed for estimating the contribution of multiple variables, enabling clinicians to adjust dosages for populations that might be more or less sensitive, contributing to model-informed precision dosing [154]. This approach combines PK and PD models with patient-specific data to individualize drug dosing. It aims to optimize the balance between efficacy and toxicity based on patient-specific criteria. However, widespread adoption of precision dosing has not occurred, remaining localized in certain academic centers. Adoption has been limited due to various challenges, including lack of needed background knowledge on the part of clinicians, a dearth of genotypic data, and cost/benefit analyses that limit interest on the part of hospitals and insurance.

State of the science: PBTK models for variability evaluation

Providing greater specificity for certain subpopulations or chemicals of interest, PBTK models have been developed for application to scenarios of particular regulatory interest, such as pregnancy, per- and polyfluoroalkyl substances (PFAS), and mixtures. Case studies for each of these are discussed below.

Time-varying parameters and models to evaluate pregnancy and gestation

Pregnancy is a particularly complex life-stage for PBTK modeling, as many parameters that can typically be described as constant are time-varying [155, 156]. Moreover, components that are not present in adults need to be accounted for, particularly for the fetus, which has its own distinct compartments and physiology, such as blood flow through the ductus venosus and foramen ovale [156]. These physiological and anatomical changes may subsequently impact metabolism and clearance both in the mother and in the fetus [157]. Examination of TK parameters during pregnancy is particularly important as fetal development presents a critical window of susceptibility that could result in life-long effects [156]. Despite the importance of and significant interest in this life-stage, there is a major dearth of data as most drugs are “off-label” for pregnancy, as pregnant women are typically excluded from clinical trials [156, 158]. Because of this, there is a critical need for development of pregnancy-specific PBTK models. Both FDA and the European Medicines Agency have issued pregnancy exposure guidelines–the European Medicines Agency for monitoring exposures throughout pregnancy and the FDA for prospective exposure registries while supporting PK/PD modeling [157, 158].

The development of pregnancy PBTK models relies on large datasets of physiological and anatomical parameters to account for the time-varying nature of this information and a wide degree of interindividual variability [158]. A number of data curation studies have been published identifying time-varying parameters and where gaps exist in the literature [155, 156, 158, 158, 159]. From these reviews, algorithms have been generated that describe anatomical and physiological changes with gestational age. Anatomical and physiological changes typically start in the first trimester, peak in the second trimester, and then stay relatively constant through birth [157, 158]. Significant physiological changes are seen such as in body weight, cardiac output, protein binding, renal clearance, respiration, levels of sex hormones, and intrauterine volumes [156,157,158,159]. Differences in TK properties, e.g., f_up and Cl_int, are seen as well [155]. However, data are still sparse for certain parameters such as hepatic blood flow, early fetal growth, and metabolic enzyme activity, making it difficult to determine how these might be altered during pregnancy. Clinical observations of certain drugs support the observed physiological changes, such as increases in plasma volume resulting in increased V_d and reductions in plasma protein levels leading to increased f_up [159].

Algorithms from such data curation studies have been used to generate pregnancy PBTK models. The PBTK models can be further refined by including separate fetal compartments and/or considering the lactation phase [133, 155, 158, 160]. Comparisons of model predictions with in vivo experimental data demonstrate that these models provide reasonable estimates of concentrations, despite a high degree of uncertainty [155, 158]. Zhang et al. performed a sensitivity analysis to quantify how fetoplacental metabolism, gestational age, and placental transport impact fetal drug exposure following maternal drug exposure [133]. The umbilical venous: maternal plasma ratio, a common measure of fetal drug exposure, was found to be a poor indicator of fetal exposure except for under relatively steady-state conditions. A time-embedding network model underscored the importance of enzyme ontogeny—changes in enzyme expression—to chemical kinetics and toxicity [161]. This time-embedding model incorporated time-varying expression profiles for 23 enzymes, metabolites, and glutathione reactivity based on data on Americans eight weeks gestation to 18 years old, to model population and life-stage variability. Unlike most models that include a fixed degree of population variability, this time-embedding model allowed for capturing changes in variability over time based on the data. Such time-embedding models will allow for greater understanding of population and life-stage variability.

Applying a pregnancy model for PFAS

Pregnancy PBTK models have been further applied to specific chemicals, such as PFAS [160]. PFAS are a class of chemicals that have been widely used for their chemical stability and properties, such as water- and oil-repellence and are thus found in a wide range of products ranging from cookware to food wrappers and furniture coatings [160, 162]. PFAS have garnered major public interest and concern in the past several years due to their widespread prevalence in humans—PFAS are detected in the blood of > 99% of individuals in the U.S.—and potential health effects: PFAS have been linked to developmental effects and alterations in lipids and immune function [160, 162]. Due to their chemical stability, PFAS typically have long half-lives in humans—over 5 years for some longer chain PFAS [153, 160, 162,163,164]. This contrasts with relatively short half-lives in rats, necessitating the importance of human-specific data [162]. The mounting concern over PFAS has led to the phasing out of several classes of PFAS, particularly perfluorooctanoic acid (PFOA) and perfluorooctane sulfonic acid (PFOS) since ~ 2000, correlating with a drop in serum levels since then [160, 162]. While PFAS are detected in nearly all individuals, there is a huge degree of interindividual variability, which can be largely attributed to geographic-specific exposures [162, 165, 166].

Variability in PFAS levels is particularly evident in pregnant and postpartum mothers due to several factors, such as parity, breastfeeding, age, and body mass index (BMI) [160]. Prior pregnancies and breastfeeding reduce the body burden of PFAS by increasing the V_d and preferential partitioning of PFAS into lipophilic tissues, including breast milk. To better understand fetal PFAS exposure, Brochot et al. developed a pregnancy and lactation PBTK model to estimate PFOA and PFOS in mothers and placental transfer to the fetus [160]. The model was built using PFAS exposure estimates derived from INMA, a prospective birth cohort study in Spain for which maternal blood spots were taken in either the first or second trimester and umbilical cord blood at delivery. The model was based on a generic maternal PBTK model to which a previously published fetal model was fit. Oral absorption was described through food and drinking water with elimination via breastfeeding and excretion, as well as blood loss during delivery. A dynamic exposure model was used to capture the drop in PFAS exposure beginning in ~ 2000. This model allowed for characterizing tissue distributions of PFOS and PFOA and showed differences in tissue distributions for the two compounds that fit with their distinct physicochemical properties. Moreover, the model showed reliability and broader applicability as it produced estimates that correlate well with levels reported in the U.S. population.

More general PBTK models for PFAS have been built as well, with many of these based on the U.S. population as this was the initial site of PFAS production [167]. Worley et al. developed and evaluated a PBTK model that incorporated renal transport mechanisms using IVIVE and MC to simulate serum PFOA concentrations following exposure to PFAS at levels reported in drinking water. Biomonitoring data were collected from individuals residing in two former major PFAS manufacturing sites, Ohio and Alabama, before and after filtration systems in public water systems were introduced, which brought PFAS water levels to below the limit of detection [162]. Although the model calibration revealed excellent agreement between the measured serum PFOA data and the predicted values, model evaluation with the biomonitoring samples noted an underestimation of the serum concentrations for a few of the communities, indicating that variability in PFAS exposure through non-drinking water sources may be involved. Additional MC simulations of inter-individual variability and variability in exposure (e.g., variations in drinking water source, distance from PFAS source) confirmed that non-drinking water sources may have a greater impact on exposures than previously thought.

Additional models exist as well. Chou et al. developed a multi-species model that was later refined as a Bayesian dose–response model to refine TK and TD uncertainties, as well as allow for performing population level probabilistic risk assessment [168]. Indeed, EFSA has used several such PFAS-specific PBTK models in their assessments of PFAS [169]. These compound-specific models provide more refined estimates than a generic PBTK model can, which can be important for particular chemicals of concern [170, 171].

Lifetime exposure modeling

Lifetime exposure modeling tends to estimate cumulative exposure to a chemical over the course of a person's lifetime. By integrating data from various sources and considering the different factors influencing exposure, lifetime exposure modeling can provide a comprehensive assessment of the potential health risks associated with long-term exposure to a chemical, as well as assess the effect of exposure during critical time windows in late-life diagnosis pathologic endpoints.

Verner et al. developed a lifetime PBTK model to estimate lifetime blood/tissue exposure levels of persistent organic pollutants during hypothesized time windows of susceptibility in breast cancer development to assess the association between persistent organic pollutants and breast cancer incidence [172]. In the PBTK model, the values of physiologic parameters (e.g., organ volume, composition, and blood flow) throughout a woman’s entire life were estimated based on data on pregnancies, height, weight, and age. The lifetime TK profile for various exposure scenarios and physiologic factors (i.e., breastfeeding, growth, pregnancy, lactation, and body weight changes) was assessed. The study revealed that lactation periods and body weight changes are the factors that had the greatest impact on lifetime TK profiles.

In another study, a PBTK model was used to simulate blood levels of polychlorinated biphenyls during specific pre- and postnatal periods, which helped evaluate the association of chemical exposure with impairment of infant behaviors (e.g., attention, activity). Specific windows of susceptibility to individual infant behaviors were identified, highlighting the importance of modeling TK profiles during these periods (e.g., within the first year of life) [173].

When considering TK variability due to lifetime exposure, in addition to those factors discussed above, factors such as occupation, lifestyle, geographic location, and specific windows of susceptibility need to be included. The consideration of these additional factors shall reduce uncertainty linked to past chemical exposure and help to identify the impact of exposure during sensitive time windows to disease.

Mixtures

Perhaps the greatest risk assessment challenge is adequate consideration of mixtures or chemical co-exposures and downstream effects. Indeed, the range of co-exposure scenarios is infinite and exposures can also occur via multiple routes. Often, study formulation is driven by understanding the impact of chemicals likely to interact with one another. PBTK models provide a useful framework to evaluate how TK can be impacted by chemical interactions. Examples are provided below to build our understanding of the impact of mixtures on TK and TD which may in turn stimulate thinking for future application.

Chemical-chemical interactions are impacted primarily by dose, administration route, and the number of compounds. Interactions can influence TK by impacting transport, binding properties, biotransformation, or formation of complexes that may impact physicochemical properties such as lipophilicity [171]. Mixture interactions can be modeled if all binary interactions are known, to give a network of interactions, though this becomes exponentially more complicated when more components are added to the system [171]. To simplify the process, it is suggested that only the interactions that most impact kinetics are modeled, and chemicals with similar physicochemical properties are lumped together. Alternatively, models that capture only specific pathways or for specific chemical groups can be modeled on their own.

Perhaps the most common example in pharmacology of how chemical interactions can influence human TK is that of CYP3A4 and grapefruit juice—an inhibitor of CYP3A4 that also inhibits the efflux transporter P-gp. Conversely, St. John’s wort is a commonly used inducer of CYP3A4 and P-gp. Being a major concern for pharmacology, Quignot et al. developed a meta-regression model based on an extensive literature search on the degree of interaction—or TK modulation—between CYP3A4 and P-gp substrates and grapefruit juice or St. John’s wort, as well as the degree of variability in metabolism [174]. The degree of interaction in kinetics (e.g., maximal concentration, AUC) was calculated based on the ratio of the binary mixture over the substrate alone. Substrate bioavailability and fraction metabolized were found to be the greatest contributors to the degree of interactions. To characterize variability, UFs were calculated for bioavailability and fraction metabolized, as well as the degree of interaction. For single compounds, UFs were all below the default. For multiple compounds, however, UFs reached 18.9 for acute exposure, and 17.1 for chronic exposure, when capturing the lowest bioavailability to highest fraction metabolized, in the presence of grapefruit juice.

A common example of mixtures of environmental compounds is pesticides, particularly pyrethroids. Quindroit et al. developed a PBTK model based on reverse dosimetry to estimate exposures to four pyrethroids—deltamethrin, permethrin, cypermethrin, and cyfluthrin, based on levels of metabolites measured in biomonitoring studies as these compounds produce several of the same metabolites and are metabolized by common pathways [175]. Effects of the individual pyrethroids were combined by dose addition, as this has been well-characterized in the literature, and interindividual variability in metabolism was assessed. Even with the incorporation of the default UF for interindividual variability, exposures were found to fall below a threshold of concern covering the 95th percentile.

Several PBTK models have been developed to model exposures to mixtures of volatile compounds and estimate human TK interindividual variability. Volatile compounds are not as well-characterized and provide additional complications but have provided insight into how exposure to mixtures can impact human TK variability at different life-stages, and the types of interactions that may occur [176, 177]. Co-exposures to drinking water contaminants, including benzene, trichloroethylene, and toluene, impacted the variability index for high exposures, with only minimal impact from low exposures. Variability generally fell below the default UF apart from high multi-route exposures at early life-stages. Sensitivity analyses, in which uncertainty and the influence of different model parameters are estimated, can be particularly informative for these complex models to determine drivers of variability. A sensitivity analysis for this model showed the blood:air partition coefficient to be the main driver of variability in the interaction model for the CYP2E1 pathway. For inhalation exposures to a similar set of compounds, TK variability was impacted based on the size of exposure, subpopulation, and substance type [177]. The number of compounds in the model did not matter at the low exposure level (20 ppm). At the high exposure level (50 ppm), the only interaction effect seen was in pregnant women where an increase in the number of compounds resulted in an increase in the maximal concentration of benzene. More sophisticated models must be developed to more comprehensively cover the range of chemical mixtures humans are exposed to, but these models provide a starting point.

Conclusion: barriers to adoption and opportunities to advance the field

Incorporation of NAMs with PBTK modeling has allowed for making great strides in estimating variability for a range of sophisticated exposure scenarios, from evaluations of specific life-stages or diseased populations to cumulative exposures, lifetime exposures, or a combination thereof [11, 29, 120, 121, 147]. Data gaps remain, particularly for establishment of more comprehensive databases that capture ontogeny data for any and all enzymes and transporters to support chemical risk assessment. Significant progress was made for enzyme ontogeny data gathered over 15 years ago, with databases containing experimental ontogeny data [9, 40, 65, 178] and in vivo clinical PK data to support model evaluation [144]. One such database comprised information for 6 polymorphic enzymes—CYP2D6, CYP2E1, ALDH2, PON1, GSTs, and NATs, chosen based on knowledge of toxicity, genotype-phenotype, and linkage to environmentally-mediated disease to allow evaluation of how polymorphism can affect metabolism [12]. Although databases exist that capture physiologic parameters and population variability across typically healthy Caucasian populations, such data are largely absent for ethnic and diseased populations and specific life-stages requiring protection. Efforts to establish open-source databases of physiological parameters haven’t proven sustainable long-term thus far [95, 179]. Fortunately, a recent effort made great strides in capturing time-varying anatomical and physiological features for use in a maternal–fetal PBTK model [155, 156].

Including information on substrate turnover, how polymorphisms influence turnover, mode of action, and population variability in enzyme function allows for estimation of CSAFs, with certain enzymes and populations flagged that may warrant further consideration given values exceeding the default UFs. As more transporter and enzyme ontogeny data have been captured since then, an update and re-evaluation of existing data is warranted.

Despite progress to date, regulatory adoption of CSAFs for decision-making is still quite limited, although barriers to implementation have been identified [63]. A fundamental need is a call to modelers to provide transparent, explicit documentation of their model framework, parameterization and code to ensure model evaluation is possible. With the release of recent OECD guidance regarding PBTK model development along with several reports advocating for best practices in PBTK model reporting and ongoing stakeholder support for outreach and communication, these efforts should work to align subject matter experts with regulators grappling with interpreting modeling findings for regulatory application.

PBTK modeling has been demonstrated to be a powerful tool to evaluate interindividual and population TK and TD variability. With the advent of NAMs, many efforts further demonstrate its potential when applied with in vitro data streams, particularly for its ability to translate in vitro data out to in vivo relevant equivalent dosages. These applications have garnered much support from regulatory stakeholders, at least for certain decision-making needs [27, 138, 180]. As more sophisticated NAMs are developed, for instance microphysiological models to inform tissue-level responses, TK NAMs and models are poised to inform target tissue dosimetry. Although some gaps and barriers do remain that limit the ability to capture variability across all populations and life-stages that require consideration in risk assessment, research activities addressing these limitations are already underway, paving the way for future progress.