Introduction
In this paper, we consider a generic one-compartment toxicokinetic (TK) model describing uptake and elimination processes within living organisms when exposed to chemical substances. From a general point-of-view, the term TK refers to the movement and fate (also called disposition) of chemical substances that are potentially toxic. In particular, this term is used when describing the time course of absorption, distribution and elimination (accounting for both biotransformation and excretion) of substances within organisms (ADME). Toxicokinetics is also related to pharmacokinetics and is often considered as the same scientific field, the main difference being linked to the type of compound (toxicants or pharmaceuticals). In comparison with pharmaceutical exposure, exposure to toxicants is often uncontrolled and variable over time predominantly at larger concentrations.1
Mathematical TK models are particularly useful to predict the adverse effects of xenobiotics or to prevent undesired residues within animal tissues from entering the human food chain.2 Indeed, TK models allow the prediction of tissue concentrations over time, involving several kinetic parameters, such as uptake and elimination rates, from which useful metrics for assessing risk may be derived.3 Additionally, TK models can be used to facilitate a better understanding of the underlying physiological mechanisms driving chemical ADME processes.4,5
A wide variety of TK models exists, chosen from both the available data and the use intended for the model (e.g., describe and/or predict); simple or more complex models may be preferred accordingly.6 Basically, we distinguish classic TK models from physiologically-based (PB) TK models, referred to as PBTK, or PBPK (standing for PB Pharmaco-Kinetic) models in the field of pharmacology. As fully detailed in Gehring et al.,7 PBPK models are a crucial tool in Environmental Risk Assessment (ERA) for regulatory bodies, such as the US Environmental Protection Agency (US EPA) and European Food Safety Authority (EFSA).8,9 In fact, PBPK models are an ethical and scientifically sound method to predict exposure to toxic xenobiotics in humans through animal-to-human extrapolation or based on human biomonitoring data. (PB)TK models are also of particular interest for ERA, when the calculation of bioaccumulation metrics is required by regulators to evaluate dossiers for substance marketing authorizations.10–12
The generic TK model we unravel in this paper is a one-compartment model that considers organisms as a whole, in which chemical compounds may enter and from which these compounds can be eliminated. More complex TK models (namely PBTK models) refine the description of contamination pathways within organisms, distinguishing organs and tissues from physiological hypotheses on potential targets of exposure compounds; PBTK models have been mainly developed for humans13,14 and big mammals.15,16 Recently, a PBTK model was proposed by Gestin et al.17 for invertebrates, pushing the immense potential of these models a step further.
Our TK model, even if not fully refined as a PBTK model, has the great advantage of accounting for several exposure sources (e.g., water, sediment, and/or food), several elimination processes (e.g., direct elimination, dilution by growth and/or biotransformation), and several potential metabolites of the parent chemical compound to which organisms are exposed.18 The next section details the dynamic system of ordinary differential equations (ODEs) and the mathematical exact solution, followed by a summary of the full set of ODE in section 3. Then, section 4 illustrates the use of model based on simulations from parameter values obtained by the study of real datasets.
One-compartment TK set of ODE
Our generic one-compartment TK model is composed of two sets of ODE: accumulation and depuration sets. The accumulation phase includes both absorption and elimination processes, as those related to the biotransformation of the parent compound into metabolites, during which organisms are exposed to a given compound. The depuration phase concerns only elimination processes, including biotransformation, during which organisms are transferred into a clean medium. The transition from one set of ODE to the other takes place at time tc corresponding to the duration of the accumulation phase (see Table 1).
Table 1Symbols, meaning, and units of all parameters and variables involved in the full set of ordinary differential equations defining the generic one-compartment toxicokinetic model
| Symbol | Meaning and unit |
|---|
| t | time (in time units) |
| tc | duration of the accumulation phase (in time units) |
| Cp(t) | internal concentration of the parent compound at t (in μg·g−1) |
| Cmℓ(t) | internal concentration of metabolite ℓ at t (in μg·g−1) |
| U=∑i=1Ikuici | sum of all uptake terms |
| ci | concentration via the exposure route i (in μg·mL−1) |
| E=∑j=1Jkej | sum of all elimination terms for the parent compound |
| M=∑ℓ=1Lkmℓ | sum of all elimination terms for metabolite ℓ |
| i | index of exposure sources, i = 1…I |
| j | index of elimination processes, j = 1…J |
| ℓ | index of metabolites, ℓ = 1…L |
| I | total number of exposure sources (#) |
| J | total number of elimination processes (#) |
| L | total number of metabolites (#) |
| kui | uptake rate of exposure source i (per time units) |
| kej | excretion rates of elimination process j (per time units) |
| keℓ | excretion rates of metabolite ℓ (per time units) |
| kmℓ | metabolization rate of metabolite ℓ (per time units) |
| R=UE+M | NA |
| Dℓ=keℓ−(E+M) | NA, ∀ℓ = 1…L |
| Q=C0−R(1−e(E+M)tc) | NA |
In practice, the exposure concentration to which organisms are exposed may vary over time in real environments. This means there is no analytical solution of the TK model, so that only a numerical solution can be obtained with an appropriate algorithm. Therefore, our work assumes that the exposure concentration remains constant over time regardless of the exposure source. Such an experimental condition can be ensured for most chemical compounds when performing laboratory experiments. In addition, this assumption provides the exact solution of the full one-compartment TK model by considering as many routes of exposure and as many elimination processes as desired, as well as an infinite number of phase I metabolites, i.e., directly derived from the parent compound to which organisms are exposed. Kuo et al.19 previously suggested a partially resolved one-compartment TK model but only for the accumulation phase and one exposure route. In a similar manner, D’Argenio et al.12 published an analytical solution for a linear multi-compartments TK model with a non-zero initial condition.
Accumulation phase (0 ≤ t ≤ tc)
The set of ODEs describing the accumulation phase (0 ≤ t ≤ tc) is as follows:
{dCp(t)dt=U−(E+M)×Cp(t) (a)dCmℓ(t)dt =kmℓ× Cp(t)−keℓ×Cmℓ(t), ∀ℓ=1 ⋯L (b)
All parameters and variables, with their meaning and units when applicable, are provided together in Table 1.
Equation (1.a) for the parent compound is a linear first-order ODE with constant coefficients and a second member. Equation (1.a) admits Cpart(t)=UE+M=R as a particular solution. Equation (1.a) without its second member is:
dCp(t)dt−(E+M)×Cp(t)=0⇔Cp(t)=K×e−(E+M)×t, K∈ℝ+
Given the initial condition C (t = 0) = C0 (C0 ≥ 0), we can obtain the full analytical solution of Equation (1.a), providing that the internal concentration of the parent compound over time during the accumulation phase (0 ≤ t ≤ tc) is as follows:
Cp(t)=(C0−R)×e−(E+M)×t+R
See Table 1 for the definition of parameter R.
Equation (1.b) is also a linear first-order ODE with constant coefficients and a second member, with the following analytical solution when the second member is removed:
dCmℓ(t)dt−keℓ×Cmℓ(t)=0⇔Cmℓ(t)=K×e−keℓ×t, K∈ℝ+, ∀ℓ=1…L
The method of variation of a constant consists of writing the general solution of Equation (1.b) as:
Cmℓ(t)=K(t)×e−keℓ×t
which is used to find function K(t) by deriving Cmℓ(t) then insert the result into Equation (5). The derivative of Equation (1.b) then becomes: dCmℓ(t)dt=dK(t)dt×e−keℓ×t−K(t)×keℓ×e−keℓ×t
Inserting the result from the previous equation into Equation (1.b) leads to:
dK(t)dt=kmℓ×((C0−R)×eDℓ×t+R×ekeℓ×t)
which integrates into: K(t)=kmℓ×(C0−RDℓ×eDℓ×t+Rkeℓ×ekeℓ×t)+C, C∈ℝ
See Table 1 for the definition of parameter Dℓ.
The general solution of Equation (1.b) finally becomes:
Cmℓ(t)=kmℓ×(C0−RDℓ×e−(E+M)×t+Rkeℓ)+C×e−keℓ×t, C∈ℝ
When considering the initial condition Cmℓ(t=0)=0 when the accumulation phase starts, organisms are only exposed to the parent compound, so that there are no metabolites. With this condition, we can finally acquire the full analytical solution of Equation (1.b), providing the internal concentration of metabolite ℓ over time during the accumulation phase:
Cmℓ(t)=kmℓ×(C0−RDℓ×(e−(E+M)×t−e−keℓ×t)+Rkeℓ×(1−e−keℓ×t))
Depuration phase (t > tc)
The set of ODEs describing the depuration phase (t > tc) is as follows:
{dCp(t)dt=−(E+M)×Cp(t) (a)dCmℓ(t)dt =kmℓ× Cp(t)−keℓ×Cmℓ(t), ∀ℓ=1 ⋯L (b)
The definitions and units (when applicable) of the parameters and variables are described together in Table 1.
Equation (11.a) is a linear first-order ODE without a second member, so that it has a general solution of the following form:
Cp(t)=K×e−(E+M)×t, K∈ℝ+
For the depuration phase and the parent compound, the initial condition comes from the calculation of the internal parent compound concentration at the end of the accumulation phase (i.e., at time t = tc) using Equation (13):
Cp(tc)=(C0−R)×e−(E+M)×tc+R
From the general analytical solution given by Equation (12), we get Cp(tc)=Ke−(E+M)tc, from which constant K in Equation (12) is determined by:
K=C0−R+R×e(E+M)×tc
Then, the final analytical solution of Equation (11.a), providing the internal concentration of the parent compound over time during the depuration phase (t > tc), becomes:
Cp(t)=(C0−R×(1−e(E+M)×tc))×e−(E+M)×t
For simplicity reasons, Equation (15) can be written as Cp(t) = Qe−(E+M)×t, where Q is defined in Table 1.
Equation (11.b) is a linear first-order ODE with constant coefficients and a second member. The analytical solution of Equation (11.b) without its second member is:
Cmℓ(t)=K×e−keℓ×t, K∈ℝ+
As previously stated, the method of the variation of a constant provides the general solution of Equation (11.b) as Cmℓ(t)=K(t)e−keℓt, from which function K(t) must be determined.
The derivative of Cmℓ(t)Cmℓ(t) from Equation (16) is:
dK(t)dt=e−keℓ×t−keℓ×K(t)×e−keℓ×t
Inserting Equation (17) into Equation (11.b) leads to:
dK(t)dt=kmℓ×Q×eDℓ×t
which integrates into: K(t)=kmℓ×QDℓ×eDℓ×t+C, C∈ℝ
finally leading to the general analytical solution of Equation (11.b): Cmℓ(t)=kmℓ×QDℓ×e−(E+M)×t+C×e−keℓ×t, C∈ℝ
Constant C can be determined from the initial condition, i.e., from the internal concentration of metabolite ℓ at t = tc both at the end of the accumulation phase and at the beginning of the depuration phase. From Equation (20), we get Cmℓ(tc)=kmℓQDℓe−(E+M)tc+Ce−keℓtc, and from Equation (10), we obtain Cmℓ(tc)=kmℓ(C0−RDℓ(e−(E+M)tc−e−keℓtc)+Rkeℓ(1−e−keℓtc)). Finally, the expression for constant C can be determined using:
C=kmℓ×(Rkeℓ×(ekeℓ×tc−1)−C0−RDℓ−RDℓ×ekeℓ×tc)
Replacing constant C in Equation (20) gives the final analytical solution of Equation (11.b), which yields the internal concentration of metabolite ℓ for the depuration phase (t > tc) as follows:
Cmℓ(t)=kmℓ×(QDℓ×e−(E+M)×t+Rkeℓ×(e−keℓ×(t−tc)−e−keℓ×t) −C0−RDℓ×e−keℓ×t−RDℓ×e−keℓ×(t−tc))
Replacing constant Q by its own expression (Table 1) leads to:
Cmℓ(t)=kmℓ×(C0−RDℓ×(e−(E+M)×t−e−keℓ×t) +Rkeℓ×(e−keℓ×(t−tc)−e−keℓ×t) +RDℓ×(e−(E+M)×(t−tc)−e−keℓ×(t−tc))
Model simulations based on real case-studies
Currently, the above exact mathematical solution of the generic one-compartment TK model is used in support of the development of the R-package rbioacc20 as well as its web interface MOSAICbioacc, which is freely available at https://mosaic.univ-lyon1.fr/bioacc .21 Both tools allow fitting the TK model to bioaccumulation data in a user-friendly way, providing all necessary outputs to check for goodness-of-fit and obtain bioaccumulation metrics, regardless of the considered species-compound combination. The database associated to the web interface is also available online at http://lbbe-shiny.univ-lyon1.fr/mosaic-bioacc/data/database/TK_database.html , with a wide collection of accumulation-depuration data in support of toxicokinetic modelling.22 From this database, we selected three real case-studies from which we got fitting results (especially the joint posterior probability distribution of all kinetic parameters) in order to perform the subsequent model simulations, thus illustrating the use of the generic TK model.
All model simulations of the generic set of mathematical solutions (Equations (3), (10), (15) and (23)) were performed using R software with the rbioacc package and function predict().20 Note that the same simulations can be performed directly online with the MOSAICbioacc platform and its prediction module (https://mosaic.univ-lyon1.fr/bioacc ). We used 500 time points in [0; tf], where tf stands for the final time of each simulation. These simulations illustrate the use of the generic TK model in three case studies where different compounds are bioaccumulated by different species. In order to proceed, the time duration of the accumulation phase (parameter tc) is required, as well as the exposure concentrations in the media and the model parameter values (Table 2); all associated numerical values were taken from fitting results performed with real data sets. Figures 1 to 5 show the simulations of internal concentrations over time for the three species-compound combinations, as subsequently described. For each case study, TK model parameters were varied one-at-a-time with 20, 50, and 80% of increase from their original value.
Table 2Sets of inputs for model simulations according to Figures 1 to 5
| Input | Units | Simple (section 4.1) | Multiple exposures (section 4.2) | Biotransformation (section 4.3) |
|---|
| cw | μg·mL−1 | 0.0044 | – | 15.53 |
| cs | μg·g−1 | – | 56.6 | – |
| cf | μg·g−1 | – | 1.46 | – |
| tc | days | 49 | 7 | 1 |
| kuw | days−1 | 10.46 | – | 16,740 |
| kus | days−1 | – | 0.071 | – |
| kuf | days−1 | – | 0.013 | – |
| km1 | days−1 | – | – | 73.27 |
| km2 | days−1 | – | – | 0.5166 |
| km3 | days−1 | – | – | 0.1957 |
| kee | days−1 | 0.04 | 0.178 | 4.164 |
| ke1 | days−1 | – | – | 561 |
| ke2 | days−1 | – | – | 0.123 |
| ke3 | days−1 | – | – | 0.7808 |
Elementary one-compartment TK model
Equations (3) and (15) were simulated for a simple case study in which fish are exposed to a highly hydrophobic chemical contaminated water and only excretion is considered (Fig. 1 and Table 2).5 The corresponding inputs are the exposure concentration cw (referring to ci and I = 1), uptake rate from water kuw (referring to kui and I = 1), and excretion rate kee (referring to kej and J = 1). When parameter kuw increases, internal concentrations become higher (e.g., blue curve in Fig. 1a) than with the original value (orange curve in Fig. 1a). Biologically speaking, the higher the uptake rate is for a given substance, the more it is bioaccumulated by organisms. Conversely, an increase in parameter kee leads to a lower internal concentration (e.g., blue curve in Fig. 1b), which is consistent with the known underlying biological mechanism: the faster a contaminant is eliminated, the quicker its concentration decreases in organisms.
One-compartment TK model with several exposure routes, no metabolites
Equations (3) and (15) are simulated for case, in which freshwater shrimp are exposed to an organic chlorine compound through contaminated sediment and food, and only natural excretion is considered as the elimination process (Fig. 2 and Table 2).23 The corresponding inputs are the exposure concentration via sediment cs and food cf (referring to ci and I = 2), the two uptake rates from sediment and food, kus and kuf (referring to kui and I = 2), and excretion rate kee (referring to kej and J = 1). When parameter kus increases, internal concentrations are higher (e.g., blue curve in Fig. 2a) than with the original value (orange curve in Fig. 2a). Biologically, the higher the uptake rate from sediment is for a given substance, the stronger it is bioaccumulated by organisms. Regarding the exposure from food, as parameter kuf is low in this case, its variation has minimal influence on the internal concentration of the contaminant (Fig. 2b). This means that the exposure via sediment is the major route of contamination for these organisms. Besides, as previously shown in section 4.1, an increase in parameter kee leads to a faster decreasing concentration (e.g., blue curve in Fig. 2c).
One-compartment TK model with one exposure route, several metabolites
Equations (3), (10), (15) and (23) are simulated for case, in which freshwater shrimp are exposed to an organic biocide by contaminated water. Three metabolites derived from the parent compound are considered24 together with the natural excretion (Figs. 3 to 5 and Table 2). The corresponding inputs are the exposure concentration via water cw (referring to ci and I = 1), the uptake rate kuw (referring to kui and I = 1), excretion rate kee (referring to kej and J = 1), three metabolization rates km1, km2 and km3 (referring to kmℓ and L = 3), and three elimination rates of the metabolites ke1, ke2 and ke3 (referring to to keℓ and L = 3). When parameter km1 increases, internal concentrations of metabolite 1 increase (e.g., blue curve in Fig. 3b) greater than with the original value (orange curve in Fig. 3b). Conversely, internal concentrations are lower for the parent compound (e.g., blue curve in Fig. 3a). In other words, the more the biotransformation rate for a given metabolite increases, the higher its concentration becomes within organisms due to the highly biotransformation of the parent compound. This leads to a lower internal concentration of the parent compound than with the original value of km1 (e.g., blue curve, Fig. 3a). An increase in km1 also induces a decrease in the internal concentrations of the other metabolites (Figs. 3c and d). Besides, an increase in parameter ke1 only affects the internal concentration of metabolite 1 (Fig. 4b). In addition, as previously viewed (sections 4.1 and 4.2), an increase in parameter kuw will induce high internal concentrations for both the parent compound and its metabolites (Fig. 5). Indeed, the more the internal concentration of the parent compound increases, the more the biotransformation process will intensify, leading to high internal concentrations for each metabolite.