Introduction

Kinetic modelling is a mathematical tool used to predict absorption, distribution, metabolism and excretion (ADME) of a substance of interest (a drug or a contaminant) in biota. Kinetic modelling is used in various research fields, as it allows the researchers to reduce animal experiments (EU 2010) and decrease costs for drug development (Knight-Schrijver et al. 2016). Moreover, pharmacokinetic (PK) modelling is a regulatory requirement for pharmaceutical companies in preclinical analysis, while risk assessment and toxicology studies implement toxicokinetic modelling (Jamei 2016). The questions solved using PK/TK modelling include prediction of local substance concentrations, translation between different exposure routes and species, dose scaling, inter-subject variability estimation (WHO and IPCS 2010) and prediction of contaminant transfer into foods of animal origin.

To establish a reliable PBTK model, experimental data are required to fit the unknown parameters, evaluate the model performance and, additionally, to measure some of the constants used in the model equations (WHO and IPCS 2010). Thus, despite one of the noble aims of PBTK modelling of substituting animal studies, experiments are often still needed to generate the model itself. Therefore, an additional effort has to be made to reduce and refine the animal experiments by, e.g., finding additional approaches to estimate the physiological model parameters. One of these parameters is the partition coefficient between model compartments for a modelled compound.

Polychlorinated dibenzodioxins (PCDDs) are dioxin-like compounds that share toxicity mechanisms with polychlorinated dibenzofurans (PCDFs) and dioxin-like polychlorinated biphenyls (dl-PCBs). Dioxins are known for their ubiquity, toxicity and bioaccumulative properties (Hites 2011). In the present study, tetrachlorodibenzo-p-dioxin (TCDD) was chosen as a model contaminant due to the highest human toxicity among the PCDDs (Van der Berg et al. 2006), availability of numerous publications on TCDD toxicity mechanisms (Denison et al. 2011), gathered information of TCDD ADME in various species (Abraham et al. 1988; Andersen et al. 1993; Emond et al. 2005; Inui et al. 2014; Kedderis et al. 1991; Molcan et al. 2017; Sakaki and Munetsuna 2010; Tritscher et al. 1992) and published PBTK models on TCDD in humans, mice and rats.

We opted for pigs as a model organism due to their significance as meat producing animal (Statistisches_Bundesamt 2018). In addition, there is an urgent need to develop PBTK models on TCDD in farm animals due to the ubiquitous distribution of TCDD in the environment and recurring inadvertent contamination of feed for farm animals with TCDD (Hites 2011). Even a background exposure to TCDD of both farm animals and humans may pose a risk to human health due to the slow excretion (Korrick et al. 2011) and has to be assessed. Another application of the here-developed PBTK model of TCDD in fattening pigs has emerged after the recent re-evaluation of the tolerable weekly intake (TWI) of dioxins and dioxin-like PCBs performed by the European Food Safety Agency (EFSA 2018). As a consequence, maximum tolerable concentrations of dioxins in edible tissues of farm animals need to be recalculated which, in turn, ideally requires the recalculation of the allowed concentrations of dioxins in feed, which is to be performed based on the predictions of PBTK models for all dioxin congeners, including the here-developed model on TCDD.

Here, we make use of the online calculator for the in silico predicted partition coefficients (Ulrich et al. 2017) and their subsequent implementation for PBTK modelling and discuss the limitations of their practical applicability. As a case example a PBTK model for orally consumed TCDD in growing pigs is developed.

Materials and methods

In silico predicted partition coefficients

Progress in environmental chemistry has made it possible to predict partition coefficient Kijk of a solute (i) among two phases (j and k) with the use of a single linear algebraic equation (Eq. S1) called polyparameter linear free energy relationship (pLFER) (Endo et al. 2013; Goss 2005). The values of these parameters (called solute descriptors and phase descriptors) can be estimated without performing animal experiments, but rather using physicochemical methods. A particular case of pLFER is the linear solvation energy relationships (LSERs). In this paper we will call Kijk in silico predicted partition coefficient of compound i between compartments j and k.

Among numerous methods to estimate partitioning (reviewed by Graham et al. 2011) of a solute between biological tissues, the online on-line LSER database was chosen due to the presence of data related to TCDD and other dioxin-like environmental pollutants, user-friendly interface of the web-calculator of partition coefficients, free availability for research aims, fast calculation time and independence of the calculation from the solvent-related in vivo experiments. Online LSER calculator of in silico predicted partition coefficients (Ulrich et al. 2017) provides as an output the value of partition coefficient between two solvents or between water and a biological tissue, where a biological tissue is approximated as a heterogeneous mixture of water, cell proteins, lipids and other biomolecules in a proportion that corresponds to this tissue (Endo et al. 2013). To obtain this value, a sufficient input for the LSER calculator is the SMILES (simplified molecular input line entry system) chemical structure of interest and the composition of a biological tissue. For the reasons of brevity, the details can be found in the supplementary section “Calculation of in silico predicted partition coefficients” and Table S1. To illustrate and to discuss the application of in silico predicted partition coefficients for modelling, we have chosen TCDD distribution in the tissues of rats, humans and pigs. Therefore, the predicted ratios of TCDD concentrations in liver, adipose tissue and blood (\({C}_{\mathrm{l}}/{C}_{\mathrm{a}}\), \({C}_{\mathrm{b}}/{C}_{\mathrm{a}}\) and \({C}_{\mathrm{b}}/{C}_{\mathrm{l}}\)) were compared with the experimentally measured values for these species.

Basic three-compartment model for TCDD kinetics in growing pigs

A basic three-compartment PBTK model of TCDD ADME was proposed based on the available experimental data on TCDD deposition in pig organs after oral consumption of TCDD (Shen et al. 2012a, b). The model compartments were blood, adipose tissue and liver. Here, the adipose tissue compartment represents all fats of the swine organism, as it constitutes more than 80% of the total lipid content of a pig (Fowler et al. 1992) and it was presumed that TCDD distributes equally into all types of fatty tissues. Moreover, the biopsy samples for TCDD concentration analysis were taken from adipose tissue, which is also reflected in the name of the compartment. Liver and adipose tissue compartments were divided into two parts (Fig. 1), where the part directly connected to the bloodstream represented the extracellular space of a tissue. The second part was the cells of the tissue. The model considered exclusively oral intake of TCDD, which upon absorption from the gastrointestinal tract (GIT) by the lymphatic route (Lakshmanan et al. 1986) entered the blood circulation and was distributed into the liver and adipose tissue. Unabsorbed TCDD was excreted. Mathematical description of TCDD absorption, distribution within the compartments and excretion consisted of a system of differential Eqs. 1–8.

Fig. 1
figure 1

Basic three-compartment PBTK model. GIT is the gastrointestinal tract. A stands for the amount of TCDD

Full size image
$$\frac{\mathrm{d}{A}_{\mathrm{T}\mathrm{C}\mathrm{D}\mathrm{D}\mathrm{c}\mathrm{o}\mathrm{n}\mathrm{s}\mathrm{u}\mathrm{m}\mathrm{e}\mathrm{d}}}{\mathrm{d}t}=\mathrm{I}\mathrm{n}\mathrm{p}\mathrm{u}\mathrm{t}\mathrm{F}\mathrm{l}\mathrm{u}\mathrm{x}$$
(1)
$$\frac{{\mathrm{d}A}_{f}}{\mathrm{d}t}={F}_{\mathrm{e}\mathrm{x}\mathrm{c}{\mathrm{r}}_{\mathrm{G}\mathrm{I}}}\times \mathrm{I}\mathrm{n}\mathrm{p}\mathrm{u}\mathrm{t}\mathrm{F}\mathrm{l}\mathrm{u}\mathrm{x} \left[+{K}_{l2f}\times {A}_{l}\right]$$
(2)
$$\frac{{\mathrm{d}A}_{\mathrm{b}}}{\mathrm{d}t}=\mathrm{I}\mathrm{n}\mathrm{p}\mathrm{u}\mathrm{t}\mathrm{F}\mathrm{l}\mathrm{u}\mathrm{x}\times \left(1-{F}_{\mathrm{e}\mathrm{x}\mathrm{c}{\mathrm{r}}_{\mathrm{G}\mathrm{I}}}\right)+{Q}_{\mathrm{a}}\times \left(\frac{{A}_{\mathrm{a}\mathrm{b}}}{{V}_{\mathrm{a}\mathrm{b}}}-\frac{{A}_{\mathrm{b}}}{{V}_{\mathrm{b}}}\right)+{Q}_{\mathrm{l}}\times \left(\frac{{A}_{\mathrm{l}\mathrm{b}}}{{V}_{\mathrm{l}\mathrm{b}}}-\frac{{A}_{\mathrm{b}}}{{V}_{\mathrm{b}}}\right)-{k}_{\mathrm{U}}\times {A}_{\mathrm{b}}$$
(3)
$$\frac{{\mathrm{d}A}_{\mathrm{u}}}{\mathrm{d}t}={k}_{\mathrm{U}}\times {A}_{\mathrm{b}}$$
(4)
$$\frac{{\mathrm{d}A}_{\mathrm{a}\mathrm{b}}}{\mathrm{d}t}={Q}_{\mathrm{a}}\times \left(\frac{{A}_{\mathrm{b}}}{{V}_{\mathrm{b}}}-\frac{{A}_{\mathrm{a}\mathrm{b}}}{{V}_{\mathrm{a}\mathrm{b}}}\right)-{\mathrm{M}\mathrm{P}\mathrm{A}}_{\mathrm{a}}\times \left(\frac{{A}_{\mathrm{a}\mathrm{b}}}{{V}_{\mathrm{a}\mathrm{b}}}-\frac{{A}_{\mathrm{a}}}{{V}_{\mathrm{a}}}\times {P}_{\mathrm{b}\mathrm{a}}\right)$$
(5)
$$\frac{{\mathrm{d}A}_{\mathrm{a}}}{\mathrm{d}t}=\mathrm{M}\mathrm{P}{\mathrm{A}}_{\mathrm{a}}\times \left(\frac{{A}_{\mathrm{a}\mathrm{b}}}{{V}_{\mathrm{a}\mathrm{b}}}-\frac{{A}_{\mathrm{a}}}{{V}_{\mathrm{a}}}\times {P}_{\mathrm{b}\mathrm{a}}\right)$$
(6)
$$\frac{{\mathrm{d}A}_{\mathrm{l}\mathrm{b}}}{\mathrm{d}t}={Q}_{\mathrm{l}}\times \left(\frac{{A}_{\mathrm{b}}}{{V}_{\mathrm{b}}}-\frac{{A}_{\mathrm{l}\mathrm{b}}}{{V}_{\mathrm{l}\mathrm{b}}}\right)-{\mathrm{M}\mathrm{P}\mathrm{A}}_{\mathrm{l}}\times \left(\frac{{A}_{\mathrm{l}\mathrm{b}}}{{V}_{\mathrm{l}\mathrm{b}}}-\frac{{A}_{\mathrm{l}}}{{V}_{\mathrm{l}}}\times {P}_{\mathrm{b}\mathrm{l}}\right)$$
(7)
$$\frac{{\mathrm{d}A}_{\mathrm{l}}}{\mathrm{d}t}={\mathrm{M}\mathrm{P}\mathrm{A}}_{\mathrm{l}}\times \left(\frac{{A}_{\mathrm{l}\mathrm{b}}}{{V}_{\mathrm{l}\mathrm{b}}}-\frac{{A}_{\mathrm{l}}}{{V}_{\mathrm{l}}}\times {P}_{\mathrm{b}\mathrm{l}}\right) \left[-{K}_{\mathrm{l}2\mathrm{f}}\times {A}_{\mathrm{l}}\right]$$
(8)

where A is the TCDD amount, Q is the blood flow rate, P is the partition coefficient, K is a rate constant, V is the compartment volume and MPA is the membrane permeability multiplied by total membrane surface area (for more details refer to Table S5). The terms in square brackets (Eqs. 2 and 8) refer to the toxicokinetic model with self-induced metabolism, which is discussed below.

Little is known about the ability of a swine organism to metabolise TCDD (Puccinelli et al. 2011; Rasmussen 2017). Here it is assumed that all possible TCDD metabolites are being formed exclusively in liver and the kinetics of their formation can be described as a first-order process dependent on CYP concentration in the liver. All TCDD that has been metabolised (chemically modified in the organism) was implicitly considered as removed from the animal.

Physiological parameters for growing pigs

The average body weight of pigs was 8.8 kg at the beginning and 127 kg at the end of the experimental period (Shen et al. 2012b). The body weight development is presented in Table S5. The growth of the pigs was modelled with a logistic equation (Eq. 9), where parameters L, k and x0 were fit (Fig. S3) with the experimental data from Table S5. The growth of the adipose tissue was also modelled with a logistic equation, based on published data for amounts of adipose tissue in relation to pig weight (Pfeiffer et al. 1984) (Table S6).

$$f\left(x\right)=\frac{L}{1+{\mathrm{e}}^{-k\times \left(x-{x}_{0}\right)}}$$
(9)

Blood weight was estimated as 6.5% of the total body weight (GV-SOLAS 2009), while liver weight was assumed to be 3% of the total body weight (Essien and Fetuga 1987). These two values were kept constant regardless of animal age, as they vary insignificantly during the lifespan. The physiological values for the kinetic models were adopted from various publications referring to pigs, which are listed in Table S9. For the parameters MPA/Qa and MPA/Ql the actual values related to pigs could not be found, so, the corresponding values for humans or rats were used as the starting values and fit.

Dose-dependent absorption (DDA) and self-induced metabolism (SIM) toxicokinetic models

In addition to the basic three-compartment model (Basic) described by Eqs. 1–8, two model extensions incorporating either dose-dependent TCDD absorption (DDA) or hepatic self-induced TCDD metabolism (SIM) were developed. The DDA model considers dose-dependent absorption of TCDD in the gastrointestinal tract. Similar to the basic model, it is described by Eqs. 1–8, but the value of parameter \({F}_{\mathrm{e}\mathrm{x}\mathrm{c}{\mathrm{r}}_{\mathrm{G}\mathrm{I}}}\) varies with the change of the TCDD exposure. The assumption of DDA is supported by previous work on TCDD kinetics in rats, while the demonstrations of the absence of enzymatic activity of CYPs on TCDD in rodents could make SIM unnecessary in the same model (Sakaki and Munetsuna 2010). The SIM model considers the specific interaction of TCDD with the Aryl hydrocarbon Receptor (AhR) (Denison et al. 2011) leading to a subsequent induction of CYP1A2 enzyme and, consequently, TCDD metabolism (Hu and Bunce 1999). The binding of TCDD to AhR is assumed to follow Michaelis–Menten saturation kinetics, so that the concentration of AhR–TCDD complex can be calculated by the following algebraic subsystem:

$$\left\{\begin{array} {l}{{C}_{\mathrm{T}\mathrm{C}\mathrm{D}\mathrm{D}}}={Conc}_{\mathrm{T}\mathrm{C}\mathrm{D}{\mathrm{D}}_{\mathrm{f}\mathrm{r}\mathrm{e}\mathrm{e}}}+{\mathrm{TCDD}}_{\mathrm{A}\mathrm{h}\mathrm{R}} \\ {\mathrm{Conc}}_{\mathrm{A}\mathrm{h}\mathrm{R}}={Conc}_{\mathrm{AhR}_{\mathrm{free}}}+TCD{\mathrm{D}}_{\mathrm{A}\mathrm{h}\mathrm{R}} \\ {K}_{\mathrm{diss}_{\mathrm{AhR}}}=\frac{{Conc}_{\mathrm{AhR}_{\mathrm{free}}}\times {Conc}_{\mathrm{T}\mathrm{C}\mathrm{D}{\mathrm{D}}_{\mathrm{f}\mathrm{r}\mathrm{e}\mathrm{e}}}} {\mathrm{T}\mathrm{C}\mathrm{D}{\mathrm{D}}_{\mathrm{A}\mathrm{h}\mathrm{R}}}\end{array}\right.$$
$$\mathrm{T}\mathrm{C}\mathrm{D}{\mathrm{D}}_{\mathrm{A}\mathrm{h}{\mathrm{R}}^{2}}-\left({C}_{\mathrm{T}\mathrm{C}\mathrm{D}\mathrm{D}}+\mathrm{C}\mathrm{o}\mathrm{n}{\mathrm{c}}_{\mathrm{A}\mathrm{h}\mathrm{R}}+{K}_{\mathrm{d}\mathrm{i}\mathrm{s}{\mathrm{s}}_{\mathrm{A}\mathrm{h}\mathrm{R}}}\right)\times \mathrm{T}\mathrm{C}\mathrm{D}{\mathrm{D}}_{\mathrm{A}\mathrm{h}\mathrm{R}}+{\mathrm{C}}_{\mathrm{T}\mathrm{C}\mathrm{D}\mathrm{D}}\times \mathrm{C}\mathrm{o}\mathrm{n}{\mathrm{c}}_{\mathrm{A}\mathrm{h}\mathrm{R}}=0$$
$$\mathrm{D}\mathrm{i}\mathrm{s}\mathrm{c}\mathrm{r}={\left(\frac{{A}_{\mathrm{l}}}{{W}_{\mathrm{l}}}+\mathrm{C}\mathrm{o}\mathrm{n}{\mathrm{c}}_{\mathrm{A}\mathrm{h}\mathrm{R}}+{K}_{\mathrm{d}\mathrm{i}\mathrm{s}{\mathrm{s}}_{\mathrm{A}\mathrm{h}\mathrm{R}}}\right)}^{2}-4\times \frac{{A}_{l}}{{W}_{l}}\times \mathrm{C}\mathrm{o}\mathrm{n}{\mathrm{c}}_{\mathrm{A}\mathrm{h}\mathrm{R}}$$
$${\mathrm{C}\mathrm{o}\mathrm{n}\mathrm{c}}_{\mathrm{A}\mathrm{h}\mathrm{R}-\mathrm{T}\mathrm{C}\mathrm{D}\mathrm{D}}=\left(\frac{{A}_{\mathrm{l}}}{{W}_{\mathrm{l}}}+\mathrm{C}\mathrm{o}\mathrm{n}{\mathrm{c}}_{\mathrm{A}\mathrm{h}\mathrm{R}}+{K}_{\mathrm{d}\mathrm{i}\mathrm{s}{\mathrm{s}}_{\mathrm{A}\mathrm{h}\mathrm{R}}}-\sqrt{\mathrm{D}\mathrm{i}\mathrm{s}\mathrm{c}\mathrm{r}}\right)/2$$
(10)

The induction of CYP1A2 enzyme in hepatocytes can be modelled with equations proposed in Dayneka et al. (1993) (Eqs. 11, 12). Similar TCDD molecular mechanism of action was used in PBTK models for rats and humans (Emond et al. 2005; Wang et al. 1997).

$$\mathrm{S}\mathrm{t}\mathrm{i}\mathrm{m}\mathrm{u}\mathrm{l}\mathrm{u}\mathrm{s}=1+\mathrm{I}\mathrm{n}A2\times \frac{{\mathrm{C}\mathrm{o}\mathrm{n}\mathrm{c}}_{\mathrm{A}\mathrm{h}\mathrm{R}-\mathrm{T}\mathrm{C}\mathrm{D}\mathrm{D}}^{h}}{{\mathrm{I}\mathrm{C}\mathrm{A}2}^{\mathrm{h}}+{\mathrm{C}\mathrm{o}\mathrm{n}\mathrm{c}}_{\mathrm{A}\mathrm{h}\mathrm{R}-\mathrm{T}\mathrm{C}\mathrm{D}\mathrm{D}}^{\mathrm{h}}}$$
(11)
$$\frac{\mathrm{d}}{\mathrm{d}t}\left({C}_{\mathrm{C}\mathrm{Y}{\mathrm{P}}_{\mathrm{i}\mathrm{n}\mathrm{d}\mathrm{u}\mathrm{c}\mathrm{e}\mathrm{d}}}\right)={k}_{\mathrm{i}\mathrm{n}}\times \mathrm{S}\mathrm{t}\mathrm{i}\mathrm{m}\mathrm{u}\mathrm{l}\mathrm{u}\mathrm{s}-{k}_{\mathrm{o}\mathrm{u}\mathrm{t}}\times {C}_{\mathrm{C}\mathrm{Y}{\mathrm{P}}_{\mathrm{i}\mathrm{n}\mathrm{d}\mathrm{u}\mathrm{c}\mathrm{e}\mathrm{d}}}$$
(12)

According to Dayneka et al. (1993), the rate of TCDD metabolism in hepatocytes is proportional to the overall concentration of the CYP1A2 enzyme, as shown in Eq. 13.

$${K}_{\mathrm{l}2\mathrm{f}}=\frac{{C}_{\mathrm{C}\mathrm{Y}{\mathrm{P}}_{\mathrm{i}\mathrm{n}\mathrm{d}\mathrm{u}\mathrm{c}\mathrm{e}\mathrm{d}}}}{{C}_{\mathrm{C}\mathrm{Y}{\mathrm{P}}_{\mathrm{b}\mathrm{a}\mathrm{s}\mathrm{a}\mathrm{l}}}}\times {K}_{\mathrm{p}\mathrm{i}\mathrm{g}}$$
(13)

All in all, the SIM model is described by Eqs. 1–13 using a constant parameter \({F}_{\mathrm{e}\mathrm{x}\mathrm{c}{\mathrm{r}}_{\mathrm{G}\mathrm{I}}}\).

Experimental data for modelling

All three models are fit to the most comprehensive set of data on TCDD ADME in growing pigs (Lenk 2007; Shen et al. 2012a, b). The authors used three groups of pigs (DG1, DG10 and DG100) fed capsules containing a mixture of PCDD/F congeners including TCDD. The experimental workflow is shown in Fig. S4. The daily consumption of TCDD by each animal group is shown in Table S7. The results of TCDD concentration measurements in subcutaneous adipose tissue and liver are presented in Table S8.

The ratio of TCDD concentrations in liver and adipose tissue on the final day of the experiment was used as the experimental values of liver–adipose tissue partition coefficient. Modelled values of partitioning between liver and adipose tissue were estimated as the ratio of fit parameters Pba and Pbl.

Parameterisation of the toxicokinetic models

The complete set of parameters for the basic, DDA and SIM models was fitted as described in “Details on the parameterisation of the toxicokinetic models” in the Supplementary Material. The performance of the DDA and SIM models was compared using the Akaike and Bayesian information criteria (AIC and BIC) (Eqs. S4 and S5).

Sensitivity analysis of the toxicokinetic models

The fit parameters were varied by 50–150% of their initial values, and the corresponding changes of the Dist value (Eq. S3) were calculated.

Evaluation of the performance of toxicokinetic models

Evaluation of model performance was carried out by comparing the results of the toxicokinetic model with the measurements during the 1999 Belgian dioxin incident.

Modelling tool

Model simulations were conducted using MATLAB 2018. The script of the DDA and SIM models is available as Supplementary Material.

Results

Applicability of calculator of partition coefficient to living systems. Comparison of in silico predicted values with experimental values

The values of the in silico predicted partition coefficients (Ulrich et al. 2017) \({K}_{\mathrm{l}\mathrm{a}}\), which describe the distribution of TCDD between liver and adipose tissue in an equilibrated state, were calculated for rats, humans and pigs and compared with the experimental data.

The required data on tissue composition are presented in Table S1. The values of in silico predicted liver–adipose tissue partition coefficients for various species are shown in Table S2 and Fig. 2.

Fig. 2
figure 2

Comparison of the in silico predicted and experimentally measured liver—adipose tissue partition coefficients. The figure illustrates in silico predicted (black line—average value, grey bars—error of indirect measurement), experimentally measured (stars for rats, circles for humans and triangles for pigs) and modelled partition coefficients (pentagons for rats and squares for pigs)

Full size image

The tissue–blood partition coefficient in humans was estimated using a PBTK model (Maruyama et al. 2002). The model was based on the experimental analysis of Japanese people’s tissues after low-level and long-term exposure to TCDD via food ingestion, which mimics the real background exposure scenario. Low-level exposure of humans to TCDD and the long elimination process of dioxins in humans are the justification to assume TCDD distribution within the human tissues to be near steady state. The work (Maruyama et al. 2002) provides the following data on partition coefficients and their standard deviations for TCDD: \({P}_{\mathrm{a}\mathrm{b}}=247\pm 78\) for partitioning between adipose tissue and blood and \({P}_{\mathrm{l}\mathrm{b}}=9.8\pm 5.7\) for partitioning between liver and blood, which results in \({P}_{\mathrm{l}\mathrm{a}}=0.039\pm 0.026\) for partitioning between liver and adipose tissue.

Partition coefficients in rats were calculated based on the experimental measurements of TCDD concentrations (Abraham et al. 1988). In this study, different TCDD dosages were injected intravenously. The resulting TCDD distribution within the tissues depends on the injected TCDD amounts and duration of the depuration phase, which determine actual TCDD concentration in liver. Data on partition coefficients for pigs shown in Fig. 2 includes in silico predicted values, experimental measurements for all three dosage groups and the values predicted by kinetic models. Similar to rats, experimentally measured \({C}_{\mathrm{l}}/{C}_{\mathrm{a}}\) depends on the TCDD concentration in liver.

The three toxicokinetic models: basic, dosage-dependent absorption (DDA) and self-induced metabolism (SIM).

A basic model describing the transfer of TCDD from feed into the tissues of growing pigs was established and parameterised by fitting to experimental values from Table S8 and as described in the Supplementary Material. With the use of a single parameter set with constant fraction of excretion Fexcr GI, it is possible to achieve an excellent fit (Dist < 1.5, eq. S3) to the experimental data of a single dose experiment only. The mode of TCDD accumulation in tissues and its clearance depends on the TCDD concentration in the feed, because the TCDD amount stored in adipose tissue by the end of the exposure period (day 91) is not linearly proportional to the amount of consumed TCDD. The basic model is thus discarded, as it cannot fit all three dosage levels simultaneously. As a consequence, additional hypotheses regarding the TCDD kinetic properties were taken into consideration, leading to more complex PBTK models.

To explain the differences in TCDD kinetics caused by variations in dose, we assumed that either the absorption of TCDD in the gastrointestinal tract of pigs is less effective for animals receiving higher doses of the contaminant (Inui et al. 2014; Kedderis et al. 1991; Sakaki and Munetsuna 2010), or that TCDD stimulates the transcription of enzymes involved in the metabolism and elimination of this chemical (Emond et al. 2004, 2005; Wang et al. 1997). The first assumption was addressed in the dose-dependent absorption model (DDA), whereas the second hypothesis was included in the model with self-induced metabolism (SIM). The best fit of the experimental data by all three PBTK models (Basic, DDA and SIM) as well as the predicted time course in blood, liver and adipose tissue are shown in Fig. 3. The basic and DDA models are effectively two-compartment models like the one for hexachlorodibenzo-p-dioxins (Adolphs et al. 2013), since blood and liver behave kinetically identically and could be potentially grouped.

Fig. 3
figure 3

Comparison of simulation results from the optimally fit basic model, dose-dependent absorption (DDA) model and self-induced metabolism (SIM) model for dose groups DG1 (0.18 ng TCDD/kg dry feed), DG10 (2.68 ng TCDD/kg dry feed) and DG100 (24.29 ng TCDD/kg dry feed). Experimental measurements are shown in red, basic model in grey, the DDA model in blue, the SIM model in green. TCDD was fed to the animals via capsule method once per day during days 15–91, which was modelled as continuous ingestion

Full size image

Fitting of experimental data with the DDA model resulted in a TCDD absorption rate in the gastrointestinal tract around 75% for low TCDD concentrations (DG1 and DG10). But the absorption dropped down to 31% when the animals ingested higher doses (DG100) (Fig. S5).

In the SIM model CYP enzymes in the liver regulate the TCDD excretion rate. The change of the CYP concentration with time is shown in Fig. S6. In the case of low levels of TCDD exposure, induction of CYP was negligible, but increased with increasing TCDD exposure (i.e. DG100).

Sensitivity analysis was performed on all parameters for the DDA and SIM models (Figs. S7 and S8). In case of the DDA model, the parameters that had the least influence on the model behaviour were MPA/Ql and MPA/Qa. For the SIM model, the least influential parameters were MPA/Ql, MPA/Qa, Fexcr GI, kU, CAhR, Kdiss AhR, InA2 and kout. Their values cannot be unambiguously determined with the current set of the experimental data.

We found that both the DDA and SIM models can precisely describe the experimental data (DistDDA = 1.44, DistSIM = 1.35). However, the relative quality of the DDA model is higher due to the lower amount of free parameters (AICDDA = − 9.4, AICSIM = − 0.2; BICDDA = − 5.5, BICSIM = 35.9).

TCDD distribution

According to both DDA and SIM models, more than 90% of the total consumed TCDD was excreted from the animals of all experimental groups (DG1, DG10 and DG100) during the 10-week-depuration period (Fig. S9). As it is reasonable for a lipophilic compound, adipose tissue is the main depot of TCDD in the swine organism: it contains more than 97% of the total TCDD stored in the animal. According to the DDA model, this result was achieved due to a fast elimination process and low absorption efficacy (Tables S11 and S12) in comparison to the SIM model, which explained fast TCDD depuration process with liver metabolism of this chemical.

TCDD half-life in pigs

The half-life of TCDD in pigs was estimated by fitting the simulated curves for relative concentrations in adipose tissue, liver, blood compartment and their sum with the exponential function (Fig. S10 and Table S13). The resulting TCDD half-life for DG1, DG10 and the apparent half-life of TCDD for DG100 were 24 days.

Discussion

To unravel the fate of TCDD during its transfer from contaminated feed into edible tissues of pigs, we developed and tested three different PBPK models (basic, DDA and SIM). Both the DDA and the SIM models were capable of describing the experimental data. Unfortunately, there is currently no experimental evidence to unambiguously choose one above the other. Comparison of AIC and BIC values for the DDA and SIM models indicate that the DDA model should be preferred. However, neither of the models can reliably predict TCDD fate immediately after the administration of the contaminant. Additional data on the TCDD concentrations in compartments during the exposure period would be required to parameterize the model and make it applicable not only for the depuration, but also for the exposure phase.

The DDA and SIM models have two common assumptions: (1) compartments relevant to the TCDD kinetics are blood, liver and adipose tissue. (2) The ratio of the adipose tissue weight to total body weight can be described with a logistic equation based on data from Pfeiffer et al. (1984). According to the DDA model, the absorption efficacy of TCDD in gastrointestinal tract is a dynamic value that depends on the TCDD concentration in feed. In contrast to this, the SIM model considered that (a) TCDD kinetics is influenced by hepatic CYPs; (b) basal CYPs concentration in the liver is constant; and (c) TCDD binding to and metabolism by liver CYP and TCDD excretion occur instantaneously (fast in comparison with the characteristic time of the system of 1 day).

In general, both the DDA and the SIM models can be used to estimate the final TCDD concentration in adipose tissue at TCDD concentrations in animal feed ranging between 0.18 and 24 ng/kg feed. For other doses, the model predictions need to be further evaluated.

Applicability of in silico predicted partition coefficients for estimation of TCDD distribution in a living organism

Partition coefficients are free parameters of the toxicokinetic models. During the fitting procedure, in silico predicted partition coefficients were chosen as the starting values. A comparison of in silico predicted partition coefficients and fit partition coefficients for DDA and SIM models is presented in Table 1. The method of calculation of partition coefficients with the use of pLFER has the standard deviation of 0.3–0.5 log units (Endo et al. 2013), so that the final standard deviation of partition coefficient between two compartments is 0.7 log units. Uncertainty of the estimation of partition coefficients caused by variability of animal tissue composition is insignificant in comparison with the input from the calculation error.

Table 1 In silico predicted partition coefficients are compared with the fit partition coefficients of the DDA and SIM models
Full size table

The in silico predicted partition coefficients coincide with the partition coefficients predicted by both toxicokinetic models (DDA and SIM). However, the DDA model predicts blood concentration of TCDD in a steady state to be one order of magnitude higher than for the SIM model. This can be partly explained by accumulation of TCDD in liver due to formation of the TCDD–AhR complex, as the value of the TCDD bound to AhR in liver constitutes 35% of the total TCDD amount in the liver for all dosage groups (derived from the constant fit values Kdiss and CAhR).

None of the models consider specific binding of TCDD to CYPs (Inui et al. 2014) or other proteins, despite the fact that basal CYP concentration in the liver is several orders of magnitude higher than TCDD concentration (Achour et al. 2011; Puccinelli et al. 2011). It is therefore possible that a significant amount of TCDD molecules form a complex with CYPs in liver, which can cause \({C}_{\mathrm{l}}/{C}_{\mathrm{a}}\) ratio to be higher than the in silico predicted value, as it is observed for the DG10 and DG100. Unfortunately, to our knowledge, there have been no measurements performed so far to evaluate Kdiss for the TCDD–CYP complex, which would allow us to estimate how much TCDD is accumulated in the liver in a bound state.

In rat, experimentally measured \({C}_{\mathrm{l}}/{C}_{\mathrm{a}}\) vary by a factor of ten and correlate positively with the exposure dose (Abraham et al. 1988). Its inconsistency with the in silico predicted values can be explained with three hypotheses. First, TCDD flows within the tissues are not in a quasi-equilibrated state. But, TCDD elimination half-life in rats is around 20 days (Andersen et al. 1993), which is very close to the one of pigs (24 days). Therefore, the quasi-equilibrated state of TCDD for a time range of one day is justified. Second, toxin-induced liver steatosis may cause redistribution of TCDD in the organism. However, there is evidence that neither chronic (Kociba et al. 1978) nor acute exposure to TCDD doses as high as 10 µg/kg bw (Wang et al. 1997) cause an increase in fat in rat hepatocytes or changes in body weight. The third hypothesis coincides with the explanation for the pigs: TCDD accumulation in the liver is caused by specific binding of TCDD to induced CYPs and/or other proteins. So, the equilibrated distribution of TCDD within rat tissues can be driven purely by partitioning processes only when the administrated dose is below 3 ng/kg bw or after a long-term depuration, so that the TCDD concentration in liver becomes less than 10 pg/g.

Availability of information related to humans allowed us to compare in silico predicted liver–adipose tissue partition coefficient with experimentally measured values. The experimental data on \({C}_{\mathrm{l}}/{C}_{\mathrm{a}}\) ratio was obtained from Japanese people exposed to a background level of TCDD, which was mainly consumed with food, as well as Seveso victims exposed to TCDD as a result of a chemical plant accident (Maruyama et al. 2002). Among all species discussed here, humans have the longest elimination time for TCDD, ranging from several months to several years. Long elimination time complemented with a long-term exposure to the contaminant causes TCDD to be in a quasi-equilibrated state. Therefore, it is reasonable to use partition coefficients to describe the distribution of this substance within the tissues. The \({C}_{\mathrm{l}}/{C}_{\mathrm{a}}\) ratio is higher for Seveso victims than for Japanese citizens (0.08 and 0.04, respectively), which can be explained by the higher exposure level. The dose of TCDD received by Seveso victims must have been sufficient to induce CYPs and cause redistribution of TCDD in the body (Hites 2011). In background-exposed Japanese people, CYPs induction had to be less influential. Japanese people exposed to the background TCDD amounts of 12.8 pg/day (Maruyama et al. 2002) had 0.13 pg TCDD/g tissue in liver (or 1.9 pg TCDD/g lipid in liver) and 3.1 pg TCDD/g tissue in fat (or 3.8 pg TCDD/g lipid in fat). The values of \({C}_{\mathrm{l}}/{C}_{\mathrm{a}}\) for both Seveso victims and Japanese people coincide with the in silico predicted values.

Taken together, the experimental data related to pigs, rats and humans indicate, that TCDD distribution within the tissues can be reliably described by in silico predicted partition coefficients from the LSER database if the TCDD concentration in liver is around 5 pg/g fat or lower. If liver tissue contains higher amounts of TCDD, then \({C}_{\mathrm{l}}/{C}_{\mathrm{a}}\) value will be higher than in silico predicted value. Two possible explanations are specific binding of TCDD to CYPs in liver and/or that TCDD is not in a quasi-equilibrated state within the tissues.

Half-life

Both experimental measurements and model simulations indicate that TCDD half-life in pigs is 24 days. This clearly contradicts a study based on a literature survey and correlation analysis predicting the half-life of TCDD in animals weighing around 100 kg, such as pigs, to be tens of years (Miniero et al. 2001). We would conclude that TCDD half-life is not an allometric feature, because the TCDD kinetic properties and the ability of liver enzymes to metabolise this contaminant are species-dependent. This makes it impossible to generalise the calculation for TCDD half-life based purely on the size of the compartments where it becomes distributed.

Evaluation of the toxicokinetic model performance using the 1999 Belgian dioxin contamination incident

We test the DDA and SIM toxicokinetic models by comparing its results to the Belgian dioxin crisis measurements. Severely dioxin-contaminated technical oil was introduced into the feed for farm animals. From January until June 1999, numerous Belgian farms were supplied with the contaminated feed, which resulted in contaminated pork products (Lok and Powell 2000). Having the information regarding the dioxin and specifically TCDD concentrations in feed, the PBTK models were able to appropriately predict contaminant concentrations in pig edible tissues.

Average values of TCDD and the total TEQ-sum of all dioxins in feed were 7 ng/kg fat (Hoogenboom et al. 2004) and 200 ng/kg fat (Bernard et al. 2002), respectively, while standard pig feed contains 4% of fat (Numata et al. 2014). We consider the case when the animals received contaminated feed during 3 months strictly before reaching the market weight. The model simulation result predicts a fat concentration of 0.5 ng/g of TCDD (both DDA and SIM models), which coincides with the experimental measurements of less than 1 ng/g fat (Bernard et al. 2002). Assuming that all dioxin congeners at 200 ng TEQ/kg fat kinetically behave similar to TCDD, the DDA and SIM models predict the contaminant concentration in adipose tissue to be 14 and 12 pg TEQ/g fat, respectively, which is comparable with the actual concentration of 6 pg TEQ/g fat for the dioxin mixture (Bernard et al. 2002) The difference between these values can be easily explained by the variation of kinetic properties among different dioxin congeners. However, the model predictions were close to the experimental values from the 1999 Belgian dioxin crisis.

Conclusion

The DDA model of TCDD transfer from contaminated feed into growing pigs derived in this work is a ready-to-use tool for risk assessment. In combination with models for other congeners, it can support risk assessment and reduce the number of animal experiments needed for decision-making.

In silico predicted partition coefficients of TCDD can reliably describe the contaminant distribution within the organism of pigs, rats and humans as long as the TCDD concentration in the liver of each species is below around 5 pg/g fat. For higher liver concentrations, the approach has severe limitations since accumulation of the contaminant into the liver becomes influential.

Porcine liver enzymes are likely involved in TCDD metabolism and dramatically speed up the elimination compared to humans. Further experimental data will help to close this gap and will aid developing more precise in silico models.