Development of a Generic PBK Model for Human Biomonitoring with an Application to Deoxynivalenol

Toxicokinetic modelling provides a powerful tool in relating internal human exposure (i.e., assessed through urinary biomarker levels) to external exposure. Chemical specific toxicokinetic models are available; however, this specificity prevents their application to similar contaminants or to other routes of exposure. For this reason, we investigated whether a generic physiological-based kinetic (PBK) model might be a suitable alternative for a biokinetic model of deoxynivalenol (DON). IndusChemFate (ICF) was selected as a generic PBK model, which could be fit for purpose. Being suited for simulating multiple routes of exposure, ICF has particularly been used to relate the inhalation and dermal exposure of industrial chemicals to their urinary excretion. For the first time, the ICF model was adapted as a generic model for the human biomonitoring of mycotoxins, thereby extending its applicability domain. For this purpose, chemical-specific data for DON and its metabolites were collected directly from the literature (distribution and metabolism) or indirectly (absorption and excretion) by fitting the ICF model to previously described urinary excretion data. The obtained results indicate that this generic model can be used to model the urinary excretion of DON and its glucuronidated metabolites following dietary exposure to DON. Additionally, the present study establishes the basis for further development of the model to include an inhalation exposure route alongside the oral exposure route.


Introduction
Mycotoxins are substances produced by fungi that can contaminate plants and plantderived products during production and storage.Consequently, the general population is predominantly exposed to mycotoxins orally via the intake of food.In addition, occupational exposure can occur alongside oral exposure in certain individuals.The main route of exposure to mycotoxins in an occupational setting is via inhalation as mycotoxins can be present in organic dust.For instance, exposure can occur during the development of working routines in different types of industries (e.g., brewing and bakeries) or during the interaction with feed [1][2][3][4][5][6].
Deoxynivalenol (DON) is usually considered one of the most prevalent mycotoxins, presenting possible adverse health effects for humans and animals.Following acute exposure to DON, acute effects on the gastrointestinal tract of humans and animals have Toxins 2023, 15, 569 3 of 14 PBK models using in vitro and in silico data.This OECD guidance document provides an assessment framework for evaluating PBK models for both PBK model developers and end-users applying the models for regulatory purposes with a focus on uncertainties that underlie the model inputs and outputs [18,20].Both HttK and ICF concur with these criteria in terms of the description of the model and availability of the model code for checking equations, mass balance and/or making adaptations of the model possible.However, HttK does not include the option to simulate the distribution of the metabolites of the parent compound, whereas ICF includes this option.In addition, ICF is user-friendly and free to download.
The PBK model ICF was first developed to predict concentrations of chemicals and metabolites in alveolar air, blood and urine following environmental or occupational exposure in several populations (including males/females and children) via different routes of exposure (oral, dermal or inhalation) [14,16]).The ICF model can be applied for different exposure durations depending on the exposure route (single bolus versus repeated chronic exposure).The model contains algorithms such as Quantitative Structure-Activity Relationships (QSAR) for calculating blood:air partition coefficients, tissue:blood coefficients and the renal excretion fraction (glomerular filtration and tubular reabsorption).Besides physicochemical input, such as log (octanol:water) partition coefficients (log (Kow)) and water solubility, the model also requires an oral absorption rate constant, resorption fraction in renal tubuli, enterohepatic removal and parameters relevant for metabolism (Km and Vmax) as input.The ICF model does not take renal secretion into consideration.
Here, we adapted the ICF model as a generic modelling tool for the human biomonitoring of mycotoxins, thereby extending its applicability domain.For this purpose, chemicalspecific data for DON and its metabolites were collected directly from the literature in order to calculate metabolism rates and partition coefficients.Chemical-specific data were also collected indirectly (absorption and excretion) by fitting the ICF model to the urinary excretion data generated by the oral biokinetic model from Mengelers et al. (2019) [1].The resulting generic PBK model can be used as a tool in risk assessments relating to internal human exposure (i.e., assessed through urinary biomarker levels) to external dietary exposure.In the future, this generic PBK model can be developed further by the addition of an inhalation route, making it more suitable for an occupational population.

Results
Firstly, the structure of the DON PBK model is described based on the ICF model (Figure 1, Table 1), followed by the calculated parameters based on the publicly available literature (Figure 2, Tables 2 and 3).Lastly, the absorption rate and renal excretion rate required as input for the adapted ICF model were fitted to the urinary excretion data generated by the oral biokinetic model described by Mengelers et al. (2019) [1] (Table 4).In Figure 3, the resulting DON PBK model (red lines) was compared with the oral kinetic model [1] (black lines).Figure 4 illustrates the DON PBK model.

The DON PBK Model
Figure 1 presents the PBK model applied to DON (parent compound) and its glucuronidated metabolites.The adapted PBK model contains compartments for the liver, kidneys, heart, lungs, adipose tissue, the gastrointestinal tract from which absorption takes place (as depicted by the stomach + intestine) and lumped compartments richly and poorly perfused organs.The model shows that DON is excreted by urine and that faecal excretion is neglectable following absorption.The mathematical equations used in the DON PBK model are similar to that for ICF and can be found in the ICF manual [21].The adaption of the lumped compartments can be found in the Section 5.1 and the model code can be found in the Supplementary Materials.The calculated allometric data values for the lumped compartments poorly and richly perfused tissues can be found in Table 1.can be found in the supplementary material.The calculated allometric data values for the lumped compartments poorly and richly perfused tissues can be found in Table 1.

Calculation and Estimation of DON-Related Parameters
To adapt the ICF model as a generic modelling tool for the human biomonitoring of mycotoxins, relevant DON-related parameters were retrieved from the existing literature (i.e., organ/blood partition coefficients and metabolism-related parameters) or estimated (i.e., absorption rate constant, and excretion-related parameters through fitting the adapted ICF model to excretion amounts that were generated by the biokinetic model from Mengelers et al. (2019) [12] (see Section 5.2 and Figure 2).

Calculation and Estimation of DON-Related Parameters
To adapt the ICF model as a generic modelling tool for the human biomonitoring of mycotoxins, relevant DON-related parameters were retrieved from the existing literature (i.e., organ/blood partition coefficients and metabolism-related parameters) or estimated (i.e., absorption rate constant, and excretion-related parameters through fitting the adapted ICF model to excretion amounts that were generated by the biokinetic model from Mengelers et al. (2019) [12] (see Section 5.2 and Figure 2).

Calculation of DON-Related Parameters Based on Data from the Literature
The DON-related parameters that were calculated and/or retrieved from the literature were the organ/blood partition coefficients and the metabolism-related parameters.The (organ-specific) partition coefficients were calculated using QSARs, based on physicochemical data of DON and its metabolites (Table 2).[22]).Table 3 summarises the resulting calculated organ-specific partition coefficients for DON and its metabolites, based on the QSAR calculations used in the ICF model and the parameter values are presented in Table 2.In the ICF model, metabolism is described using the parameters Vmax (maximum velocity of metabolism) and Km (Michaelis-Menten constant).These parameters were calculated as follows.The metabolic rate retrieved from the literature was the intrinsic clearance of DON of 0.39 L/hour × kg body weight [23].In combination with the ratio of the two metabolisation rate values from Mengelers et al. (2019) [12] for the parallel generation of DON-3-GlcA and DON-15-GlcA, a metabolic rate of 0.07 and 0.32 L/hour × kg body weight, respectively, was calculated.The Vmax parameter value was estimated (as described in the Section 5.2.1.and Supplementary S1), given the large value of the Km parameter for each substance.In the last calculation step, it was assumed that the metabolism process in our model was almost linear with the ratio Vmax/Km being the metabolic rate value, in accordance with the model described in Mengelers et al. (2019) [12].

Estimation of DON-Related Parameters by Fitting
Excretion in the urine is modelled in the ICF model as a linear combination of the Glomerular Filtration Rate, the fraction of the substance in arterial blood that is dissolved in water and the fraction of the dissolved substance that is excreted in the urine, the latter depending on the Kow.However, fitting the DON PBK model with ICF's default Kow dependent renal excretion fraction value did not result in a good fit.Therefore, it was necessary to re-calibrate the renal excretion fraction.As a result, two parameters still needed to be estimated, the absorption rate constant and the renal excretion fraction.
The absorption rate constant and the renal excretion fraction were calculated by fitting the PBK model to the results of the biokinetic model developed by Mengelers et al. (2019) [12].Their values are shown in Table 4.It should be noted that, since we started with a bolus intake of value 1, and the model is completely linear, the excretion amounts can be interpreted as fractions.Remarkably, the excretion fraction of the glucuronides is above 1, indicating a renal secretion step.In Figure 3A, the excreted amounts per 1 h time interval are shown after a single oral dose, starting with the first excretion measurement after 3 h.This means that the presented values show the excreted amounts at voiding time points with 1 h time intervals, after an initial 3 h time interval.In addition, in Figure 3B, the cumulative excreted amounts over time for all absorption time points, i.e., the cumulative excretion from the starting time point t = 0, are shown.Note that the model was fitted on the interval-wise excretion amount values as depicted in Figure 3A.The latter results show that all interval-wise excreted amounts calculated this way resulted in overall cumulative excretion amounts for DON-3-GlcA and DON-15-GlcA that differed only slightly from the values that resulted from the biokinetic model from Mengelers et al. (2019) [12].Moreover, as seen in Figure 3B, the fitting procedure resulted in good fits for the lower reporting time points, where most of the DON intake was excreted as expected.Toxins 2023, 15, 569 8 of 14

Illustration of the DON PBK Model
The DON PBK model can be used to develop a human biomonitoring sampling strategy, resulting from real-life exposure situations.To illustrate this, we assumed three oral exposures during the day (e.g., during breakfast, lunch and dinner).Each of these exposures was modelled as a bolus of 1 mg DON. Figure 4 shows the blood concentration of DON and its metabolites following these real-life exposure scenarios.As expected, the blood concentrations of DON and its metabolites showed rapid changes during the day following oral intake, indicating that monitoring this time course requires an accurate sampling strategy.Mengelers et al. (2019)  [12] (black lines) and after the fitting of the model (red lines).

Illustration of the DON PBK Model
The DON PBK model can be used to develop a human biomonitoring sampling strategy, resulting from real-life exposure situations.To illustrate this, we assumed three oral exposures during the day (e.g., during breakfast, lunch and dinner).Each of these exposures was modelled as a bolus of 1 mg DON. Figure 4 shows the blood concentration of DON and its metabolites following these real-life exposure scenarios.As expected, the blood concentrations of DON and its metabolites showed rapid changes during the day following oral intake, indicating that monitoring this time course requires an accurate sampling strategy.

Discussion
Similar to the original ICF model, the PBK model that we developed for DON does not include a fractional absorption parameter (fabs).That means that the model assumes that 100% of the intake also enters the liver compartment, and thus the human body.Previous studies [12,13] have shown that only a fraction (i.e., 0.44-0.56) of DON is absorbed by the human body and that this fraction is highly heterogeneous in human populations.Therefore, to make the model useful for human biomonitoring applications, a stochastic fractional absorption parameter should be included in future modelling exercises.

Discussion
Similar to the original ICF model, the PBK model that we developed for DON does not include a fractional absorption parameter (fabs).That means that the model assumes that 100% of the intake also enters the liver compartment, and thus the human body.Previous studies [12,13] have shown that only a fraction (i.e., 0.44-0.56) of DON is absorbed by the human body and that this fraction is highly heterogeneous in human populations.Therefore, to make the model useful for human biomonitoring applications, a stochastic fractional absorption parameter should be included in future modelling exercises.
When comparing different PBK models, several similarities and dissimilarities become apparent.For example, several organs are consistently included in these models because of their specific functionality, e.g., kidney for excretion, liver for metabolism, gastrointestinal tract for absorption following oral intake, lungs for inhalation and exhalation (when applicable), and skin for dermal absorption (when applicable).PBK models differ with respect to those organs (compartments) that are relevant for internal storage or for the Toxins 2023, 15, 569 9 of 14 specific organ toxicity of the substance in the human body, e.g., brain and muscles.For these reasons, the ICF model and the Httk PBK model include a much larger set of specific organs compared to other models.However, since the in-flow and out-flow of substances in the ICF model are modelled as a linear process (i.e., proportional to the blood flow), compartments can easily be combined.The parameters to be changed were the compartment volumes and blood flows, and the partition coefficients.In cases of aggregating organs into one compartment, volumes and flow values must be added, and the partition coefficient of the new compartment equals the weighted sum of the partition coefficients of the organs.
In our study, unknown model parameters were estimated either by calculating them directly from the data retrieved from the literature, e.g., the partition coefficients and the metabolism rates, or by calculating them indirectly by fitting the adjusted ICF model, the DON PBK model, to the biokinetic model by Mengelers et al. (2019) [12] using the calculated excretion amounts (model calibration).The advantage of direct estimation is the one to one relation between the data and parameter(s), whereas a disadvantage may be the different contexts; the data value may come from an animal experiment or an in vitro environment, whereas the application involves a human in vivo situation.The model calibration process guarantees that the environments of both models are the same.Since the biokinetic model from Mengelers et al. (2019) [12] is fitted to the collected human biomonitoring data, both the biokinetic and PBK-model describe these data.The disadvantage of model calibration is the issue of parameter identifiability.It is not automatically guaranteed that the parameter values calculated by model calibration are realistic.This issue occurs when the parameters are highly correlated.In that case, rather different combinations of parameter values may lead to very similar values of the model output variables.In mathematics, this is known as the model identifiability problem.In our analyses, the problem occurred when we tried to simultaneously calculate the absorption and metabolic rates, and the excretion fractions through model calibration.The procedure was unsuccessful due to the correlations between these parameters.Unsuccessful means that we were mathematically not able to assess clear parameter values, only relatively large intervals of parameter values that all resulted in similar model fits.For this reason, we were forced to calculate one of the three unknown parameters directly from the data retrieved from the literature, allowing the other two parameters to be identified by model calibration.In our case, the metabolic rate was calculated by combining in vitro data from Faeste et al. (2018) [23] and parameter values from the biokinetic model from Mengelers et al. (2019) [12].
Two important aspects can be distinguished with respect to the parameters of the PBK-model: heterogeneity and uncertainty.Heterogeneity means that parameter values can be different over time within individuals, and also between individuals.These forms of variability are known as within-variability and between-variability, respectively.Pletz et al. (2020) [15] showed that the latter form is more important.Previously, it was shown that the between-variability of the fraction absorbed is large in human populations for DON [12,13].Regarding uncertainty, each parameter's value is uncertain, due to imperfect measurements.Pletz et al. (2020) [15] showed that errors may occur in the calculation of in vivo model parameters from results of in vitro experiments.These errors may consist of the following two terms: firstly, a variance term, indicating imprecision, and secondly, a bias term, indicating a systematic over-or underestimation.
An example of bias is the modelling of the renal excretion fraction, described by the glomerular filtration rate and tubular resorption fraction [14,21].Our results and those from Pletz et al. (2020) [15] indicate that this way of modelling underestimates the excretion of highly soluble substances such as DON.This resulted in modelling in terms of the presence of active renal secretion in our study.Another example of (possible) bias relates to the validity of the partition coefficients using the QSAR equation applied in the ICF model, which uses the log (Kow) value as an input variable.The log (Kow) values of DON and its metabolites are small, resulting in partition coefficients that are close to 1, which show that DON is highly soluble.This raises the question whether the QSAR equation is applicable for such small log (Kow) values.
As previously mentioned, heterogeneity can be distinguished between within-person changes over time, and between-person differences.The parameter values of the compartmental DON model were partly based on the parameter values of the ICF model [21], partition coefficients and metabolisation rates from the literature.Conditional to these parameter values, the absorption and excretion rates were calculated by calibration of the biokinetic model from Mengelers et al. (2019) [12].The ICF model deals with heterogeneity by providing allometric parameter values for different types of humans, e.g., normal weight versus obese, or men versus women.However, the partition coefficients and metabolisation rates from the literature and the biokinetic model state-transition rates were only presented as population mean values with 95% confidence bounds.Therefore, as a result, it is possible to describe the population heterogeneity in the DON PBK model in terms of different allometric parameter values, whereas the other parameter values are assumed to be homogenous instead of heterogeneous.In future studies, it is possible to reduce the heterogeneity by looking at sub-populations.
DON belongs to a group of mycotoxins that are called 'trichothecenes'.Chemically, trichothecenes share a tetracyclic sesquiterpenoid 12,13-epoxy-trichothec-9-ene ring system and are divided into four groups (A-D) according to different functional groups.Unlike types C and D, types A and B are frequently found as contaminants in cereals and other commodities [1].DON is a type B trichothecene and it would be worthwhile to investigate whether this generic model can also be used to study the kinetic behaviour of type A trichothecenes, like T-2/HT-2 toxin, in humans.For example, there is an important similarity between DON and T-2/HT-2 toxins regarding their elimination in mammals, as both types of trichothecenes are subject to glucuronidation and renal excretion [7,24].
The PBK model we developed based on the ICF model has been fitted for oral exposure.This model is based on the ICF model, which is suitable for several exposure routes [14], including the intake of DON through inhalation.The basis is already there, since the original ICF model already includes absorption via inhalation.However, modelling the inhalation of DON is considerably more complex.Exposure to DON through contaminated flour dust particles was already reported [5,9].DON flour dust particles are relatively large, and are in essence non-gaseous.Consequently, the ICF's model assumption of the substance being in a gaseous state where the exposure can be described by a gaseous concentration does not apply to DON.Instead, this exposure should be described in a more complex way.DON molecules could be attached to flour dust particles that differ in size.The mean flour dust concentrations differ between occupations, e.g., dough-makers and oven-workers, which is also a factor to be properly considered [25].Moreover, the average size of the particles is about 80 microns (PM 80), with the smallest about 10 microns (PM 10) and the largest 400 microns [26,27].In cases of sizes that are smaller than 80 microns [25], the particles enter the lungs through inhalation, are deposited in the lungs, and enter the human body through alveolar deposition, or trachea-bronchial deposition [27,28].In cases of larger sizes, the particles cannot enter the lungs, but instead are being coughed up, and then partly enter the human body via gastrointestinal absorption after swallowing.In addition, during risk assessment including this inhalation exposure route, toxicity of the fungi-producing DON should be taken into account next to toxicity caused by DON, since the fungi (spores) themselves can also cause health effects, such as infection, allergy and inflammation [29,30].
The current PBK-model for DON heavily relies on the biokinetic model described by Mengelers et al. (2019) [12], which was built using data from a human intervention study.This results in questions on the range of application of the PBK model.For instance, the absorption may differ in the case of DON intake through day-to-day consumption instead of an experimental dose.Yet, modelling the relation between DON intake and excretion in an everyday situation with multiple consumption moments in Norwegian volunteers revealed similar excretion and time-to-excretion parameters as obtained in the human intervention study [13].Therefore, further model validation using empirical data is necessary.

Conclusions
DON is usually considered one of the most prevalent mycotoxins, presenting possible adverse health effects for humans and animals.Consequently, robust data should be produced to properly assess the associated risk of DON exposure, such as toxicokinetic data.Modelling approaches are usually used in this context as valuable tools to obtain such data.In the present study, we have applied, for the first time, a generic PBK model by adapting the ICF model.This way, generated data could be used in the context of human biomonitoring and health risk assessment.Our results showed that for DON and its glucuronidated metabolites, a PBK model was successfully developed and its parameters estimated.Despite the associated identified limitations, this PBK model can be used for the biomonitoring of the dietary DON intake.Recognising that exposure routes other than oral are possible, such as inhalation, constitutes an important aspect, particularly in the context of occupational exposure.The present study establishes the basis for the further development of the model to include an inhalation exposure route alongside an oral exposure route.

Materials and Methods
A PBK model was developed based on the already existing PBK model, IndusChem-Fate (ICF).To calibrate the model, literature data and simulated human biomonitoring data, resulting from the biokinetic model described by Mengelers et al. (2019) [12], were used (Figure 2).All models were written in the programming language R. The source code is provided in the Supplementary S2 and the source codes for the model fitting procedures are available upon request.

General Description of the PBK Model
The ICF model [14,21] was used as a starting point.This model simulates the absorption, distribution, metabolism and excretion of chemicals.The main input data are allometric data, e.g., blood flows and organ volumes, and substance-specific data, e.g., partition coefficients and metabolism rates.The multicompartmental ICF model contains compartments for the liver, kidneys, heart, lungs, adipose tissue and a compartment for the stomach and intestine.As shown in Figure 1 (Section 2.1), the ICF model was simplified by dividing the remaining organs in two compartments, the poorly perfused (bone, muscle and skin) and the richly perfused organs (brain and bone marrow).For these combined compartments, the combined organ volumes and flows were calculated as the sum of the separate ICF volumes and flows, respectively.According to the ICF model [14,21], the same QSAR equation was used to calculate the partition coefficients for each of the organs/compartments separately.The resulting PBK model, and its kinetic parameters, for the oral intake of DON by humans are given in the Section 2 (Figure 1, and Tables 2 and 3).

DON-Related Parameters
The unknown DON-related parameters to be calculated were the organ/blood partition coefficients, the absorption rate constant, the metabolism-related parameters, and the excretion-related parameters.In Figure 2 (Section 2.2), a schematic overview is given of the parameter values that were directly derived from the literature (partition coefficients and metabolism rates) or indirectly (absorption rate constant and excretion fraction) by fitting the ICF model to the urinary excretion data as simulated by the specific biokinetic model described by Mengelers et al. (2019) [1].

Deriving DON-Related Parameters from the Literature
The partition coefficients were calculated from the Know value for DON reported in the literature [22] using the QSAR equations mentioned in Section 5.1 above.
Metabolism was described in the ICF model as a restricted process including two parameters, the Vmax parameter (the maximum velocity of metabolism) and Km (the Michaelis-Menten constant).However, in the biokinetic model [12], the metabolism of

Figure 1 .
Figure 1.Schematic overview of the DON PBK model containing the compartments liver, kidney, heart, lungs, adipose and stomach and intestine similar to the ICF model and the lumped compartments in the poorly perfused (bone, muscle and skin) and the richly perfused organs (brain, bone marrow).

Figure 1 .
Figure 1.Schematic overview of the DON PBK model containing the compartments liver, kidney, heart, lungs, adipose and stomach and intestine similar to the ICF model and the lumped compartments in the poorly perfused (bone, muscle and skin) and the richly perfused organs (brain, bone marrow).

Figure 2 .
Figure 2. Schematic overview of parameter values for DON and its metabolites that were collected directly from the literature (partition coefficients and metabolism rates) or indirectly (absorption rate constant and excretion fraction) by fitting the ICF model to the urinary excretion data as simulated by the specific biokinetic model developed by Mengelers et al. (2019) [12,14,21-23].

Toxins 2023 , 14 A:Figure 3 .
Figure 3. (A) Interval-wise renally excreted amounts since the previous reporting time are shown and (B) cumulative renally excreted amounts over time, from the starting time point t = 0 for DON and its metabolites, DON-3-GlcA (hyphened line) and DON-15-GlcA (dotted line), based on simulated data from Mengelers et al. (2019) [12] (black lines) and after the fitting of the model (red lines).

Figure 3 .
Figure 3. (A) Interval-wise renally excreted amounts since the previous reporting time are shown and (B) cumulative renally excreted amounts over time, from the starting time point t = 0 for DON and its metabolites, DON-3-GlcA (hyphened line) and DON-15-GlcA (dotted line), based on simulated data fromMengelers et al. (2019) [12] (black lines) and after the fitting of the model (red lines).

Table 1 .
Allometric data values for the poorly and richly perfused tissues compartments.

Table 1 .
Allometric data values for the poorly and richly perfused tissues compartments.

Table 2 .
Maresca et al. (2013) parameter values for DON and its main metabolites used in the model (retrieved fromMaresca et al. (2013)

Table 3 .
Calculated organ-specific partition coefficients1for DON and its metabolites.Note: In ICF, at low log (Kow), organ:blood partition coefficients are determined by the water/lipid content of the blood and the organs.Hence, all calculated partition coefficients are identical for all forms of DON.Furthermore, in ICF, the adipose tissue:blood ratio has a minimum value of 0.1. 1

Table 4 .
Estimated parameter values after fitting the PBK model to the results of the biokinetic model from Mengelers et al. (2019) [12].