Toxics, Vol. 10, Pages 700: A Physiologically Based Pharmacokinetic (PBPK) Modeling Framework for Mixtures of Dioxin-like Compounds

1. IntroductionDioxin-like compounds (DLCs) are a group of structurally similar persistent organic pollutants (POPs), including compounds in the families of polychlorinated dibenzo-p-dioxins (PCDDs), polychlorinated dibenzofurans (PCDFs), and polychlorinated biphenyls (PCBs). In general, DLCs contain seven 2,3,7,8–substituted PCDDs, ten 2,3,7,8–substituted PCDFs, and twelve non-ortho- or mono-ortho-substituted PCBs, which share a common mechanism of toxicities [1,2]. PCDD/Fs are produced in a number of industrial and natural processes, including steel smelting, automobile exhaust, incineration, such as municipal waste and forest burning, paper production, and as byproduct during the synthesis of organic compounds, such as herbicides [3]. PCDD/Fs can be dispersed into the soil, water, and air, and then spread to the general population. PCBs are produced for usage in coolants, electrical insulators, capacitors, and carbonless copy papers [4]. These DLCs have been at the top of man-made emissions and identified in at least 126 National Priorities List (NPL) hazardous waste sites in the US [5].DLCs are highly toxic and can cause a variety of adverse health outcomes, including cancer, developmental defects, immunotoxicity, metabolic and endocrine disorders [1]. 2,3,7,8-Tetrachlorodibenzo-p-dioxin (TCDD) is the most toxic congener, and the remaining DLCs are generally less toxic. By activating the aryl hydrocarbon receptor (AHR), DLCs share a common mode of action. Humans are exposed to environmental pollutants existing in mixtures. While it is important to understand and estimate the health risk of individual DLCs, such as the reference compound TCDD or the most prevalent congeners, it is the mixtures in the environment that ultimately drive the AHR-mediated health outcomes. Environmental exposures to DLC mixtures have been linked to many health endpoints. For instance, analyzing the data from the 1999–2004 cycles of the National Health and Nutrition Examination Survey (NHANES), Lin et al. showed a significant positive association between all-cause mortality risk and the toxic equivalent (TEQ) values of serum DLC mixtures including 7 PCDDs, 10 PCDFs, and 9 PCBs [6]. Similarly, using the 1999–2002 cycles of NHANES, significant negative associations were found between serum total T4 and TEQ of 7 PCDDs, 5 PCDFs, and 4 PCBs, especially in older women [7]. In the longitudinal Russian Children’s Study, the serum TEQ of seven serum PCDDs was found to be associated with lower sperm concentration, sperm count, and total motile sperm count in young men [8]. Higher cancer-related mortality was suggested for residents near a Waelz plant where the air PCDD/Fs levels were high [9]. Therefore, developing the capability of making accurate predictions of the health risk of exposure to DLC mixtures has significant public health implications.The compositions of DLC mixtures in the environment are highly variable. For instance, those at the Superfund sites vary greatly in the types of constituent congeners and concentrations in the soil sediments within and between polluted sites [10,11]. The compositions of DLC mixtures that humans are exposed to, mainly through dietary intake by consuming fish and fatty meats, are also highly variable [12,13,14,15,16]. While concentrations in the sources or doses in the exposures provide some preliminary information on the toxicity potentials of DLC mixtures, they are still far removed from the apical health endpoint. More proximal to the biological effects, human blood or tissue concentrations, available through biomonitoring [9,17,18], autopsy [19] or in silico predictions, are the preferred metrics used for estimating the health risk of mixtures.Physiologically based pharmacokinetic (PBPK) models can quantitatively predict the relationship between the exposures and tissue concentrations by mathematically describing the physiological and biochemical processes involved in the absorption, distribution, metabolism, and elimination (ADME) of chemicals. Over the past three decades, there has been a rich literature on dioxin PBPK modeling for species including rodents and humans on various physiological stages including developmental, lactational, lifetime, and obesity [15,20,21,22,23,24,25,26,27,28,29]. The vast majority of these studies were focused on TCDD, and the four-compartment model developed by Emond et al. has been adopted by US EPA for dioxin risk assessment [30]. Nevertheless, these models are for single-congener exposures only, even in studies where a number of congeners were covered [15,29]. In general, adapting an existing PBPK model developed for one specific compound to a structurally similar compound in the same family is straightforward, since the same model structure can be used and for the most part it involves estimating chemical-specific parameters for the new compound. If at the exposure levels of concern there are no congener–congener interactions, e.g., the congeners do not compete for binding to kinetically important proteins or enzymes or they do not cross-induce common metabolic enzymes, then modeling the mixtures can be carried out simply by simulating each congener independently.A dioxin mixture PBPK model does not fall into this case. The pharmacokinetics of DLCs is nonlinear and dose-dependent, primarily due to processes in the liver that are specific to dioxins [31,32]. First, a dioxin compound can bind to the CYP1A2 protein in the liver such that its free concentration is greatly reduced due to the sequestration by the protein. Second, through binding to and thus activating AHR, a dioxin compound can induce the gene transcription of CYP1A2, which further enhances the sequestration of the compound and also the clearance of dioxin where CYP1A2 acts as the main metabolic enzyme [33]. As a result, for a DLC mixture exposure, different congeners will compete for binding to CYP1A2 and also collectively induce it to sequester and metabolize one another. Therefore, significant nonlinear congener–congener interactions could exist in the toxicokinetic process, and simulating each congener individually and subsequently combining the results linearly may not make accurate predictions of the mixture effect. In this situation, all congeners in a mixture have to be simulated simultaneously in an integrated mixture PBPK model.Traditionally, developing a mixture model would require manually extending the existing PBPK model constructed for a single compound to include new ordinary differentiation equations (ODEs) that would track the concentrations of new compounds added in the mixture and adjusting the existing ODEs for common variables, such as AHR and CYP1A2 in the case of DLC mixture, accordingly. This type of manual programming is tedious and error-prone and is not flexible as more compounds are considered. The lack of flexibility to efficiently scale up for multi-compound exposures constitutes a major technical barrier toward large-scale mixture simulations. The main goal of the present study is to develop a coding scheme or framework for DLC mixture PBPK models that can flexibly accommodate an arbitrary number of congeners without requiring manual modification of the model. The mixture model was constructed by adapting a TCDD model previously developed by Emond et al. [23,28].To estimate the toxicity potencies of DLCs relative to TCDD, toxic equivalency factors (TEFs) have been established for the congeners based on expert evaluations of a multitude of in vivo and in vitro experimental studies [34]. With the TEF of TCDD at 1, the TEFs currently adopted for other congeners range from 0.00003 to 1 with half order-of-magnitude increments [2,35]. To quantify the toxicity potentials of a DLC mixture in exposure, TEQ values are calculated using the TEF-weighted sum of DLC doses. Although the TEQ approach was developed with the intention to estimate the toxicity potential of oral exposure to DLC mixtures, the approach is also used often on the source matrix at polluted sites to compare the degrees of contamination at different locations [11].It is important to note that a potential pitfall in the TEQ approach is that most of the “intake” TEF values were derived primarily based on in vivo rodent studies and thus their direct applications to humans are not without uncertainties given very different rodent vs. human toxicokinetics of DLCs including half-lives [36]. Moreover, with the assumption of concentration addition, these TEFs were established based on single exposure experiments, not mixtures [2]. Therefore, using exposure or intake doses to calculate TEQ only treats these compounds independently, without taking into consideration the potential toxicokinetic inter-dependence or interactions of DLCs. It has been recommended that calculating body burden TEQ, by developing and applying “systemic” TEFs based on in vitro assays to biomonitoring data or model-predicted plasma or tissue concentrations of DLCs, will help in reducing the inter-species uncertainty in toxicokinetics and improving the accuracy of the TEQ approach for human risk assessment [2,37]. Therefore, developing a human DLC mixture PBPK model that can predict the internal doses of congeners will help toward this goal. 4. DiscussionPBPK modeling has been an important in silico tool in chemical safety assessment. Since humans are exposed to environmental pollutants in mixtures, it calls for PBPK models that can simulate co-exposure to multiple compounds simultaneously. When developing mixture models, it is common to start with PBPK models of single compounds and then combine them based on the known or best hypothesized compound interaction mechanisms [49,50]. Many mixture PBPK models, primarily for binary and some for higher-order mixtures, have been developed for both environmental and pharmaceutical compounds [51,52,53,54,55,56,57]. Constructing mixture PBPK models is time consuming, especially for mixtures containing diverse compounds, where pair-wise interactions can be compounds-specific. To simulate a high-order mixture, all binary interactions of the compounds in the mixture need to be characterized and appropriately captured mathematically [50,55,56,58]. Therefore, scaling up mixture PBPK models is a complex and challenging task.Yet, for compounds belonging to the same family and sharing common toxicokinetic and toxicodynamic mechanisms, mixture PBPK modeling may be simplified, as we did with the modeling framework for DLC mixtures here. The rationale hinges on the fact that DLC congeners can induce the CYP1A2 protein through activating a common nuclear receptor AHR, and CYP1A2 then sequesters and metabolizes the congeners. The DLC mixture model was based on an early human PBPK model for single TCDD exposure [25,28]. A major modification made to the Emond model is on how the DLC binding to AHR and CYP1A2 is modeled. In early days of PBPK modeling when the computation resource and power was limiting, a common practice was that fast processes, such as binding between chemicals and proteins, were often assumed at quasi-steady state relatively to other slower ADME processes, such that analytically solvable algebraic equations rather than ODEs were used to track the fast-reacting species to speed up the numerical integration. While this analytical approach is manageable for modeling one or two congeners, it quickly becomes intractable as more congeners are added to the mixture, which all compete for binding to AHR and CYP1A2. To circumvent this problem, we resorted to modeling these reversible binding processes by describing both the association and dissociation steps explicitly. While this modification may slow down the simulation to some extent, it provides the freedom to include an arbitrary number of DLC congeners into the mixture model.For mixture PBPK modeling, the amount of coding generally increases linearly with the number of compounds in the mixture since the tissue concentrations of each compound need to be tracked separately. In the present study, we utilize a vector to hold the ODEs describing the concentrations of all DLC congeners in a particular tissue rather than hard coding each congener separately. In addition, as detailed in the ODEs in Table S1, vectors are used conveniently to sum for the total amounts of DLC-bound AHR and DLC-bound CYP1A2, and for the total rates of association and dissociation for the binding processes, as done in the ODEs of common state variables AHR and CYP1A2. This vector-based approach not only saves the extra error-prone coding work, but also provides the flexibility to scale with the number of congeners in the mixture through varying the vector length accordingly. Although our model is implemented in Octave, it can be easily adapted to run in any other modern programming languages that can handle matrices, including R, Python, and Julia. Sasso et al. proposed a similar matrix-based approach for metals and nonmetal mixture modeling using a general model structure [59]. However, the chemical interactions in the mixtures still need to be handled on a case-by-case basis.The biggest uncertainty in the DLC mixture model lies in the chemical-specific parameters associated with the non-TCDD congeners, which have no validated, preexisting PBPK models similar to the TCDD model. Following the underlying principle in the TEQ approach, we assumed that the DLC congeners differ in their potency of AHR activation (i.e., KdAHR) but are equally efficacious in their transcriptional activity to induce CYP1A2. While this is likely true for most of the congeners, it may not apply to all congeners [60]. Nevertheless, this assumption of equal efficacy reduces the number of chemical-specific parameters and simplifies the coding for AHR-mediated induction of CYP1A2. The KdAHR value is not available for every congener. While some were calculated based on the REP values of receptor binding assays, many had to be estimated based on other in vitro assays, as annotated in Table S3. Another uncertainty is with CYP1A2 binding. Without prior information, we had to assume that every congener has the same binding affinity for CYP1A2 as TCDD. This assumption will clearly affect their hepatic burden predicted by the PBPK model. Finally, the hepatic elimination rate constant of each congener was estimated based on their reported half-lives and scaled with TCDD. Despite these uncertainties, it is worth noting that constructing a rigorously parameterized DLC mixture model is not the main purpose of the present study, rather, we aim to develop a mixture modeling framework for DLCs which can be further refined with better informed congener-specific parameter values in the future.As a starting point, some of the assumed or uncertain congener-specific parameters, such as binding affinities for AHR and CYP1A2 may be predicted using in silico tools, such as quantitative structure-activity relationship (QSAR) modeling followed by confirmation with measurements using binding assays [61,62]. Human cells-based in vitro assays can be utilized to measure hepatic enzymatic activities that can be extrapolated to in vivo conditions to parameterize liver clearance for the PBPK model [63,64]. Yet, validating the PBPK model for each non-TCDD DLC would require comprehensive datasets including (1) longitudinal human blood and tissue concentrations of DLCs, such as in fat and liver, and (2) related exposure data or estimation, albeit simultaneously obtaining both on the same individuals can be challenging. Structurally, the DLC PBPK model can be improved by including dose-independent pathways of elimination and reducing the dose-dependent elimination as suggested in [30]. As lipophilic compounds, the bound vs. unbound fractions of TCDD and other DLCs in the plasma can be significant, but they are not distinguished in the Emond and other dioxin models. Given the important role of plasma unbound fraction in chemical toxicokinetics, future iteration of the model should consider these factors explicitly. The present modeling framework is still steps away from application for DLC mixture risk assessment. In addition to rigorously optimizing and validating the model using parameter, exposure, and tissue burden data as described above, inter-individual variabilities need to be considered. These variabilities will include physiological variations in body weight growth and fat fractions, polymorphism in AHR abundance and affinity, and CYP1A2 induction and enzymatic activities [48,65,66]. Finally, a hierarchical Bayesian PBPK modeling approach can be utilized to integrate the population and individual-level parameter as well as exposure variabilities and uncertainties to better characterize the tissue TEQ burdens for mixture risk assessment [67].

Utilizing the mixture PBPK model, we found that due to the cross-induction of CYP1A2, the extrahepatic tissue burden of a congener in mixture exposure is always lower than in single exposure. In contrast, the hepatic burden can increase or decrease depending on the compositions of the mixture. If the two congeners are both strong CYP1A2 inducers, such as TCDD and 1,2,3,7,8-PeCDD, their hepatic burdens tend to increase due to more sequestration by CYP1A2. If one congener is significantly stronger than the other in inducing CYP1A2, such as TCDD and 1,2,3,4,6,7,8,9-OCDF, then the latter’s hepatic burden tends to decrease since more is metabolized by the induced CYP1A2 despite the increased sequestration capacity. In summary, the extrahepatic disease risk of mixtures may be lower than estimated based on linear summation of single exposures.

TEQ is an important metric that is used to evaluate the toxicity potentials and contributions of DLCs and their mixtures. Traditionally, the TEQ values are calculated for oral exposures using TEF assigned to each DLC congener. Since the TEF values were mostly obtained based on in vivo rodent assays, inter-species uncertainties exist when applied to humans. Rodents and humans are quite different in the toxicokinetics of dioxins, where dioxins are metabolized significantly slower in humans. For mixture exposures, congener–congener interactions can also play a significant role in the tissue burden and overall toxicity. Therefore, it is recommended that if human plasma or tissue concentrations are available, as in biomonitoring data, the body burden TEQ can be calculated directly. However, applying traditional “intake” TEF values to blood or tissue concentrations can lead to miscalculated mixture risk, due to the toxicokinetic differences between species and potential toxicokinetic interactions between congeners [37,68,69,70,71]. Therefore, using systemic TEF values, obtained from in vitro assays, on blood or tissue concentrations is more appropriate [37]. For example, HpCDD and a few other DLCs have significantly different systemic TEF values than the WHO-TEF values. Using the mixture PBPK model, we showed that the contributions of individual congeners to tissue TEQ can be quite different to intake TEQ, depending on the toxicokinetic interactions of the mixture. The relative contributions also vary between blood and different tissues. A congener’s contribution to the TEQ of an extrahepatic tissue, such as fat or muscle, is positively correlated with its partition coefficient for the tissue. In contrast, the binding affinity for CYP1A2 determines the contribution of a congener to the hepatic TEQ. In our mixture model, since all congeners share the same binding affinity for CYP1A2, their relative contribution to hepatic TEQ are nearly identical to blood TEQ since the amount of a congener sequestered in the liver is proportional to its plasma concentration.

In conclusion, we have developed a coding framework to facilitate PBPK modeling for DLC mixtures. The modeling tool is freely available, and as the chemical-specific parameters of congeners are better informed, the refined model can be used for DLC mixture risk assessment in the future.

留言 (0)

沒有登入
gif