Introduction
Polychlorinated dibenzo-p-dioxins and dibenzofurans (PCDD/Fs, collectively and colloquially called ‘dioxins’) as well as polychlorinated biphenyls (PCBs) are persistent, bioaccumulative and toxic environmental contaminants. These substances may enter the animal food chain and have in the past led to feed and food contamination incidents affecting cattle(Reference Brambilla, Fochi and Falce1–Reference Piskorska-Pliszczynska, Maszewski and Mikolajczyk3), causing elevated PCDD/F and PCB levels in milk. At the same time, up to 50% of the PCDD/F and PCB human exposure, especially in infants and toddlers, can be attributed to consumption of milk and milk products(Reference Fürst, Beck and Theelen4,Reference Knutsen, Alexander and Barregård5) . Part I of this review covered the state of knowledge on data and transfer parameters from over 50 years of experimental studies; likewise, part I stressed the large variability and uncertainty found in the data and transfer parameters, explaining it in terms of factors stemming from the cow’s metabolic state and factors stemming from the contaminants physico-chemical properties(Reference Krause, Moenning and Lamp6). Based on data from experimental studies and further in silico tools, predictive toxicokinetic models are generated as an aid to modern quantitative risk assessment and risk management. The present Review focuses on providing an overview of the models that have been developed to predict the concentration of PCDD/Fs and PCBs in cow’s milk on the basis of the oral exposure of cows.
Modelling and simulation approaches have been used for a long time to describe the fate of xenobiotics (chemicals foreign to the body) across species(Reference Bonate7). Toxicokinetic models for bovines perform predictive estimations on the basis of mathematical equations that reflect the contaminants’ fate and the cow’s physiological processes. Toxicokinetic models are often based on particular animal feeding experimental data (in vivo), but unlike the raw data, they can be used to extrapolate to conditions different from the experiment. Models can in turn make use of in silico predictions of transfer subprocesses or model parameters (based on e.g. physico-chemical properties) as well as in vitro and ex vivo laboratory models(Reference Klevenhusen, Sudekum and Breves8) (providing information on tissue distribution and metabolism); this may become necessary to predict the transfer of substances for which little or no animal experimental data are available. Toxicokinetic models thus extract transfer information from those data and allow extrapolation to describe other situations of interest with a relatively small amount of additional data(Reference Numata9). Furthermore, toxicokinetic models can be used in research as a basis to verify scientific hypotheses by implementing them into models or to predict processes that cannot be captured experimentally, providing deeper insight into the fate of contaminants in the organism.
Once the model is properly parametrised and validated, it allows a user to simulate contamination scenarios and predict their outcome, either in the form of transfer parameters or using easy-to-use implementations with graphical user interfaces such as EFSA TKPlate(Reference Wiecek, Quignot and Amzal10), RIVM/WFSR FeedFoodTransfer.nl and BfR ConTrans(Reference Siemen, Schafft and Mueller-Graf11). They are used in the contexts of human and animal health risk analysis (risk assessment and risk management) as well as in ecotoxicological and environmental risk analysis. The quantitative model predictions can help risk managers simulate courses of action and make informed decisions to ensure consumer health; likewise, model predictions help risk managers decide whether it is justifiable to preserve the affected livestock. Reliable data and models can help improve risk analysis in terms of consumer protection, financial repercussions and animal welfare.
Toxicokinetic models simulate the absorption, distribution, metabolism and excretion (ADME) of a toxic substance in an animal organism. The simplest models for cows are non-physiological and make predictions for transfer parameters without attempting to explicitly mimic the transport of a substance inside the animal tissues. Others use one to two compartments (bundling many tissues) and can, despite their simplicity, be quite successful in reproducing milk concentration data. The more complex models have as many as eight compartments mimicking (groups of) tissues such as the gastrointestinal tract, liver, udder, etc.
Many models are based on feeding experiments that yield a limited number of data points for a small subset of individuals. The fate of a contaminant in one particular cow depends not only on the chemical nature of the (mix of) contaminants, but also on the factors extensively discussed in part I as influences on transfer parameters, including the metabolic state of the cow, body fat content, milk yield and matrix of the contamination source(Reference Krause, Moenning and Lamp6). A herd consists of many individuals, each of which may be in a different lactation cycle timepoint or metabolic state and have a different milk yield or body fat content(Reference Brambilla, Fochi and Falce1,Reference Fries12) . After calving, cows reduce their fat deposits and can increase the flow of contaminants into the milk, causing milk concentrations to increase up to four times the levels determined during periods of maximum weight gain. However, for non-seasonal calving herds, variability averages out these individualities, so that the contamination of such herd milk may depend more on contaminant input than on the physiology of each individual cow(Reference Tremolada, Guazzoni and Parolini13). This supports in principle the use of simpler models. At the same time, there are trends in the dairy industry (e.g. higher milk yields) that have a systematic effect on the properties of the herd. To generate and parametrise models that will be useful in future conditions, it is advisable to avoid oversimplifying the cow’s physiology. This suggests the use of more complex models that capture the physiology of the cow more closely, so that these effects may be explicitly used in predictions.
Toxicokinetic models to predict PCDD/F and PCB transfer from feed into milk can thus widely vary in their complexity. Famous is the phrase attributed to George Box: ‘All models are wrong, but some are useful’(Reference Box, Launer and Wilkinson14). Table 1 provides an aid to balance model complexity and usefulness and to choose among the different available modelling approaches that we summarise in this review. Since different purposes require different models, this review provides a summary of the available approaches and their use.
Table 1. Summary of models discussed and their respective strengths and limitations
The main focus of the review is to
-
• evaluate the availability of toxicokinetic models for all toxicologically relevant congeners (seven PCDDs, ten PCDFs and twelve dl-PCBs) as well as the indicator ndl-PCBs in terms of their applicability for risk assessment;
-
• appraise the available toxicokinetic models with respect to their capacity to make quantitative predictions for risk assessment and management.
We begin by introducing the mathematical tools used in kinetic modelling, starting with three key quantities to describe transfer: transfer rate (TR), transfer factor (TF) and biotransfer factor (BTF). These are discussed thoroughly in part I of this review in the chapter on kinetic parameters to characterise the feed-to-milk transfer behaviour(Reference Krause, Moenning and Lamp6), and we recall their mathematical definitions below in eqns. (1–3). The transfer rate (TR) describes the percentage of congener intake with the diet (mass or mole) that is excreted with the milk,
$${{\rm{TR}}\left[ \% \right] = {{{\rm{DailyExcretionViaMilk}}\left[ {{{{\rm{ng}}} \over {\rm{d}}}} \right]} \over {{\rm{DailyIntakeViaFeed}}\left[ {{{{\rm{ng}}} \over {\rm{d}}}} \right]}} \times 100\%} $$
While TRs can be calculated for any given time period during an experiment or an incident, they reach a maximum when a steady state between constant intake and output is reached. The transfer factor (TF), also known as bioconcentration factor (BCF), is a dimensionless quantity describing the ratio of the congener concentration in milk (fat) to its concentration in the feed,
$${{\rm{TF}} = {{{\rm{ConcentrationInMilkfat}}\left[ {{{{\rm{ng}}} \over {{\rm{kg}}}}} \right]} \over {{\rm{ConcentrationInFeed}}\left[ {{{{\rm{ng}}} \over {{\rm{kg}}}}} \right]}}}\;.\;$$
Lastly, the biotransfer factor (BTF) is calculated on a whole milk basis instead of milk fat, deviating from the standard for TF. Moreover, the BTF is not dimensionless and has units of time/mass, such as [d/kg], and is not restricted to an exposure from a single source (e.g. feed) but can also account for contamination through multiple pathways
$${{\rm{BTF}}\left[ {{{\rm{d}} \over {{\rm{kg}}}}} \right] = {{{\rm{ConcentrationInMilk}}\left[ {{{{\rm{ng}}} \over {{\rm{kg}}}}} \right]} \over {{\rm{TotalDailyIntake}}\left[ {{{{\rm{ng}}} \over {\rm{d}}}} \right]}}}.$$
One- and two-compartment models: mathematical motivation
In general, during contamination incidents or feeding studies with a more or less constant exposure amount or dose D [ng/d], the concentration in milk
${C_{{\rm{Milk}}}}\left[ {{\rm{ng}}/{\rm{L}}} \right]$
(usually in milk fat basis) will constantly increase and asymptotically converge towards a steady-state concentration
${C_{{\rm{max}}}}\left[ {{\rm{ng}}/{\rm{L}}} \right]\;$
(Fig. 1). This kinetic behaviour can be most simply described with a one-compartment model (Fig. 2) as done in MacLachlan (2009)(Reference MacLachlan and Bhula15), which mathematically corresponds to the differential equation
$${\rm{d}}{A_{{\rm{Cow}}}}/{\rm{d}}t = {F_{{\rm{abs}}}}D - {k_{{\rm{Milk}}}}{A_{{\rm{Cow}}}}\left( t \right)$$
Fig. 1. Hypothetical plot of the assimilation phase of a one compartment model. The system starts the assimilation phase with an initial contamination of
${C_0}\;$
and converges asymptotically against its steady state
${C_{{\rm{max}}}}$
.
Fig. 2. The one-compartment model. Here it is assumed that the cow consumes a constant amount D of a contaminant, of which
${F_{{\rm{abs}}}}$
portion gets absorbed into the ‘cow’ compartment. Finally, the cow shows a continuous excretion of the contaminant into milk at the rate
${k_{{\rm{Milk}}}}$
.
with the concentration in milk (hereafter, specifically in milk fat) thus given by
$${C_{{\rm{Milk}}}}\left( t \right) = {k_{{\rm{Milk}}}}{A_{{\rm{Cow}}}}\left( t \right)/{V_{{\rm{Milk}}}},$$
where
${k_{{\rm{Milk}}\;}}$
[1/d] is the milk excretion rate constant,
${V_{{\rm{Milk}}}}\;$
[L/d] is the milk fat yield,
${A_{{\rm{Cow}}}}\left[ {{\rm{ng}}} \right]\;$
is the amount of contaminant in the cow and
${F_{{\rm{abs}}}}$
[unitless] is mainly the fraction of dose absorbed into the cow but also accounts for all non-milk routes of elimination.
${F_{{\rm{abs}}}}$
can depend on multiple factors, such as the source of the contaminant (e.g. soil, grass, gelatine capsule) but also on the concentration itself as shown for pigs in Savvateeva et. al (2020)(Reference Savvateeva, Numata and Pieper16). From eqns. (4) and (5), one can solve for the concentration in milk fat as
$$C_{{\rm{Milk}}}(t) = {C_0}{e^{\left( { - {k_{{\rm{Milk\;}}}}t} \right)}} + {C_{{\rm{max}}}}\left( {1 - {e^{\left( { - {k_{{\rm{Milk}}\;}}t} \right)}}} \right),$$
where
${C_{0\;}}\left[ {{\rm{ng}}/{\rm{L}}} \right]$
is the initial concentration at time t = 0. Equation (6) represents the typical monoexponential behaviour of growing towards the asymptote
${C_{{\rm{max}}}} = {F_{{\rm{abs}}}}D/{V_{{\rm{Milk}}}}$
(Fig. 1) corresponding to accumulation until equilibrium in the cow. This suggests using the experimentally obtained steady-state concentration (or as an approximation, the maximum experimentally observed concentration) to estimate
${C_{{\rm{max}}}}$
, as was done in, for example, ref. (17). The value of
${C_{{\rm{max}}}}$
is the result of the dynamic equilibrium between input and elimination. The transfer factor TF can be obtained using
${\rm{TF}} = {C_{{\rm{max}}}}/{C_{{\rm{feed}}}} = {F_{{\rm{abs}}}}{W_{{\rm{Feed}}}}/\left( {{\rho _{{\rm{Milk}}}}{V_{{\rm{Milk}}}}} \right)$
, where
${W_{{\rm{Feed}}}}\left[ {{\rm{kg}}/{\rm{d}}} \right]$
is the feeding rate and
${\rho _{{\rm{Milk}}}}$
[kg/L] is the density of milk fat. Additionally, the transfer rate is given by
${\rm{TR}} = {F_{{\rm{abs}}}} \times 100\% $
.
The depuration phase commences after removing the daily exposure to contaminants and is characterised by a decrease in the amount of contaminants in the cow with a concomitant decrease in their concentration in milk over time. For the one-compartment model, the depuration behaviour of
${A_{{\rm{Cow}}}}$
is an exponential decrease to an asymptote
${C_0}$
, with the same rate
${k_{{\rm{Milk}}}}$
. For PCDD/Fs and PCBs, the depuration phase is experimentally characterised by an initial fast depuration during the first few days and a slower second depuration over several weeks and months. This biphasic behaviour is a signature of the presence of a peripheral compartment (body fat, i.e. adipose tissue) that stores contaminants and releases them slowly. During the initial fast depuration phase, mainly the portion of the contaminant in tissue that is in rapid exchange with blood is excreted via milk fat. As a result, the equilibrium of contaminants between blood and body fat is disturbed, leading to a slow remobilisation and elimination of contaminants from body fat tissues into blood and therefore into milk fat. The biphasic nature of the depuration indicates that a single rate constant is not sufficient to capture the necessary behaviour. The simplest mathematical description such a biphasic depuration phase is a two-compartment model, as shown in Fig. 3 and corresponding to the differential equation system
$${{{\rm{d}}{A_{{\rm{Cent}}}}\left( t \right)} \over {{\rm{d}}t}} = {k_{{\rm{Fat}} - {\rm{Cent}}}}{A_{{\rm{Fat}}}}\left( t \right) - \left( {{k_{{\rm{Milk}}}} + {k_{{\rm{Cent}} - {\rm{Fat}}}}} \right){A_{{\rm{Cent}}}}\left( t \right),\;$$
$${{{\rm{d}}{A_{{\rm{Fat}}}}\left( t \right)} \over {{\rm{d}}t}} = {k_{{\rm{Cent}} - {\rm{Fat}}}}{A_{{\rm{Cent}}}}\left( t \right) - {k_{{\rm{Fat}} - {\rm{Cent}}}}{A_{{\rm{Fat}}}}\left( {t} \right),$$
Fig. 3. A two-compartment model with input set to 0 and only a single output via milk. Here
${k_{{\rm{Cent}} - {\rm{Fat}}}},{k_{{\rm{Fat}} - {\rm{Cen}}}}$
are the flow rates between the compartments, and
${k_{{\rm{Milk}}}}$
is the excretion rate via milk, which is assumed to happen continuously.
with again the concentration in milk fat given by
$${C_{{\rm{Milk}}}}\left( t \right) = {k_{{\rm{Milk}}}}{A_{{\rm{Cent}}}}\left( t \right)/{V_{{\rm{Milk}}}},$$
where
${A_{{\rm{Cent}}}}$
and
${A_{{\rm{Fat}}}}$
[ng] are the amounts in the central and fat compartments and
${k_{{\rm{Fat}} - {\rm{Cent}}}},{k_{{\rm{Cent}} - {\rm{Fat}}}},{k_{{\rm{Milk}}}}\;\left[ {1/{\rm{d}}} \right]$
are the respective transition rates. During depuration phase the explicit solution for the concentration in milk fat
${C_{{\rm{Milk}}}}$
therefore has the form
$${C_{{\rm{Milk}}}}\left( t \right) = {C_A} \cdot {e^{\left( {\alpha \cdot t} \right)}} + {C_B} \cdot {e^{\left( {\beta \cdot t} \right)}},$$
where C
A
+ C
B
[ng/L] is the concentration at the beginning of the depuration phase and α and β are the elimination rate constants, which are always negative. This is the well-known biexponential decay, that is, there are two half-lives that describe the time until the concentration in the milk fat is halved in the respective phase of elimination. Inspecting eqn. (10) suggests a simple method to obtain half-lives from experimental depuration data: plot the depuration phase on a semilogarithmic scale (ln
$\left( {{C_{{\rm{Milk}}}}\left( t \right)} \right)$
) and estimate the initial slope (
$\alpha $
) and terminal (
$\beta $
) slopes (Fig. 4). This simple method has been used, for example, by Fries et al. (1973) and Brambilla et al. (2008)(Reference Brambilla, Fochi and Falce1,Reference Fries, Marrow and Gordon18) . More formally, the elimination rate constants can now be analytically determined, as they are equal to the eigenvalues of the induced transformation matrix (Supplementary Material Chapters 1 and 2), that is,
\begin{align}\alpha &= {1 \over 2}( - \sqrt {{{({k_{{\rm{Milk}}}} + {k_{{\rm{Cent}} - {\rm{Fat}}}} + {k_{{\rm{Fat}} - {\rm{Cent}}}})}^2} - 4{k_{{\rm{Milk}}}}{k_{{\rm{Fat}} - {\rm{Cent}}}}}\\
&\quad - \left( {{k_{{\rm{Milk}}}} + {k_{{\rm{Cent}} - {\rm{Fat}}}} + {k_{{\rm{Fat}} - {\rm{Cent}}}}} \right))\end{align}
Fig. 4. Hypothetical plot of the depuration phase of a two-compartment model with a linear y-axis scale (left) and logarithmic y-axis scale (right). The system starts the depuration phase with initial contaminant concentration
${C_0} + {C_1}$
and decreases double exponentially towards 0. Thereby it transitions from an almost monoexponential α depuration phase to an almost monoexponential β depuration phase.
and
\begin{align}\beta &= {1 \over 2}(\sqrt {{{({k_{{\rm{Milk}}}} + {k_{{\rm{Cent}} - {\rm{Fat}}}} + {k_{{\rm{Fat}} - {\rm{Cent}}}})}^2} - 4{k_{{\rm{Milk}}}}{k_{{\rm{Fat}} - {\rm{Cent}}}}}\\
&\quad - \left( {{k_{{\rm{Milk}}}} + {k_{{\rm{Cent}} - {\rm{Fat}}}} + {k_{{\rm{Fat}} - {\rm{Cent}}}}} \right)).\end{align}
Thus, the elimination half-lives (τ 1/2 [d]) for the depuration phase can be calculated as
$${\tau _{{1 \over 2}\alpha }} = {{\ln \left( 2 \right)} \over { - \alpha }};\;{\tau _{{1 \over 2}\beta }} = {{\ln \left( 2 \right)} \over { - \beta }}$$
Here
${\tau _{{1 \over 2}\alpha }}$
is the initial fast half-life, or ‘α-half-life’, of the contaminant, and is the result of the initial elimination from the central compartment at the start of depuration;
${\tau _{{1 \over 2}\beta }}$
is the second slower half-life of the contaminant, which is often called
$'{\rm{\beta }}$
-half-life’ or terminal-half-life, as it describes the latter and final phase of continuous elimination of the remobilised contaminant (e.g. Toutain and Bousquet-Mélou (2004)(Reference Toutain and Bousquet-Mélou19)).
Often models are proposed that comprise more than two compartments, which technically results in more than two half-lives. These additional compartments are introduced to reproduce the kinetics more precisely. However, the additional half-lives have a negligible effect on the shape of the concentration-time curve, effectively resulting in a biphasic behaviour that can be well described using only α- and
${\rm{\beta }}$
-half-lives.
Non-physiological approaches for calculating transfer parameters
Firstly, it should be noted that all three transfer parameters TR or eqn. (1), TF or eqn. (2) and BTF or eqn. (3) are conceptually similar, as all of them relate the input to the output of the contaminant (often in steady state) using different measurements of the contaminant (total amount, concentration in milk fat or concentration in milk). Therefore, it is possible to interconvert between them, as shown in part I of the review(Reference Krause, Moenning and Lamp6).
While TR, TF and BTF can be derived from experimental feeding studies or estimated from field observational data, there have been multiple attempts at predicting them for a contaminant using data from lactating cows that have not reached the steady state. One common strategy is to use experimental data from feeding studies where the cows did not reach steady-state conditions and estimate the steady-state concentration with the help of a non-physiologically based one-compartment model (Fig. 2) as presented by, for example, Connet and Webster (1987)(Reference Connett and Webster20). For this purpose, they note that in such a model the concentration in milk fat (
${C_{{\rm{Milk}}}}$
) for a given constant concentration in feed (
${C_{{\rm{Feed}}}}$
) can be described by the differential equation
$${{{\rm{d}}{C_{{\rm{Milk}}}}\left( t \right)} \over {{\rm{d}}t}} = {k_{{\rm{ass}}}}{C_{{\rm{Feed}}}} - {k_{{\rm{eli}}}}{C_{{\rm{Milk}}}}\left( t \right)$$
with the rate constants
${k_{{\rm{ass}}}},{k_{{\rm{eli}}}}$
, which can be derived from the one-compartment model (eqns. 4 and 5). These are then fitted to the experimental data (Research Triangle Institute (RTI), 2005, Appendix A(17) for more details on fitting the data). The steady-state concentration is then given by
$\displaystyle{{{k_{{\rm{ass}}}}} \over {{k_{{\rm{eli}}}}}} \cdot {C_{{\rm{Feed}}}}$
and subsequently TF =
$\displaystyle{{{k_{{\rm{ass}}}}} \over {{k_{{\rm{eli}}}}}}$
.
Other approaches are not based on animal experimental data, but rather on physical chemistry, such as Travis and Arms (1988)(Reference Travis and Arms21), who proposed a relation between BTF and the octanol-water partition coefficient
${K_{{\rm{ow}}}}$
(see also the chapter on degree of chlorination and partitions coefficients in part I of the review(Reference Krause, Moenning and Lamp6)) using linear least-squares fitting such that
$${\log _{10}}\left( {{\rm{BTF}}} \right) = - 8.085 + 0.992\;{\log _{10}}\left( {{K_{{\rm{ow}}}}} \right).$$
A geometric mean approach was discussed by Birak et al. (2001)(Reference Birak, Yurk and Adeshina22) as an alternative to the linear square approach. The idea of Travis and Arms was further developed in RTI(17) as this method becomes increasingly inaccurate for higher values of
${\log _{10}}\left( {{K_{{\rm{ow}}}}} \right)$
, which is especially relevant here, as PCDD/Fs and PCBs have rather high
${\log _{10}}\left( {{K_{{\rm{ow}}}}} \right)$
values. Therefore, they fitted the BTF data with help of a second-order polynomial resulting in
$${\log _{10}}\left( {{\rm{BTF}}} \right) = - 3.56 + 1.07{\log _{10}}\left( {{K_{{\rm{ow}}}}} \right) - 0.099{\log _{10}}\left( {{K_{{\rm{ow}}}}} \right)_{}^2.$$
For a more in-depth comparison of such methods using some form of fitting of
${\log _{10}}\left( {{K_{{\rm{ow}}}}} \right)$
, see Takaki et al. (2015)(Reference Takaki, Wade and Collins23). Dowdy et al. (1996)(Reference Dowdy, McKone and Hsieh24) took a slightly different approach, as they noted that the experimentally derived
${\log _{10}}\left( {{K_{{\rm{ow}}}}} \right)\;$
values for the same contaminant can vary widely depending on the method used(Reference Mackay, Shiu and Lee25) and furthermore that the metabolisation rate of the contaminant should also be taken into account. Therefore, they developed a quantitative structure–activity relationship (QSAR) method, based on the Randic branching index(Reference Randic26) of a given contaminant’s molecular structure to derive the ‘normal path first-order Molecular Connectivity Index’,
${}_{}^1{\chi _{pc}}$
. They presumed
${}_{}^1{\chi _{pc}}$
determines the lipophilicity and the metabolic stability of the contaminant. Hence, they effectively used
${}_{}^1{\chi _{pc}}\;$
instead of
${\log _{10}}\left( {{K_{{\rm{ow}}}}} \right)$
for linear square fitting, resulting in a formula that depends only on
${}_{}^1{\chi _{pc}}$
to predict the BTF of a contaminant, that is,
$${\log _{10}}\left( {{\rm{BTF}}} \right) = - 5.879 + 0.421{}_{}^1{\chi _{pc}}.\;$$
Models based on physiological approaches
The first physiologically based pharmacokinetic/toxicokinetic models (PBPK/TK) for PCDD/Fs and PCBs in lactating cows were published by Derks et al. in 1994(Reference Derks, Berende and Olling27) and McLachlan as early as 1992(Reference McLachlan28). They used different modelling approaches, both of which are still in use today. Additional models for lactating cows focusing on general lipophilic/hydrophobic contaminants with similar physico-chemical properties for PCDD/Fs and PCBs were proposed by different authors and have since been used for PCDD/Fs and PCBs. These models are discussed below.
The classical PBTK approach by Derks
The most prominent model was published by Derks et al. (1994)(Reference Derks, Berende and Olling27). It is a classical physiologically based toxicokinetic (PBTK) model that describes the ADME processes of a contaminant in an organism while taking into account various physiological and physico-chemical factors of an individual lactating cow. In a classical PBTK approach, the contaminant is distributed from one compartment to another, whereby the concentration-driven rate terms depend on several characteristics of the animal and contaminant, as well as on the compartments themselves. All the rate terms are combined into a system of mass balance equations that describes the amount of contaminant in each compartment over time, as well as the outflow in the form of metabolised contaminant and milk excretion. The PBTK model of Derks et al. (1994)(Reference Derks, Berende and Olling27) consists of six compartments (Fig. 5): blood, which connects all compartments; liver, in which metabolic degradation occurs; udder (represented only by udder fat), from which continuous excretion via milk fat occurs; body fat as peripheral storage compartment; and the remaining organs, which are divided into slowly (e.g. muscle, skin, bones) and richly blood-perfused (main internal organs except liver, e.g. kidney and gastrointestinal tract). The substance enters the system via the liver, so this model takes first-pass kinetics into account. The distribution between blood and each tissue compartment depends on three variables: the blood flow
${Q_i}$
[L/d], the compartment volume
${V_i}$
[L] (both of which depend on the physiology of the individual cow) and the partition coefficient
${P_i}$
[unitless], which reflects the physico-chemical properties of the contaminant by describing the tissue–blood ratio of the contaminant in equilibrium. In addition, it is assumed that all transitions between the compartments are blood flow limited, except for the fat compartment, which is diffusion limited and is taken into account by multiplying the blood flow
${Q_{\rm{F}}}\;\left[ {{\rm{L}}/{\rm{d}}} \right]$
by a constant
${F_Q}$
≤ 1. Blood flow limited means that it is assumed that the amount of blood flow into the tissue is the limiting factor in the exchange of substances, that is, the blood within the tissue is immediately in steady state with the tissue. Diffusion limited means that we assume the limiting factor is the exchange of contaminant from blood to tissue and is not instantaneous (and by definition not instantly in steady state). The liver metabolism is accounted for with a first-order rate constant
${k_{{\rm{met}}}}$
[1/d] and the proportion that is absorbed from the GIT via first-pass into the liver is accounted for by a redefined
${F_{{\rm{abs}}}}$
[unitless]. The milk fat yield is now labelled
${\rm{C}}{{\rm{L}}_{{\rm{Milk}}}}$
[L/d], instead of the synonymous
${V_{{\rm{Milk}}}}$
from previous models to underline the fact that this variable serves the function of kinetic clearance of contaminant through the removal of udder fat (identical in concentration to milk fat). The assumption of continuous lactation throughout the day is made. The resulting differential equation system is
$${{{\rm{d}}{A_{{\rm{Blood}}}}} \over {{\rm{d}}t}} = \sum_{i\epsilon T} ({{Q_{i}}{A_{i}} \over {{V_{i}}{P_{i}}}} - {{Q_{i}}{A_{Blood}} \over {{V_{Blood}}}}) + {F_Q}{{{Q_{{\rm{Fat}}}}{A_{{\rm{Fat}}}}} \over {{V_{{\rm{Fat}}}}{P_{{\rm{Fat}}}}}} - {F_Q}{Q_{{\rm{Fat}}}}{{{A_{{\rm{Blood}}}}} \over {{V_{{\rm{Blood}}}}}}$$
Fig. 5. Schematic depiction of the original six-compartment model derived in Derks et al. (1994)(Reference Derks, Berende and Olling27). Here,
${Q_i}_{}$
[L/d] stand for the blood flow rate into/out of the compartment i,
${P_i}$
[unitless] is the (compartment i)/blood partition coefficient and
${V_i}$
[L] is the volume of compartment i. The compartments i are liver, richly perfused tissues, slowly perfused tissues, udder, fat and blood. For fat we have an additional constant
${F_Q}$
[unitless] accounting for the fact that this compartment is diffusion limited. The input into this model happens continuously through the liver with D [ng/d] being the dose of contaminant fed to the cow daily and
${F_{abs}}$
the fraction absorbed into the system. Metabolism of the contaminant takes place in the liver at the rate
${k_{met}}$
[1/d]. Additionally, the contaminant is excreted in the udder via milk proportional to the amount of milk fat excreted
$C{L_{Milk}}\left[ {L/d} \right]$
.
with T = {Slow, Rich, Udder, Liver},
$${{{\rm{d}}{A_{{\rm{Fat}}}}} \over {{\rm{d}}t}} = {F_Q}{Q_{{\rm{Fat}}}}{{{A_{{\rm{Blood}}}}} \over {{V_{{\rm{Blood}}}}}} - {F_Q}{{{Q_{{\rm{Fat}}}}{A_{{\rm{Fat}}}}} \over {{V_{{\rm{Fat}}}}{P_{{\rm{Fat}}}}}},$$
$${{{\rm{d}}{A_{{\rm{Liver}}}}} \over {{\rm{d}}t}} = {Q_{{\rm{Liver}}}}{{{A_{{\rm{Blood}}}}} \over {{V_{{\rm{Blood}}}}}} + {F_{{\rm{abs}}}}D - {{{Q_{{\rm{Liver}}}}{A_{{\rm{Liver}}}}} \over {{V_{{\rm{Liver}}}}{P_{{\rm{Liver}}}}}} - {k_{{\rm{met}}}}{A_{{\rm{Liver}}}},$$
$${{{\rm{d}}{A_{{\rm{Rich}}}}} \over {{\rm{d}}t}} = {Q_{{\rm{Rich}}}}{{{A_{{\rm{Blood}}}}} \over {{V_{{\rm{Blood}}}}}} - {{{Q_{{\rm{Rich}}}}{A_{{\rm{Rich}}}}} \over {{V_{{\rm{Rich}}}}{P_{{\rm{Rich}}}}}},$$
$${{{\rm{d}}{A_{{\rm{Slow}}}}} \over {{\rm{d}}t}} = {Q_{{\rm{Slow}}}}{{{A_{{\rm{Blood}}}}} \over {{V_{{\rm{Blood}}}}}} - {{{Q_{{\rm{Slow}}}}{A_{{\rm{Slow}}}}} \over {{V_{{\rm{Slow}}}}{P_{{\rm{Slow}}}}}},$$
$${{{\rm{d}}{A_{{\rm{Udder}}}}} \over {{\rm{d}}t}} = {Q_{{\rm{Udder}}}}{{{A_{{\rm{Blood}}}}} \over {{V_{{\rm{Blood}}}}}} - {{{Q_{{\rm{Udder}}}}{A_{{\rm{Udder}}}}} \over {{V_{{\rm{Udder}}}}{P_{{\rm{Udder}}}}}} - {{{\rm{C}}{{\rm{L}}_{{\rm{Milk}}}}{A_{{\rm{Udder}}}}} \over {{V_{{\rm{Udder}}}}}}.$$
The concentration in milk fat is thus given by
$${C_{{\rm{Milk}}}} = {{{A_{{\rm{Udder}}}}} \over {{V_{{\rm{Udder}}}}}}.$$
Different methods have been used to obtain model parameters. Especially notable is the calculation of the partition coefficients
${P_i}$
, which was discussed in detail by Derks (1994)(Reference Derks, Berende and Olling27) and van Eijkeren (1998)(Reference van Eijkeren, Jager and Sips29), as we summarise below. Blood flow and organ volume were directly derived from experimental data and
${k_{{\rm{met}}}},\;{F_Q}$
and
${F_{{\rm{abs}}}}$
fitted to experimental data with numerical methods.
The determination of partition coefficients
${P_i}$
was done differently in Derks (1994)(Reference Derks, Berende and Olling27) and van Eijkeren (1998)(Reference van Eijkeren, Jager and Sips29). While Derks estimated the partition coefficient
${P_i}$
by dividing the tissue concentration of the contaminant by the blood concentration at the end of the study, van Eijkeren et al. (1998) estimated the partition coefficients using the
${K_{{\rm{ow}}}}$
of the contaminant and various generic tissue component fractions(Reference van Eijkeren, Jager and Sips29). But in MacLachlan 2009(Reference MacLachlan30) it is noted that the latter method produces almost indistinguishable values for contaminants with log
$\left( {{K_{{\rm{ow}}}}} \right)$
> 3; since all PCDD/Fs and PCBs fulfill this property, the method incorrectly predicts the same partition coefficients and therefore almost identical distribution for each congener among the compartments. It is thus recommended to use better methods to predict the partition coefficient for PCDD/Fs and PCBs, for example, Graham et al. (2011)(Reference Graham, Walker and Jones31) or Endo et al. (2013)(Reference Endo, Brown and Goss32).
Derks et al. originally used their model to describe the dynamics of 2,3,7,8-TCDD in lactating cows(Reference Derks, Berende and Olling27), and other authors have since adapted it to describe other lipophilic contaminants. More recent studies(Reference Hoogenboom, Zeilmaker and van Eijkeren2,Reference van Eijkeren, Jager and Sips29,Reference Freijer, van Eijkeren and Sips33) combined the udder fat and blood compartments into one blood compartment (Fig. 6), as the udder has a high blood flow
${Q_{{\rm{Udder}}}}$
compared with its small volume
${V_{{\rm{Udder}}}}$
, and therefore is almost instantly in equilibrium with the blood(Reference van Eijkeren, Jager and Sips29); this modification introduces a milk/blood partition coefficient,
${P_{{\rm{Milk}}}}$
, which is conceptually similar to the now missing compartment udder/blood partition coefficient, that is, the concentration in milk fat is then given by
$${C_{{\rm{Milk}}}} = {{{A_{{\rm{Blood}}}}} \over {{V_{{\rm{Blood}}}}}}{P_{{\rm{Milk}}}}.$$
Fig. 6. Schematic depiction of the modified Derks (1994)(Reference Derks, Berende and Olling27) model with the udder included in the blood compartment. Here
${Q_i}_{}$
[L/d] stands for the blood flow rate into/out of the compartment i,
${P_i}\left[ {unitless} \right]$
is the partition coefficient between blood and compartment i and
${V_i}$
[L] is the volume of compartment i. The compartments i are liver, richly perfused tissues, slowly perfused tissues, body fat, blood and milk. For body fat, there is an additional constant
${F_Q}$
[unitless] accounting for the fact that this compartment is diffusion limited. The input into this model happens continuously through liver with daily contaminant dose D [ng/d] and fraction absorbed
${F_{{\rm{abs}}}}$
[unitless]. Metabolism of the contaminant takes place in the liver at the rate
${k_{{\rm{met}}}}$
[1/d]. Additionally, the contaminant from the blood can be excreted via milk proportional to the amount of milk fat excreted
${\rm{C}}{{\rm{L}}_{{\rm{Milk}}}}\;\left[ {{\rm{L}}/{\rm{d}}} \right]$
.
This only changes the equation system slightly (Supplementary Material Chapters 1–4 and Equation S10). Additionally, it is possible to use this model for beef cattle or calves (non-lactating) by also removing the udder compartment and setting
${\rm{C}}{{\rm{L}}_{{\rm{Milk}}}} = 0$
and therefore having no milk excretion(Reference Freijer, van Eijkeren and Sips33,Reference Bogdal, Züst and Schmid34) . Such a model without milk excretion had already been used by Leung et al. (1990)(Reference Leung, Paustenbach and Murray35) for the description of TCDD kinetics in rats.
The fugacity approach by McLachlan
A different approach was proposed in McLachlan (1992)(Reference McLachlan28): a fugacity model to describe the dynamics of hydrophobic contaminants in a lactating cow; this was further developed in Rosenbaum et al. (2009)(Reference Rosenbaum, McKone and Jolliet36) and Tremolada et al. (2014)(Reference Tremolada, Guazzoni and Parolini13). Such models are based on more general multimedia fugacity models (MFM) from environmental chemistry(Reference Parnis and Mackay37). MFMs are often used to describe the fate of chemical contaminants across whole environmental compartments, and specifically the rates at which they move between phases. The transfer rate is proportional to the fugacity difference between the source and destination phases. The basis of the model is the mass balance equations for each phase including fugacities, fluxes and amounts, in this case, applied to a single organism with inputs and outputs. The fugacity (
$f$
) has units of pressure [Pa].
A key concept is the fugacity capacity (
${Z_m}$
) [mol/(m3Pa)], which is conceptually the capacity of compartment m (a phase) to absorb a solute (contaminant). The fugacity capacities
${Z_m}$
are calculated with the equilibrium partition coefficients of the chemicals, Henry’s law and other physico-chemical equations. The concentration
${C_m}$
of a chemical in compartment m is given by
$${C_m} = {Z_m}{f_m}.$$
Note that conceptually
${Z_m}$
is similar to the partition coefficient of the classical PBTK approach in the sense that
$${{{Z_m}} \over {{Z_i}}} = {P_{mi}}\left( { = {{{C_{m,ss}}} \over {{C_{i,ss}}}}} \right)$$
as in equilibrium among compartments
${f_{m,ss}} = {f_{i,ss}}$
holds true.
The transport coefficients
$D$
[mol/(Pa·d)] describe processes, such as advective transport (of a substance by bulk motion, e.g. the ingestion of a contaminant with feed), transformation (e.g. metabolisation) and diffusion.
$D$
is defined for advective processes as the product of a volume flow rate [m3/d] and a fugacity capacity
$Z$
[mol/(m3Pa)];
$D$
is defined for diffusive processes as the product of a conductance [m/d], an interface area [m2] and a fugacity capacity; and for transformation
$D$
is defined as the product of a rate constant [1/d], a compartment volume
$V$
[m3] and a fugacity capacity [mol/(
${m^3}$
Pa)](Reference McLachlan38). One conceptual core difference to the classical PBTK approach is that blood flow is not considered a limiting factor for the distribution of the contaminant, that is, purely diffusion-limited kinetics are assumed.
The MFM proposed by McLachlan consists of three compartments (Fig. 7): the digestive tract as the entry point into the system; the blood, which distributes the substance throughout the body; and finally, body fat as the storage compartment. The substance can be excreted either from the digestive tract via the faeces or from the blood via milk. In addition, the substance can also be metabolised in the blood compartment or the digestive system.
Fig. 7. Schematic depiction of the fugacity model proposed by McLachlan (1994)(Reference McLachlan38). Here
${D_{{\rm{Dig}} - {\rm{Blood}}}}$
[mol/(Pa·d)] and
${D_{{\rm{Blood}} - {\rm{Fat}}}}$
[mol/(Pa·d)] are the transport coefficients between the compartments. The input into the system is given by dose [mol/d] into the digestive tract. Excretion can happen via faeces out of the digestive tract or via milk out of the blood with transport coefficients
${D_{{\rm{Exc}}}}$
[mol/(Pa·d)] and
${D_{{\rm{Milk}}}}$
[mol/(Pa·d)], respectively. Additionally, in both these compartments, the contaminant can be metabolised with transport coefficients
${D_{{\rm{Dig}} - {\rm{Meta}}}}$
[mol/(Pa·d)] and
${D_{{\rm{Blood}} - {\rm{Meta}}}}$
[mol/(Pa·d)], respectively.
An additional assumption is made, namely that the system is always in a ‘pseudo-equilibrium’, that is, from the knowledge of the fugacity in one compartment, all other fugacities can be calculated; importantly, only the fat compartment acts dynamically. This results in a mass balance equation system of the form
$${\rm{Dose}} = {D_{{\rm{Exe}}}}{f_{Dig}} + {D_{{\rm{Dig}} - {\rm{Blood}}}}\left( {{f_{{\rm{Dig}}}} - {f_{{\rm{Blood}}}}} \right) + {D_{{\rm{Dig}} - {\rm{Meta}}}}{f_{{\rm{Dig}}}},$$
\begin{align}{D_{{\rm{Dig}} - {\rm{Blood}}}}\left( {{f_{{\rm{Dig}}}} - {f_{{\rm{Blood}}}}} \right) &= {D_{{\rm{Milk}}}}{f_{{\rm{Blood}}}} + {D_{{\rm{Blood}} - {\rm{Fat}}}}\left( {{f_{{\rm{Blood}}}} - {f_{{\rm{Fat}}}}} \right)\\
&\quad + {D_{{\rm{Blood}} - {\rm{Meta}}}}{f_{{\rm{Blood}}}},\end{align}
$${D_{{\rm{Blood}} - {\rm{Fat}}}}\left( {{f_{{\rm{Blood}}}} - {f_{{\rm{Fat}}}}} \right) = {{{\rm{d}}\left( {{V_{{\rm{Fat}}}}{Z_{{\rm{Fat}}}}{f_{{\rm{Fat}}}}} \right)} \over {{\rm{d}}t}}.$$
And therefore the concentration in milk fat is given by
$${C_{{\rm{Milk}}}} = {{{D_{{\rm{Milk}}}}{f_{{\rm{Blood}}}}} \over {{\rm{C}}{{\rm{L}}_{{\rm{Milk}}}}}},$$
where
${\rm{C}}{{\rm{L}}_{{\rm{Milk\;}}}}\left[ {{\rm{mol}}/{\rm{d}}} \right]$
is the amount of milk fat excreted each day.
Owing to the pseudo-equilibrium assumption, there is only one linear differential equation, so the McLachlan (1994)(Reference McLachlan38) model mathematically behaves as a one-compartment model, thereby inducing only one half-life (no biphasic behaviour). With the help of various data sets, McLachlan was able to create formulas for all non-metabolic transport coefficients that depend only on the
${K_{{\rm{ow}}}}$
value and Henry’s law H of the contaminant. To do that, it was assumed that the contaminant has to pass through a water and lipid layer to change from one compartment to another. For the metabolic transport coefficients
${D_{{\rm{Blood}} - {\rm{Meta}}}}$
and
${D_{{\rm{Dig}} - {\rm{Meta}}}}$
, no satisfactory data were available and the respective factors were set to 0 in the simulations.
A similar approach with the same three compartments was later used in Tremolada et al. (2014)(Reference Tremolada, Guazzoni and Parolini13). Here, the pseudo-equilibrium assumption was dropped so that a biexponential behaviour can be reproduced; the volumes of all three compartments (and not only the volume
${V_{{\rm{Fat}}}}$
of the fat compartment) were additionally considered. Furthermore, the input parameter Dose is also described in terms of fugacity, that is, Dose =
${D_{{\rm{Grass}}}}{f_{{\rm{Grass}}}} + {D_{{\rm{Feed}}}}{f_{{\rm{Feed}}}} + \;{D_{{\rm{Soil}}}}{f_{{\rm{Soil}}}}$
. This results in the differential equation system
\begin{align}&
{{{\rm{d}}{f_{{\rm{Dig}}}}} \over {{\rm{d}}t}}{\rm{ }} =\\
&\qquad{\rm{ }}{{{D_{{\rm{Grass}}}}{f_{{\rm{Grass}}}} + {D_{{\rm{Feed}}}}{f_{{\rm{Feed}}}} + \;{D_{{\rm{Soil}}}}{f_{{\rm{Soil}}}} + {D_{{\rm{Blood}} - {\rm{Dig}}}}{f_{{\rm{Blood}}}} - \left( {{D_{{\rm{Exc}}}} + {D_{{\rm{Dig}} - {\rm{Met}}}}} \right){f_{{\rm{Dig}}}}^,} \over {{V_{{\rm{Dig}}}}{Z_{{\rm{Dig}}}}}}\end{align}
\begin{align}&
{{{\rm{d}}{f_{{\rm{Blood}}}}} \over {{\rm{d}}t}}{\rm{ }} =\\
&\qquad {\rm{ }}{{{D_{{\rm{Blood}} - {\rm{Dig}}}}\left( {{f_{{\rm{Dig}}}} - {f_{{\rm{Blood}}}}} \right) + {D_{{\rm{Blood}} - {\rm{Fat}}}}({f_{{\rm{Fat}}}} - {f_{{\rm{Blood}}}}) - \left( {{D_{{\rm{Milk}}}} + {D_{{\rm{Blood}} - {\rm{Meta}}}}} \right){f_{{\rm{Blood}}}}^,} \over {{V_{{\rm{Blood}}}}{Z_{{\rm{Blood}}}}}}\end{align}
$${{{\rm{d}}{f_{{\rm{Fat}}}}} \over {{\rm{d}}t}} = {{{D_{{\rm{Blood}}}}\left( {{f_{{\rm{Blood}}}} - {f_{{\rm{Fat}}}}} \right)} \over {{V_{{\rm{Fat}}}}{Z_{{\rm{Fat}}}}}}.$$
And the concentration in milk fat is again given by
$${C_{{\rm{Milk}}}} = {{{D_{{\rm{Milk}}}}{f_{{\rm{Blood}}}}} \over {{\rm{C}}{{\rm{L}}_{{\rm{Milk}}}}}}.$$
Here, the transport coefficients
${D_i}$
were derived similarly as in McLachlan (1994)(Reference McLachlan38). Additionally, the metabolic rate constants were calculated under the assumption that they are the sole reason for the discrepancy between measured excretion via milk + faeces and input of contaminants. Furthermore, it is assumed that the metabolic rate is also proportional to the lipid volume of the compartment and its fugacity capacity, that is,
${D_{i - {\rm{Meta}}}} = {k_i}{V_i}{Z_{{\rm{oct}}}}$
for a fitted
${k_i}$
, where
${Z_{{\rm{oct}}}}\;$
is the fugacity capacity of octanol.
In this context, the CKow dynamic model of transfer to meat and milk for lipophilic contaminants proposed by Rosenbaum et al. (2009)(Reference Rosenbaum, McKone and Jolliet36) should be mentioned. At its core, CKow is a three-compartment model of the same structure as McLachlan (1994)(Reference McLachlan38), where the transition terms between the compartments are also derived similarly to McLachlan’s, but instead of the fugacities of each compartment, they work with concentration of the contaminant, thereby eliminating the need of transforming fugacities into concentration in practical applications.
Generalised models for the transfer of lipophilic contaminants into milk
Generalised models for the transfer of lipophilic contaminants into cow’s milk can also be used for PCDD/Fs and PCBs. One such generalised model was developed in MacLachlan (2009)(Reference MacLachlan30). This is a classic PBTK model with eight compartments (Fig. 8), which is similar in structure to the model developed by Derks in 1994, but with two major differences. The first difference is that the remaining tissues are not divided into poorly and richly perfused, but into muscle, kidney and other tissue compartments. The other difference is the addition of a rumen compartment, which creates a gradual passage (exponentially distributed input) to the intestine following first-pass kinetics via the liver; thereafter, the contaminant follows liver first-pass metabolism. While these generalisations make the model widely applicable, for PCDD/Fs and PCBs (because of their long half-lives), rumen lag before liver first-pass effect may not be so important to model explicitly(Reference MacLachlan30).
Fig. 8. Schematic depiction of the general eight-compartment model derived in MacLachlan (2009)(Reference MacLachlan30). Here
${Q_i}$
[L/d] stand for the blood flow rates into/out of the compartment i,
${P_i}$
[unitless] is the (compartment i)/blood partition coefficient,
${V_i}$
[L] is the volume of compartment i. The compartments i are liver, richly perfused tissues, slowly perfused tissues, udder, body fat, blood and milk. For body fat there is an additional constant
${F_Q}$
accounting for the fact that this compartment is diffusion limited. The input into this model happens continuously into the rumen, with D [ng/d] being the dose of contaminants in feed. From the rumen, the fraction
${F_{{\rm{abs}}}}$
[unitless] of contaminant gets absorbed at the rate
${k_a}\;\left[ {1/{\rm{d}}} \right]\;$
into the main part of the system; the rest is excreted via the faeces. Metabolism of the contaminant takes place in the liver with the clearance
${\rm{C}}{{\rm{L}}_{{\rm{Liver}}}}$
[1/d]. Additionally, the contaminant can be excreted from the udder via milk, proportional to the amount of milk fat excreted
${\rm{C}}{{\rm{L}}_{{\rm{Milk}}}}\;\left[ {{\rm{L}}/{\rm{d}}} \right]$
.
The model can be written as the differential equation system
$${{{\rm{d}}{A_{{\rm{Rumen}}}}} \over {{\rm{d}}t}} = \left( {{F_{{\rm{abs}}}} - 1} \right){k_a}{A_{{\rm{Rumen}}}} + {\rm{Dose}},$$
\begin{align}{{{\rm{d}}{A_{{\rm{Liver}}}}} \over {{\rm{d}}t}}{\rm{ }} &= {Q_{{\rm{Liver}}}}{{{A_{{\rm{Blood}}}}} \over {{V_{{\rm{Blood}}}}}} + {F_{{\rm{abs}}}}{k_a}{A_{{\rm{Rumen}}}} - {{{Q_{{\rm{Liver}}}}{A_{{\rm{Liver}}}}} \over {{V_{{\rm{Liver}}}}{P_{{\rm{Liver}}}}}}\\
&\quad - {{{k_{{\rm{met}}}}{A_{{\rm{Liver}}}}} \over {{P_{{\rm{Liver}}}}}},\end{align}
$${{{\rm{d}}{A_{{\rm{Blood}}}}} \over {{\rm{d}}t}} = \sum_{i \epsilon T}\left({{{Q_{i}{{\it{A_{i}}}}}} \over {{V_{i}{{\it{P_{i}}}}}}} - {Q_{i}A_{Blood}\over{V_{Blood}}}\right) + {F_{Q}{{Q_{{\it{Fat}}}}{A_{{\rm{Fat}}}}} \over {{V_{{\rm{Fat}}}}{P_{{\rm{Fat}}}}}} - F_{Q}Q_{Fat}{{{A_{{\rm{Blood}}}}} \over {{V_{{\rm{Blood}}}}}} $$
with T = {Kidney, Muscle, Rest, Udder, Liver},
$${{d{A_{{\rm{Fat}}}}} \over {{\rm{d}}t}} = {F_Q}{Q_{{\rm{Fat}}}}{{{A_{{\rm{Blood}}}}} \over {{V_{{\rm{Blood}}}}}} - {{{F_Q}{Q_{{\rm{Fat}}}}{A_{{\rm{Fat}}}}} \over {{V_{{\rm{Fat}}}}{P_{{\rm{Fat}}}}}},$$
$${{{\rm{d}}{A_{{\rm{Kidney}}}}} \over {{\rm{d}}t}} = {Q_{{\rm{Kidney}}}}{{{A_{{\rm{Blood}}}}} \over {{V_{{\rm{Blood}}}}}} - {{{Q_{{\rm{Kidney}}}}{A_{{\rm{Kidney}}}}} \over {{V_{{\rm{Kidney}}}}{P_{{\rm{Kidney}}}}}},$$
$${{{\rm{d}}{A_{{\rm{Muscle}}}}} \over {{\rm{d}}t}} = {Q_{{\rm{Muscle}}}}{{{A_{{\rm{Blood}}}}} \over {{V_{{\rm{Blood}}}}}} - {{{Q_{{\rm{Muscle}}}}{A_{{\rm{Muscle}}}}} \over {{V_{{\rm{Muscle}}}}{P_{{\rm{Muscle}}}}}},$$
$${{{\rm{d}}{A_{{\rm{Udder}}}}} \over {{\rm{d}}t}} = {Q_{{\rm{Udder}}}}{{{A_{{\rm{Blood}}}}} \over {{V_{{\rm{Blood}}}}}} - {{{Q_{{\rm{Udder}}}}{A_{{\rm{Udder}}}}} \over {{V_{{\rm{Udder}}}}{P_{{\rm{Udder}}}}}} - {{{\rm{C}}{{\rm{L}}_{{\rm{Milk}}}}{A_{{\rm{Udder}}}}} \over {{V_{{\rm{Udder}}}}}}.$$
Thus, the concentration in milk fat is again given by
$${C_{{\rm{Milk}}}} = {{{A_{{\rm{Udder}}}}} \over {{V_{{\rm{Udder}}}}}}.$$
Similar to the Derks model(Reference Derks, Berende and Olling27), the transition between each compartment depends on the blood flows
${Q_i}$
[L/d], the compartment volumes
${V_i}$
[L/d] (both of which depend on the properties of the individual cow) and the partition coefficient
${P_i}$
[unitless], which reflects the physico-chemical properties of the contaminant by describing the tissue–blood ratio of the contaminant in the stationary state. Additionally, the milk excretion model is the same as in Derks, that is, proportional to the amount of milk fat excreted
${\rm{C}}{{\rm{L}}_{{\rm{Milk}}}}\;$
[L/d]; likewise, the metabolism follows linear kinetics with rate
${k_{{\rm{met}}}} = {\rm{C}}{{\rm{L}}_{{\rm{Liver}}}}/{P_{{\rm{Liver}}}}\;$
[1/d], where
${\rm{C}}{{\rm{L}}_{{\rm{Liver}}}}$
[1/d] is the liver clearance.
The parameters should all be taken from the literature, except for the partition coefficient they proposed, which can be calculated using the contaminant’s log
$\left( {{K_{{\rm{ow}}}}} \right)$
value if no further information is available. But as mentioned in the classical PBTK approach by Derks, such a method suffers from prediction problems for PCDD/Fs and PCBs. An alternative would be to predict partition coefficients with other methods (see e.g. Graham et al. (2011)(Reference Graham, Walker and Jones31) or Endo et al. (2013)(Reference Endo, Brown and Goss32)).
An even more general model that considers multiple trophic levels for several kinds of contaminants was developed by Hendriks et al. (2001)(Reference Hendriks, van der Linde and Cornelissen39). It was later adapted to cattle by Hendriks et al. (2007) to calculate the BTF of various contaminants into milk and beef(Reference Hendriks, Smítková and Huijbregts40). For lactating cows, this latter model essentially boils down to a one-compartment model with multiple input and output sources (Fig. 9), yielding a differential equation of the form
\begin{align}{{{\rm{d}}{C_{{\rm{Cow}}}}\left( t \right)} \over {{\rm{d}}t}}{\rm{ }} &= {k_{{\rm{in}},n}}{C_{{\rm{Feed}}}} + {k_{{\rm{in}},w}}{C_{{\rm{Water}}}}\\
&\quad - \left( {{k_{{\rm{out}},n}} + {k_{{\rm{out}},w}} + {k_p} + {k_{{\rm{met}}}}} \right){C_{{\rm{Cow}}}}\left( t \right).\end{align}
Fig. 9. Schematic description of the multitrophic level model of Hendriks et al. (2001)(Reference Hendriks, van der Linde and Cornelissen39), adapted to the lactating cow(Reference Hendriks, Smítková and Huijbregts40). The source of contamination could be feed, divided into water and lipid, or just water. The absorption rate of both,
${k_{{\rm{in}},i}}$
[1/d], is derived assuming that these contaminants must first pass through both water and lipid layers to enter the cow. The excretion of contaminants is divided into urinal excretion represented as water in the model on the one hand, and biomass excretion on the other (e.g. milk), which is further divided into water and lipid. The excretion rates
${k_{{\rm{out}},i}}$
[1/d] from the system are influenced by a water and lipid layer, as was the case for absorption. In addition, the reduction of the contaminant concentration in the cow’s body can occur via metabolism or dilution of the biomass with the rate constants
${k_{{\rm{met}}}}$
[1/d] or
${k_p}$
[1/d].
and the concentration in milk fat is thus given by
$${C_{{\rm{Milk}}}}\left( t \right) = {C_{{\rm{Cow}}}}\left( t \right){{{V_{{\rm{Cow}}}}} \over {{V_{{\rm{Milk}}}}}}{k_{{\rm{Milk}}}}$$
Here,
${k_{{\rm{in}},n}}$
and
${k_{{\rm{out}},n}}$
[1/d] are the input and output rates via feed, where
${k_{{\rm{out}},n}}$
includes the excretion with milk fat at rate
${k_{{\rm{Milk}}}}\;\left[ {1/{\rm{d}}} \right]$
;
${k_{{\rm{in}},w}}$
and
${k_{{\rm{out}},w}}$
[1/d] are the input and output rates via water (irrelevant for highly hydrophobic contaminants such as PCDD/Fs and PCBs). Additionally, elimination of the substance can happen via metabolism/transformation with rate constant
${k_{{\rm{met}}}}$
, and dilution of biomass (e.g. growth) with rate constant
${k_p}$
. The concentration in food and water are given by
${C_{{\rm{Feed}}}}$
[ng/L] and
${C_{{\rm{Water}}}}$
[ng/L]. Finally, the volumes of the cow and its daily milk fat yield is given by
${V_{{\rm{Cow}}}}\;\left[ {\rm{L}} \right]$
and
${V_{{\rm{Milk}}}}\;\left[ {{\rm{L}}/{\rm{d}}} \right]$
, respectively.
One of the main focuses of Hendriks (2001) was to show how to calculate the rate constants, especially
${k_{{\rm{in}}}}$
and
${k_{{\rm{out}}}}$
(Reference Hendriks, van der Linde and Cornelissen39). For these, it was assumed that the contaminant moves in a path through lipid and water layers upon both entering and leaving the animal via feed or water, similarly to the approach by McLachlan (1994)(Reference McLachlan38). From this, they derived formulas describing
${k_{{\rm{in}}}}$
and
${k_{{\rm{out}}}}$
only depending on the
${K_{{\rm{ow}}}}$
value of the contaminant and the weight of the animal. For the dilution of biomass constant
${k_p}$
, they assume it also scales with the weight of the animal. Lastly, for the elimination via metabolism, the model has to be fitted using experimental data.
As an aside, we note that models related to Hendriks’ have been developed for broader applications. For example, the model for transfer from feed into cow’s milk is only one part of a larger model for PCDD/Fs and PCBs along the human food chain (e.g. ACC-Human)(Reference Czub and McLachlan41).
Calculating transfer parameters from toxicokinetic models
The compartment models described in this review can be used to calculate transfer parameters, such as congener-specific elimination half-lives and transfer rates mentioned in part I of the review chapter on Kinetic parameters to characterise the feed-to-milk transfer behaviour(Reference Krause, Moenning and Lamp6). While we always recommend using a full model in risk analysis instead of transfer parameters, calculating them allows for easy comparison among congeners, among mathematically diverse models and against experimental data; it also provides measures of transfer that are more intuitive to communicate. To calculate transfer parameters, we assume that the model parameters are constant over time (i.e. compartment values, input vector, etc.). To illustrate the present discussion, we can rewrite all these models in standard linear algebraic notation, that is,
$${{{\rm{d}}A\left( t \right)} \over {{\rm{d}}t}} = {\rm{MA}}\left( t \right) + I,$$
where
$A\left( t \right)$
is the time-dependent amount vector containing the amount of contaminant in each compartment at time
$t$
;
$M$
is the transition matrix containing in its elements the transition rates between the compartments and I is the input vector containing the amount added into each compartment from outside, that is, feed; these are the model parameters are assumed to be independent of time. For a more detailed description, see Supplementary Material Chapters 1–8.
Calculating TR, TF and BTF for multicompartment models
Given a multicompartment model with a constant invertible transfer matrix
$M$
and input vector
$I$
(Reference Klevenhusen, Sudekum and Breves8), we first need to calculate the steady-state solution of this system. This is accomplished by inserting both into the formula
$${A_{{\rm{ss}}}} = {M^{ - 1}}I$$
or in the case of fugacity models
$${f_{{\rm{ss}}}} = {M^{ - 1}}I$$
Here
${M^{ - 1}}$
is the inverse of the transfer matrix
$M$
, which can be calculated with numerical methods. Then
${A_{{\rm{ss}}}}$
is the amount vector in steady state, that is, the quantity of contaminant in each compartment, and
${f_{{\rm{ss}}}}$
is the fugacity vector in steady state, respectively. In the case of the one compartment model by Hendriks et al. (2001)(Reference Hendriks, van der Linde and Cornelissen39), the steady state
${C_{{\rm{Cow}},{\rm{ss}}}}\;$
concentration can be directly calculated as
$${C_{{\rm{Cow}},{\it{ss}}}} = {{{k_{{\rm{in}},n}}{C_{{\rm{Feed}}}}} \over {{k_{{\rm{out}},w}} + {k_{{\rm{out}},n}} + {k_p} + {k_{{\rm{met}}}}}}.$$
Here we assume that there is no input via water into the system (
${k_{{\rm{in}},w}} = 0$
), as we consider only the transfer from feed. The transfer parameters discussed in the chapter on kinetic parameters to characterise the feed-to-milk transfer behaviour from part I of this review(Reference Krause, Moenning and Lamp6) can now be calculated for each compartment model type presented here using the formulas in Table 2.
Table 2. Formulas for calculating the transfer parameters discussed in part I of the review chapter Kinetic parameters to characterise the feed-to-milk transfer behaviour. Here
${A_{ \cdot ,{\rm{ss}}}}$
[ng] and
${f_{ \cdot ,{\rm{ss}}}}$
[Pa] are the steady-state amounts and fugacities respectively in the respective compartment for each model. Additional
${V_ \cdot }$
[L] is the volume of the respective compartment;
${\rm{C}}{{\rm{L}}_{{\rm{Milk\;}}}}$
[ng/d] is the amount of milk fat excreted each day;
${\rm{Dose}}\;$
[ng/d] or [mol/d] in the fugacity context is the amount of contaminant given to the animal each day;
${\rm{\;Feed}}$
[kg/d] is the amount feed given to cow each day;
${C_{{\rm{Milkfat}}}}$
[unitless] is the milk fat concentration;
${P_ \cdot }$
[unitless] is the partition coefficient for respective compartment and blood; finally
${D_{{\rm{Milk}}}}$
[mol/(Pa·d)] is the milk transport coefficient of the fugacity models;
${k_ \cdot }$
[1/d] are the respective transition rates in Hendriks’ model(Reference Hendriks, Smítková and Huijbregts40)
Calculating the elimination half-lives for multicompartment models
For a given n-compartment model, the half-lives can be also calculated from the n eigenvalues
${\lambda _i}$
of the transition matrix M. For this, we can use numerical algorithms, as a symbolic evaluation becomes involved for transition matrices of models with more than two compartments. Knowing the eigenvalues, the half-lives are
$${\tau _i} = {{{\rm{ln}}\left( 2 \right)} \over { - {\lambda _i}}}{\rm{with \ {\it{i}} \ {\it{in}} \{ 1,}} \ldots {\rm{,n\} }}{\rm{.}}$$
As already mentioned, there are usually more than two half-lives, but most of them are either too short to be relevant for risk assessment or are almost identical to each other. This effectively leaves us with only two of the
${\tau _i}$
’s being truly different practical observable half-lives: the shorter one (the
${\rm{\alpha }}$
half-life) at the start of the depuration and the longer one (
${\rm{\beta }}$
) at the end.
Conclusions
In this review, we examined a wide range of toxicokinetic models developed to predict the transfer of PCDD/Fs and PCBs from feed to milk. These models vary in complexity, ranging from black-box approaches to others that closely mimic cow physiology and fugacity models based on thermodynamic equations. An overview of the strengths and limitations of each approach is summarised in Table 1. Because transfer parameters such as TR, TF, BTF, and half-lives are important to understand and compare models and congeners, we have also provided a guide for extracting these parameters from each toxicokinetic model discussed in this review.
What is the ideal model for risk assessors to use for predicting PCDD/F and PCB transfer into milk as a consequence of oral exposure? An ideal model has been validated with multiple datasets(Reference Savvateeva, Ohlhoff, Hoogenboom, Pieper and Numata42) and can predict the complete congener-specific spectrum of substances in question. Furthermore, it should include proper physiological modeling to allow extrapolation according to a specific cow (herd) metabolic and health status, such as body weight, body fat, milk yield and milk fat yield. Unfortunately, we have to report that no model currently satisfies all these criteria simultaneously.
For the fugacity approach, non-steady-state validation has only been performed in the work of McLachlan (1992) for PCB-138, but only the elimination phase used for calibration could be accurately described(Reference McLachlan28). The newer versions of the fugacity approach were only evaluated at a near-steady state(Reference Tremolada, Guazzoni and Parolini13,Reference Rosenbaum, McKone and Jolliet36) . While we currently cannot recommend these fugacity approaches for dynamic prediction of content in milk fat owing to the lack of validation, the approach can be used alternatively to approaches presented in the chapter on non-physiological approaches for calculating transfer parameters to predict the TR, TF or BTF as already shown by Rosenbaum (2009)(Reference Rosenbaum, McKone and Jolliet36).
The classic PBTK approach of Derks was applied and calibrated to data published by Derks et al. (1994) for TCDD(Reference Derks, Berende and Olling27) and by Hoogenboom et al. (2010) for a mixture (PCDD/F WHO2005 TEQ)(Reference Hoogenboom, Zeilmaker and van Eijkeren2), with both parametrisations showing good performance against their respective datasets. For this reason, we would currently recommend the use of these models for TCDD and PCDD/F WHO2005 TEQ, respectively, although they do not fulfil all criteria mentioned above. An implementation of the Hoogenboom et al. (2010) model(Reference Hoogenboom, Zeilmaker and van Eijkeren2) can be found in the RIVM/WFSR tool www.FeedFoodTransfer.nl. For other congeners, there are only theoretically parametrised approaches that have not yet been sufficiently validated(Reference Tremolada, Guazzoni and Parolini13,Reference van Eijkeren, Jager and Sips29,Reference MacLachlan30,Reference Bogdal, Züst and Schmid34,Reference Rosenbaum, McKone and Jolliet36) . We recommend caution when employing them and encourage the community to perform additional validation. It would be beneficial if the models presented here were to be further validated for all congeners using independent datasets to assess predictive accuracy. This is also true for the Derks (1994) and Hoogenboom (2010) models since the validation dataset was also used for calibration. In addition, it would be interesting to see how well these models can predict changes in the excretion of these congeners caused by differences in cow (herd) metabolic and health status. The question of upscaling simulations to reflect whole herds is also not trivial, which is never directly addressed. It was only indirectly addressed in Hendriks (2007) by taking dilution biomass as a parameter into the model(Reference Hendriks, Smítková and Huijbregts40). This was also done in models, which deal with a much broader context, that are not discussed here such as ACC-Human(Reference Czub and McLachlan41).
Future model developers are well advised to follow the guidelines from Lautz et al.(Reference Lautz, Oldenkamp and Dorne43,Reference Lautz44) , which include basing them on generic and flexible model structures and incorporating tools to assess model performance. We encourage the community of modellers to pursue congener-specific, physiologically based models that can be extrapolated, used for herds, and have been developed and validated with a multiplicity of independent datasets.
Acknowledgements
We would like to thank Dr. Matthew Salewski for the technical and linguistic proofreading of this manuscript.
H.S. and P.F. acknowledge the support from the Maria Sibylla Merian Fellowship; the Fellowship had no role in the design, analysis or writing of this article. All authors received no specific grant from any funding agency, commercial or not-for-profit sectors.
There are no conflicts of interest.
Jan-Louis Moenning: literature search, modelling aspects, writing (original draft), figures
Torsten Krause: literature search, data compilation and analysis, writing (review and editing), figures
Julika Lamp: literature search, data gathering, introduction, animal health and physiology aspects, figures
Ronald Maul: conceptualisation, writing (review and editing)
Hans Schenkel: conceptualisation, literature search
Peter Fürst: conceptualisation, literature search
Robert Pieper: animal health and physiology aspects, conceptualisation, writing (review and editing)
Jorge Numata: supervision, conceptualisation, modelling aspects, writing (original draft, review and editing), figures
Supplementary material
To view supplementary material for this article, please visit https://doi.org/10.1017/S0954422422000208.
Open access
