Generic Solving of Physiologically-based Kinetic Models in Support of Next Generation Risk Assessment Due to Chemicals

doi:10.14218/JERP.2022.00043

Original Article

OPEN ACCESS

Download PDF

Generic Solving of Physiologically-based Kinetic Models in Support of Next Generation Risk Assessment Due to Chemicals

Sandrine Charles^1,* ,
Ophelia Gestin^1,2,3,
Jérémie Bruset¹,
Dominique Lamonica¹,
Virgile Baudrot⁴,
Arnaud Chaumot²,
Olivier Geffard²,
Thomas Lacoue-Labarthe³ and
Christelle Lopes¹

Journal of Exploratory Research in Pharmacology 2023;8(2):140-154

doi: 10.14218/JERP.2022.00043

Received: May 16, 2022

Revised: June 19, 2022

Accepted: July 6, 2022

Published online: September 29, 2022

Author information

Citation: Charles S, Gestin O, Bruset J, Lamonica D, Baudrot V, Chaumot A, et al. Generic Solving of Physiologically-based Kinetic Models in Support of Next Generation Risk Assessment Due to Chemicals. J Explor Res Pharmacol. 2023;8(2):140-154. doi: 10.14218/JERP.2022.00043.

Abstract

Background and objectives

Increasing confidence in using in vitro and in silico model-based data to aid the chemical risk assessment process is one of the most significant challenges faced by regulatory authorities. A crucial concern is taking full advantage of scientifically valid physiologically-based kinetic (PBK) models. The present study aims to present a very innovative solution of a fully generic PBK model written as a set of ordinary differential equations (ODE).

Methods

This study proposes an innovative and unified modeling framework for writing PBK equations as matrix ODE and their solutions, expressed with matrix products. This generic PBK solution considers as many state variables as needed to quantify chemical absorption, distribution, metabolism, and excretion processes within living organisms when exposed to chemical substances.

Results

We first introduce our PBK modeling framework, with all the intermediate steps from the matrix ODE to the exact solution. Then we apply this framework to bioaccumulation testing before illustrating its concrete use through complementary case studies regarding species, compounds, and model complexity.

Conclusions

This generic PBK model makes it possible for any compartmentalization to be considered, as well as all appropriate interconnections between compartments and with the external medium.

Keywords

Pharmacokinetics, Toxicokinetics, General unified modelling framework, Bioaccumulation, Chemical exposure

Introduction

Physiologically-based kinetic (PBK) models encompass both physiologically-based pharmaco-kinetic and physiologically-based (eco-)toxico-kinetic models depending on the context.^1,2 Hence, PBK models can refer either to therapeutic drug development³ or environmental risk assessment.^4–6 All PBK models are compartment models employing ordinary differential equations (ODE) to quantify chemical absorption, distribution, metabolism (i.e., biotransformation), and excretion processes within living organisms when exposed to chemical substances.^7,8 One fundamental aspect of PBK model complexity is the degree of compartmentalization (i.e., differentiation of an organism into various tissues or organs).^9,10 This complexity translates into a high number of parameters usually valued from literature information or expert knowledge. Consequently, most PBK models have purely predictive usages, for example, for human health assessments where novel experiments are very limited, even impossible.¹¹ Nevertheless, recent advances in the development of Bayesian inference tools make it possible to deal with complex models^12–15 even with sparse data sets, which opens up new possibilities for PBK models to obtain parameter estimates associated with their uncertainties and to propagate this information to model-informed predictions. Also, as part of new approach methodologies , PBK models as well as in vitro-in vivo extrapolation approaches can be combined with bioactivity data, all of this helping prioritize thousands of “data poor” chemicals in human health risk assessment.^10,16

Relying a priori on the anatomical and physiological structure of the body, PBK model compartments usually correspond to target organs or tissues, possibly all interconnected and related to the external exposure medium. However, many PBK models only relate organs or tissues to blood or lymph, with only one or two organs or tissues linked to the external medium.^17,18 Unfortunately, such simplifications of PBK models are often justified based on technical bases rather than physiological ones. All PBK models are written as coupled ODE whose parameters correspond to fluxes between compartments, for which information is partly available in the scientific literature. The general tendency to develop more complex models to include as many physiological processes as possible, and the current computational capacities, may give the impression that the most complex models will be the most efficient. Nevertheless, with increasing parameters to calibrate, such complex models may be avoided in favor of balancing complexity and simplicity. Current regulatory documents for assessing the bioaccumulation of chemicals in organisms usually only require simple one-compartment models,^19–23 even if it is now recognized that it is necessary to also consider internal concentrations within target organs in order to fully capture the chemical bioaccumulation behavior, the specific role of organs, and the dynamics of toxic effects.²⁴

Multi-compartment models sometimes reveal the necessary, for example, to finely decipher the internal contamination routes of specific chemical compounds causing damage to only specific organs.^25–27 Additionally, PBK models can be crucial to predicting organ-level concentration-time profiles in a situation where animal testing is now prohibited, using PBK model information from one chemical substance to inform the development or evaluation of a PBK model for a similar chemical substance.³ In the perspective to enlarge and facilitate the use of PBK models, namely to include more compartments, to better estimate parameter values from data, and to better support a fine deciphering of underlying contamination processes after chemical exposure, there is today a clear need for user-friendly tools. From an automatized implementation, fully transparent and reproducible, such tools should simplify the use of any PBK models, preventing users from investing in technicalities, whatever the required number of compartments to consider physiologically, whatever the number of connections to account for between compartments in pairs or between compartments and the external medium, and whatever the species-compound combination of interest. Such tools seem the only way to gradually achieve greater acceptability of complex PBK models, even in a regulatory context.^10,28

Capitalizing on recent publications on TK models,^15,22 we present in this paper a very innovative solution of a fully generic PBK model written as a set of ODE. Benefiting from an exact solution is a tremendous advantage in numerical implementation. Indeed, it avoids discretizing ODE as there is no more numerical approximation. The solution is directly used for the inference process and the subsequent simulations. The gain in calculation time is enormous (more than 100-fold), associated with fair use of computer resources. Moreover, the new modeling framework we propose makes it possible to account for an infinite number of compartments, with all possible connections between pairs of compartments and between compartments and the exposure medium, independently of the investigated species or chemical substance. Indeed, we found a particularly condensed way to write a linear ODE system, which is typical of PBK models, thus allowing us to fully and exactly solve the ODE system to write an exact generic solution. This exact solution is fascinating when estimating many parameters related to many state variables, for which experimental data may be sparse, with few replicates and high variability.

In the methods, we first detail our generic modeling framework, together with notations of parameters and variables, providing the generic solution at the end of section 2. Then, the generic modeling framework is applied to the particular context of bioaccumulating chemical substances within organisms. We detail how to write the ODE system for both accumulation and depuration phases of standard bioaccumulation tests and then how to get the final generic solution to simulate internal concentrations over time from parameter estimates. Results presents four different situations where simulations were useful to predict what happens within organs and/or tissues according to the species under consideration and the chemical substance to which it is exposed. These case studies were chosen from the literature to be diverse and complementary in terms of questions that a PBK model can help to investigate. We thus present two case studies for the species Gammarus fossarum exposed to cadmium, for which internal concentrations have been measured within four organs: one case study with one-compartment PBK models for each organ considered separately; the other case study with a four-compartment PBK model. These case studies illustrate the added value of considering one single four-compartment model rather than four one-compartment models for each organ. The third case study concerns the sea cucumber exposed to six different antibiotics. Furthermore, the fourth case study concerns the species Danio rerio exposed to arsenic, whose bioaccumulation process is described with a six organ-based compartment PBK model.

Materials and methods

Generic modeling framework

Biologists often expect an exhaustive description of the phenomenon they are studying. In the same way, mathematicians will want to use the most sophisticated methods they know. However, all models are inherently wrong; only some of them will prove helpful.²⁹ As a consequence, the modeler should position between these two points of view to be efficient. Such a position is known as the parsimony principle by which the simplest model that adequately explains the data should be used; it was proposed in the 14^th century by William of Ockam, an English Franciscan friar, scholastic philosopher, and theologian.³⁰ In its general form, the parsimony principle, also referred to as Occam’s Razor, states that the simplest of competing explanations is the most likely to be correct. In model fitting, the simplest model providing a good fit will be preferred over a more complex one. The compromise is thus between the good description of the observed data and simplicity.

In this spirit, the most generic PBK models will rely on simplifying hypotheses to decipher enough internal mechanisms under chemical exposure and remain mathematically reasonable to be easily manipulated. Below is a non-prioritized short list of the most current hypotheses: (i) the exposure concentration is assumed constant over time; (ii) there can be any number of compartments in direct relation with any number of tissues and organs that are needed to consider on a biological point of view; (iii) all compartments can be connected two-by-two to all the others; (iv) the exposure contaminant can enter within each compartment; and additionally, (v) all compartments can be theoretically connected to the external medium, the final choice to be based on biological expertise.

From these hypotheses, Figure 1 gives the general schematic representation of exchanges between the external medium and compartments and between compartments themselves. Table 1 gathers all variables and parameters involved in the generic writing of the PBK model and is used in Figure 1.

Generic scheme of a multi-compartments physiologically-based kinetic model connecting <italic>n</italic> compartments two-by-two and each of them to the external medium at exposure concentration <italic>c<sub>x</sub></italic> (<italic>i, j</italic> ∊ [1;<italic>n</italic>]).

Fig. 1 Generic scheme of a multi-compartments physiologically-based kinetic model connecting n compartments two-by-two and each of them to the external medium at exposure concentration c_x (i, j ∊ [1;n]).

Refer to Table 1 for names, meanings, and units of variables and parameters.

Table 1

Variable and parameter names, meanings, and units used within the generic PBK model all along this paper

Names	Meaning	Unit
t	time	[t]
n	total number of compartments	#
i, j	Compartment numbers	i, j ∊ [1;n]
c_x	exposure concentration in the external medium	mass per volume
c_i(t)	internal concentration in compartment i at time t	mass per weight
k_u,i	uptake rate from the external medium to compartment i	[t]⁻¹
k_e,i	elimination rate from compartment i to the external medium	[t]⁻¹
k_i,j	input rate from compartment j to compartment i	[t]⁻¹
k_j,i	output rate from compartment j to compartment j	[t]⁻¹

The symbol # means dimensionless; [t] stands for a time unit.

Mathematical equations of the PBK model

From Figure 1, we can derive the entire system of the ODE describing the dynamics of the multi-compartments model when organisms are exposed to an external constant concentration c_x:

dci(t)dt=ku,icx−ke,ici(t)+∑j≠iki,jCj(t) −∑j≠ikj,ici(t) ∀i,j∈[1;n](1)

Names, meanings, and units of variables and parameters are provided in Table 1.

This full system of ODE for n compartments all related by pairs (Equation 1) can equivalently be written in a matrix way as follows:

dC(t)dt=Ucx+EC(t)(2)

where vector C(t) gathers all internal concentrations in compartments i at time t, i ∊ [1;n]:

C(t)=(c1(t)c2(t)…cn(t))T(3)

Vector U contains all uptake rates from the external medium at exposure concentration c_x:

U=(ku,1ku,2…ku,n)T(4)

Matrix E gathers both input and output rates between compartments two-by-two, together with the elimination rates from each compartment i, i ∊ [1;n]:

E=[eij]i,j∈[1;n] with (eii=−ke,i−∑j≠ikj,ieij=ki,j(5)

Equation (2) is a matrix ODE system, which is linear with a second member. It can be solved in two steps, as detailed below.

Generic solving the PBK model

The first step is to solve the matrix system (2) without its second member. Then, the second step consists in finding the final general solution using the method of the variation of constant.

Solving the ODE system without the second member

Removing the second member from the matrix ODE system (2) leads to the following system to solve:

dCwosm(t)dt=ECwosm(t)(6)

With C_wosm(t) the desired solution of Equation (6) without a second member (abbreviated by index wosm). Using matrix exponential immediately provides the solution:

Cwosm(t)=etEΩ1(7)

With Ω₁ a vector integration constant (∊ℝⁿ), and the following definition for the matrix exponential:

etE=∑k=0∞1k!(tE)k(8)

Getting the generic solution for the ODE system

For the second step, including the second member, we used the method of the variation of the constant, starting from the assumption that the final general solution with the second member (abbreviated by index wsm) can be written as follows:

Cwosm=etEΩ1(t)(9)

with a vector function Ω₁(t) to be determined.

Given that the exposure concentration is assumed constant (equal to c_x) in this paper, deriving Equation (9) and replacing terms in Equation (2) leads to the following result:

dΩ1(t)dt=e−tEUcx⇔Ω1(t)=(∫0te−τEdτ)Ucx+Ω2(10)

The final generic solution of Equation (2) will thus write as:

Cwsm(t)=(∫0te(t−τ)Edτ)Ucx+etEΩ2(11)

With Ω₂ ∊ ℝⁿ a constant to be determined.

From an initial condition C_wsm(t = 0) = C₀, we finally get Ω₂ = C₀, which leads to the following final particular solution of Equation (2):

Cwsm(t)=(∫0te(t−τ)Edτ)Ucx+etEC0(12)

It remains to calculate the matrix integral to achieve the final solution of the matrix ODE system (2).

Final expression of the PBK solution

As detailed in the Supplementary File 1, and using the definition of a matrix exponential from Equation (8), the matrix integral in Equation (12) can be calculated with the following expression:

∫0te(t−τ)Edτ=(etE−I)E−1(13)

where matrix I is the identity matrix, i.e. the (n × n) square matrix with ones on the main diagonal and zeros elsewhere, and E⁻¹ the inverse matrix of E.

It can immediately be deduced that the final solution of Equation (12) simplifies as follows:

Cwsm(t)=etE(E−1Ucx+C0)−E−1Ucx(14)

In the following sections, we employed this generic expression to go beyond and build generic physiologically-based kinetic models that we applied to different case-studies in the field of toxicology.

Application to bioaccumulation testing

Accumulation and depuration phases

Bioaccumulation is defined as an increase in contaminant concentrations inside living organisms following uptake from the surrounding medium (living media, food, even workplace for humans). Bioaccumulation results from dynamical processes of uptake and elimination that can be modeled with the above ODE system (Equation (2)). The extent to which bioaccumulation occurs within a given species determines the subsequent toxic effects. Hence, a better knowledge of bioaccumulation enables us to assess the risk of exposure to chemicals and to evaluate our ability to control their use and emissions in the field.³¹ Bioaccumulation is thus the net result of all uptake and elimination processes by egestion, passive diffusion, metabolization, excretion, and maternal transfer. Concomitantly, the organism’s growth modulates the bioaccumulation by diluting chemical quantities in increasing body or organ mass (Fig. 2).

Absorption, distribution, metabolism, and excretion processes and their relationships with effects and responses within living organisms, leading to toxicity or efficacy depending on the chemical substance they are exposed to.

Fig. 2 Absorption, distribution, metabolism, and excretion processes and their relationships with effects and responses within living organisms, leading to toxicity or efficacy depending on the chemical substance they are exposed to.

Bioaccumulation tests are usually mid-to-long term laboratory experiments designed to identify all the potential uptake pathways, including food and waterborne exposure routes.³² Bioaccumulation tests commonly comprise an accumulation phase followed by a depuration phase.^33,33 During the accumulation phase, organisms are exposed to a chemical substance of interest. After a specific time (for t ∊ [0;t_c]), with t_c fixed by the experimental design, organisms are transferred to a clean medium for the depuration phase (for t > t_c). The concentration of the chemical substance (and of its potential metabolites) within organisms is then measured internally at regular time points during both phases. From an ERA perspective, such data can be used finely to estimate bioaccumulation metrics.²²

If the bioaccumulation within organisms is widely studied for humans and large animals, namely, fish, birds, and farm animals,^34,35 this is less the case for invertebrates.³⁶ However, it is equally essential to decipher internal processes at the target organ level in invertebrates. Indeed, profoundly unraveling internal routes of chemical substances between organs after they enter the body is of great interest to better understand mechanisms implied in the subsequent effects on fitness, a phenomenon known as organotropism.^37–39 Among invertebrate species of interest, crustacean amphipods are already recognized as particularly relevant as aquatic biomonitors of trace metals.^36,40–43

Generic modeling of bioaccumulation

Within this context, the generic ODE system (Equation (2)) may be fully applied to describe, simulate and predict what happens within organs (in terms of internal concentration over time) and between organs (in terms of uptake, elimination, and exchange rates) when an organism is exposed to a given chemical substance. To this end, each organ can be associated with one model compartment, leading to the following equations for both accumulation and depuration phases:

dCA(t)dt=Ucx+ECA(t)(15)

Equation (15) is identical to Equation (2), denoting C_A(t) the internal concentration at time t during the accumulation phase.

dCD(t)dt=ECD(t)(16)

Variable C_D(t) is the internal concentration at time t during the depuration phase. Parameters and variables have the same meaning as given in Table 1.

Regarding the accumulation phase, Equation (15) has a solution directly given by Equation (14), whatever the initial condition equal to C_A(t = 0) = C₀.

As a consequence, the generic solution for the accumulation phase writes as follows:

CA(t)=(etE−I)E−1Ucx+etEC0(17)

Regarding the depuration phase, Equation (16) is similar to Equation (6), with the corresponding generic solution given by Equation (7). The constant vector Ω₁ can be determined from the initial condition of the depuration phase that corresponds to the internal concentration reached at t = t_c at the end of the accumulation phase. We must therefore solve the following equation:

CA(tc)=CD(tc)(18)

Given solutions from Equations (14) and (7), we get:

CA(tc)=(etcE−I)E−1Ucx+etcEC0 and CD(tc)=etcEΩ1(19)

Hence, the constant vector Ω₁ derives from the following equation:

(etcE−I)E−1Ucx+etEC0=etcEΩ1 ⇔Ω1=(I−e−tcE)E−1Ucx+C0(20)

The generic solution for the depuration phase then writes as follows:

CD(t)=(etE−e(t−tc)E)E−1Ucx+etEC0(21)

Results

This section presents four case studies with different numbers of compartments. Two case studies concern the species Gammarus fossarum exposed to cadmium (Cd).^27,44 The third one concerns Apostichopus japonicus sea cucumber exposed to six different antibiotics;⁴⁵ the fourth one concerns species Danio rerio exposed to arsenic (As).⁴⁶ Gestin and colleagues^27,44 used several one-compartment PBK models compared with one four-compartment PBK model to gain knowledge on the accumulation and fate dynamic of Cd in and between gammarids’ organs. Subsections give the generic solutions with median parameter values for this particular case study (Table 2).

Table 2

Medians of parameters estimated by Bayesian inference from TK one-compartment models separately fitted to each organ of Gammarus fossarum exposed to dissolved Cd at 11.1 μg·L⁻¹ for seven days before being placed for 14 days under depuration conditions

Process	Organ	Parameter	Median value (in [t]⁻¹)
Uptake	Intestines	k_u,₁	1917
	Caeca	k_u,₂	1571
	Cephalons	k_u,₃	91.1
	Remaining tissues	k_u,₄	135
Elimination	Intestines	k_e,₁	0.506
	Caeca	k_e,₂	0.053
	Cephalons	k_e,₃	0.060
	Remaining tissues	k_e,₄	0.026

Cd, Cadmium; TK, toxicokinetics.

Case study with one compartment

Applying Equation (2) with only one compartment leads to two single equations:

{dcA,i(t)dt=ku,1cx−ke,1cA,i(t) when 0≤t≤tc (accumulation phase) (a)dcD,i(t)dt=−ke,1cD,i(t) when t>tc (depuration phase) (b)(22)

System of Equations (22a) and (22b) can easily be solved with the method of the separation of variables (also known as the Fourier method), leading to the following system of solutions for both accumulation and depuration phases:

{ cA,i(t)=ku,ike,icx(1−e−ke,it) when 0≤t≤tc (accumulation phase) (a)cD,i(t)=ku,ike,icx(eke,i(tc−t)−e−ke,it) when t>tc (depuration phase) (b)(23)

Getting the set of solutions (23) directly from the generic expressions given by the combination of both Equations (17) and (21) leads precisely to the same result. Indeed, with one compartment, matrix E = −k_e,i and vector U = k_u,i.

Inspired from Gestin et al.,^27,44 when considering only solid black arrows, Figure 3 highlights the target organs that can correspond to one compartment according to i: intestine (i = 1); cephalon (i = 2); caeca (i = 3); remaining tissues (i = 4). Each compartment has its own parameter pair for uptake (k_u,i) and elimination (k_e,i) rates (Table 2). Model parameters were estimated under a unified Bayesian framework.²⁰ In particular, parameters of PBK one-compartment models were fitted separately for each organ of G. fossarum exposed to dissolved Cd at 11.1 μg·L⁻¹ for seven days before being placed for 14 days under depuration conditions. Getting median parameter values as given in Table 2 allows simulating what happens within the intestines when it is connected to all other organs, for example (see Supplementary File 1 for more details).

General scheme of the multi-compartment physiologically-based kinetic model used by at the initial modeling stage when all compartments were connected to each other.

Fig. 3 General scheme of the multi-compartment physiologically-based kinetic model used by ¹ at the initial modeling stage when all compartments were connected to each other.

Parameter values as given in Table 3.

Case study with four compartments

Applying the general matrix ODE system given by the set of Equations (15) and (16) to the particular case of four compartments connected by pairs (Fig. 3) leads to the following writing:

{dCA(t)dt=Ucx+ECA(t) when 0≤t≤tc (accumulation phase) (a)dCD(t)dt=ECD(t) when t>tc (depuration phase) (b)(24)

with vectors and matrices defined as follows:

CA(t)=CD(t)=(c1(t)c2(t)c3(t)c4(t)c5(t)) U=(ku,1ku,2ku,3ku,4ku,5) E=(e1,1k1,2k1,3k1,4k1,5k2,1e2,2k2,3k2,4k2,5k3,1k3,2e3,3k3,4k3,5k4,1k4,2k4,3e4,4k4,5k5,1k5,2k5,3e5,4e5,5)(25)

and diagonal elements of matrix E defined by:

{e1,1=−ke,1−(k2,1+k3,1+k4,1) (a)e2,2=−ke,2−(k1,2+k3,2+k4,2) (b)e3,3=−ke,3−(k1,3+k2,3+k4,3) (c)e4,4=−ke,4−(k1,4+k2,4+k3,4) (d)(26)

The exact solution of the matrix ODE system (24) can be directly deduced from Equations (17) and (21):

{CA(t)=(etE−I)E−1Ucx when 0≤t≤tc (accumulation phase) (a)CD(t)=(etE−e(t−tc)E)E−1Ucx when t>tc (depuration phase) (b)(27)

Developing the matrix Equations (27a) and (27b) finally provides the following sets of four equations for both accumulation and depuration phases:

For the accumulation phase ( 0 ≤ t ≤ t_c):

{dc1(t)dt=ku,1cx−ke,1c1(t)+k1,2c2(t)+k1,3c3(t) +k1,4c4(t)−(k2,1+k3,1+k4,1)c1(t) (a)dc2(t)dt=ku,2cx−ke,2c2(t)+k2,1c1(t)+k2,3c3(t) +k2,4c4(t)−(k1,2+k3,2+k4,2)c2(t) (b)dc3(t)dt=ku,3cx−ke,3c3(t)+k3,1c1(t)+k3,2c2(t) +k3,4c4(t)−(k1,3+k2,3+k4,3)c3(t) (c)dc4(t)dt=ku,4cx−ke,4c4(t)+k4,1c1(t)+k4,2c2(t) +k4,3c3(t)−(k1,4+k2,4+k3,4)c4(t) (d)(28)

For the depuration phase (t > t_c):

{dc1(t)dt=−ke,1c1(t)+k1,2c2(t)+k1,3c3(t) +k1,4c4(t)−(k2,1+k3,1+k4,1)c1(t) (a)dc2(t)dt=−ke,2c2(t)+k2,1c,(t)+k2,3c3(t) +k2,4c4(t)−(k1,2+k3,2+k4,2)c2(t) (b)dc3(t)dt=−ke,3c3(t)+k3,1c1(t)+k3,2c2(t) +k3,4c4(t)−(k1,3+k2,3+k4,3)c3(t) (c)dc4(t)dt=−ke,4c4(t)+k4,1c1(t)+k4,2c2(t) +k4,3c3(t)−(k1,4+k2,4+k3,4)c4(t) (d)(29)

The four-compartment model in Equations (27a) and (27b) based on the matrix ODE system in Equations (24a) and (24b), with its exact solution in (28) and (29), assumes that all compartments are connected to each other by pairs and with the external medium (Fig. 3). This means that the model is considering all incoming and outgoing arrows from all compartments. This model comprises a total of 20 parameters, plus the c_x value for the exposure concentration. As given in Table 3, median parameter values were used to simulate what happens within each organ in terms of internal concentration dynamic and compared to the previous results with the four independent one-compartment PBK models. See Supplementary File 1 for details.

Table 3

Parameter estimates (expressed as medians and 95% uncertainty intervals) of the four-compartment model corresponding to solid black arrows in Figure 3 as provided by¹ in Table S6

Organ-Connection	Parameter	Median	Q_2.5%	Q_97.5%
Intestines (uptake)	k_u,₁	3342	2720	3707
Intestines (elimination)	k_e,₁	0.54	0.415	1.402
Intestines-Caeca	k₂₁	0.873	0.603	1.739
Caeca-Intestines	k₁₂	0.218	0.132	0.376
Intestines-Cephalons	k₃₁	0.059	0.034	0.124
Cephalons-Intestines	k₁₃	0.262	0.124	0.871
Intestines-remaining tissues	k₄₁	0.069	0.049	0.126
Remaining tissues-Intestines	k₁₄	0.14	0.086	0.238
Intestines	σ₁	8.974	6.469	15.28
Caeca	σ₂	17.94	13.07	26.84
Cephalons	σ₃	1.223	0.863	1.818
Remaining tissues	σ₄	1.468	1.06	2.242

The first column stands for connected organs, either to water or to the other organs; the second column is for parameter names; the following three columns are for medians, lower and upper quantiles of parameter estimates when G. fossarum was exposed to Cd = 11.1 μg·L⁻¹.

As illustrated above, our generic modeling framework allows simulating complex situations involving several compartments, their connections in pairs and/or with the exposure media. Let us now relate what we did for simulations of four one-compartment models for each organ of G. fossarum exposed to Cd – with the four-compartment model developed by Gestion et al.²⁷ Indeed, they showed that G. fossarum takes up and eliminates Cd rather quickly, with the intestines and the caeca accumulating and depurating the most compared to the cephalon and the remaining tissues. Gestin et al. also proved that a four-compartment model better describes the Cd internal contamination route than the single one-compartment model for each organ. Furthermore, they finally highlighted that the most parsimonious multi-compartments model corresponds to the solid black arrows in Figure 3.

Such a situation corresponds to a nested model within the four-compartment ODE system as given by Equations (24a) and (24b). This model thus comprises only 12 parameters whose values are listed in Table 3. Figure 4 shows the simulated kinetics within the four organs. Our four curves exactly superimpose to the four median curves provided in Figure 3 by Gestin et al. Benefiting from this exact match between our exact generic solution and what was numerically integrated by these authors before the implementation of their inference process ultimately strengthens both approaches: the numerical integration (a basic Euler integration scheme with a time-step equal to 1/10 day); and the curve plotting from the exact solution. Nevertheless, to infer parameter values from observed data, there is absolutely no doubt that the exact solution will provide much better computational performance for implementing the Monte Carlo Markov Chain simulations needed to use Bayesian inference. Readers who would like to convince themselves of this added value of our generic solving can refer to our dedicated R-package named ‘rPBK,’ officially available from the official R CRAN website (https://CRAN.R-project.org/package=rPBK ) together with a very recent related paper.⁴⁷ Thanks to the R-package ’rPBK’, we already experienced that using this exact generic solution divided at least by 100 the computation time of the whole process. Our first inference implementation required numerically integrating the original ODE system and running the MCMC algorithm to fit the numerical solution to the accumulation-depuration data sets corresponding to the different organs and tissues. The numerical integration step was highly computationally demanding. Avoiding this numerical step much faster now delivers relevant and precise parameter estimates.

Simulations of the internal concentrations within the different organs or tissues of <italic>Gammarus fossarum</italic> when exposed to an external cadmium concentration equal to 11.1 <italic>μg</italic>·<italic>L</italic><sup>−1</sup>.

Fig. 4 Simulations of the internal concentrations within the different organs or tissues of Gammarus fossarum when exposed to an external cadmium concentration equal to 11.1 μg·L⁻¹.

The four-compartment PBK solution used for simulations is given in Equation (27). Parameter values are those given in Table 2.

Case study with five compartments

Zhu and colleagues⁴⁵ studied the effect of six different antibiotics on the sea cucumber Apostichopus japonicus: sulfadiazine, trimethoprim, enrofloxacin, ofloxacin, clarithromycin, and azithromycin. All compartments (blood or organs) are internally related to the coelomic fluid. Except for the digestive tract, all compartments are externally related to seawater (Fig. 5). Based on the original paper,⁴⁵ biological parameter values were extracted to calculate the model’s coefficients used for the simulations. All these parameters are provided in Table 4.

Schematic representation of the five-compartment kinetic model used to simulate the effects of waterborne antibiotics in sea cucumbers: the coleomic fluid (<italic>CF</italic>), the body wall (<italic>BW</italic>), the mouth (<italic>MH</italic>), the respiratory tree (<italic>RT</italic>) and the digestive tract (<italic>DT</italic>).

Fig. 5 Schematic representation of the five-compartment kinetic model used to simulate the effects of waterborne antibiotics in sea cucumbers: the coleomic fluid (CF), the body wall (BW), the mouth (MH), the respiratory tree (RT) and the digestive tract (DT).

Parameters k_u,i and k_e,i, ∀i = 1,4, stand for uptake and elimination rates from or towards seawater, respectively; parameters k_i,j represent from and back exchanges between compartments i and j, ∀i, j = 1,5.

Table 4

Parameter values for the five-compartment physiologically-based kinetic model on sea cucumbers, after recalculation from the initial parameters as given by⁴⁵ in their supplementary information

Parameter	Sulfadiazine	Trimethoprim	Enrofloxacin	Ofloxacin	Clarithromycin	Azithromycin
k_u_,1	68.73	68.73	68.73	68.73	68.73	68.73
k_e_,1	68.73	5.61	4.01	4.89	5.35	5.64
k_u_,2	0.53	34.97	21.10	36.77	24.10	24.10
k_e_,12	0.06	285.76	234.84	229.46	177.87	209.72
k_u_,3	2.86	0.00	0.00	0.00	0.00	0.00
k_e_,3	0.45	68.73	68.73	68.73	68.73	68.73
k_u_,4	19.51	0.11	0.07	0.13	0.07	0.04
k_e_,4	4.11	0.51	0.59	0.88	0.22	0.14
k_u_,5	0.00	3.88	11.91	8.44	1.98	1.55
k_e_,t	0.00	0.00	0.00	0.00	0.00	0.00
k_1,2	0.06	0.08	0.09	0.11	0.07	0.04
k_2,1	0.50	0.40	0.92	0.70	0.26	0.20
k_1,3	0.43	3.27	10.75	9.29	2.20	1.93
k_3,1	2.70	0.40	1.16	1.06	0.48	0.25
k_1,4	4.11	4.00	4.94	4.17	5.52	5.50
k_4,1	19.51	27.22	32.72	29.37	28.15	35.42
k_1,5	0.48	241.35	211.96	252.56	198.06	261.31
k_5,1	2.46	29.69	29.03	26.96	25.34	25.65

Mathematical relationships between initial and recalculated parameters are detailed in the main text.

Below are the model equations written within our new generic mathematical formalism, with the following correspondence between index i and the different compartments: coleomic fluid (CF, i = 1); the body wall (BW, i = 2), the mouth (MH, i = 3), the respiratory tree (RT, i = 4) and the digestive tract (DT, i = 5).

{dCA(t)dt=Ucx+ECA(t) when 0≤t≤tc (accumulation phase) (a)dCD(t)dt=ECD(t) when t>tc (depuration phase) (b)(30)

with vectors and matrices defined as follows:

CA(t)=CD(t)=(c1(t)c2(t)c3(t)c4(t)c5(t)) U=(ku,1ku,2ku,3ku,4ku,5)E=(e1,1k1,2k1,3k1,4k1,5k2,1e2,2k2,3k2,4k2,5k3,1k3,2e3,3k3,4k3,5k4,1k4,2k4,3e4,4k4,5k5,1k5,2k5,3e5,4e5,5)(31)

and diagonal elements of matrix E defined by:

{e1,1=−ke,1−(k2,1+k3,1+k4,1+k5,1) (a)e2,2=−ke,2−k1,2 (b)e3,3=−ke,3−k1,3 (c)e4,4=−ke,4−k1,4 (d)e5,5=−k1,5 (e)(32)

The exact solution of the matrix ODE system (30) can be directly deduced from Equations (17) and (21):

{CA(t)=(etE−I)E−1Ucx when 0≤t≤tc (accumulation phase) (a)CD(t)=(etE−e(t−tc)E)E−1Ucx when t>tc (depuration phase) (b)(33)

The relationships between the initial and the recalculated parameters are the following, with the index SW referring to seawater:

ku,i=Di,SWVi and ke,i=Di,SWViPi,SW ∀i∈CF,BW,MH,RT,DTkCF,i=DCF,iVCF and ki,CF=DCF,iVCFPi,CF ∀i∈BW,MH,RT,DT

This leads to the first ODE of the five-compartment system for the sea cucumber:

dCCF(t)dt=ku,CFcw−(ke,CF+∑i≠CFki,CF) CCF(t)+∑i≠CF(kCF,iCi(t)) ∀i∈BW,MH,RT,DT(34)

with c_w = 10 μg·L⁻¹ the concentration in the seawater to which they are exposed.

kCF,i=DCF,iViPi,CF and ki,CF=DCF,iVj ∀i∈BW,MH,RT,DT

Finally, we get the four complementary equations of the system for the sea cucumber:

dCi(t)dt=ku,icw−(ke,i+kCF,i)Ci(t)+ki,CFCCF(t)(35)

Given parameter values in Table 4, we performed simulations over time for the six antibiotics and the four internal organs. As shown in Figure 6, we reproduced the same median curves as those of ⁴⁵ in Figure S10 of their supplementary information.

Predicted time course of the median internal concentrations in antibiotics (given on the top legend of each graph) under a constant exposure concentration of 10 <italic>μg</italic>·<italic>L</italic><sup>−1</sup>: in yellow the body wall (<italic>BW</italic>), in blue the mouth (<italic>MH</italic>), in red the respiratory tree (<italic>RT</italic>) and in green the digestive tract (<italic>DT</italic>).

Fig. 6 Predicted time course of the median internal concentrations in antibiotics (given on the top legend of each graph) under a constant exposure concentration of 10 μg·L⁻¹: in yellow the body wall (BW), in blue the mouth (MH), in red the respiratory tree (RT) and in green the digestive tract (DT).

Case study with six compartments

Here is a final illustration of the usefulness of our generic solution to simulate any PBK model. We identically reproduced all the simulations provided by⁴⁶ concerning exposure of zebrafish (Danio rerio) to arsenic (As). To this end, Zhang et al.⁴⁶ proposed a six-compartments model with five compartments corresponding to organs: gills (i = 2), intestine (i = 3), liver (i = 4), head (i = 5) and carcass (i = 6). The sixth compartment corresponds to blood (i = 1). The five organs were connected to blood, while the gills and intestines were also connected to the external medium (contaminated water). These assumptions are translated in Figure 7, together with the different parameters used for the corresponding PBK model (Table 5).

Fig. 7 General scheme of the six-compartment kinetic model for zebrafish exposed to arsenic (Adapted from ³).

Parameters, between-organ and/or with-water connections, numerical values, and units are given in Table 5.

Table 5

Parameter estimates for arsenic distribution in zebrafish (from⁴⁶)

Connection	ℳℒ	Value	Unit	Connection	ℳℒ	Value	Unit
Water to gills	k_u_,2	5.28 10⁻⁵	L.d⁻¹	Gill to water	k_e_,2	0.152	d⁻¹
Water to intestine	k_u_,2	1.52 10⁻⁴	L.d⁻¹	Intestine to water	k_e_,2	0.672	d⁻¹
Blood to gills	k_2,1	76.0	d⁻¹	Gill to blood	k_1,2	58.2	d⁻¹
Blood to intestine	k_3,1	27.8	d⁻¹	Intestine to blood	k_1,3	3.94	d⁻¹
Blood to liver	k_4,1	38.5	d⁻¹	Liver to blood	k_1,4	23.6	d⁻¹
Blood to head	k_5,1	97.5	d⁻¹	head to blood	k_1,5	8.48	d⁻¹
Blood to carcass	k_6,1	7.20	d⁻¹	Carcass to blood	k_1,6	0.317	d⁻¹

Notation ℳℒ stands for the estimated maximum likelihood value of the parameter.

From the generic matrix form of a PBK model we present in this paper, the model writes as follows:

Cwhole(t)=1∑i=16wi¯(∑i=16wi¯ci(t))(36)

with vectors and matrices defined as follows:

C(t)=(c1(t)c2(t)c3(t)c4(t)c5(t)c6(t))T(37)

where c_i(t), ∀i = 1,6, are the variables corresponding to internal concentrations to be simulated. Variables c_i(t) are equal to w_iC_i(t), ∀i = 1,6, with w_i the mean wet weight of blood volume (when i = 1) or fish organs (∀i = 2,6). Variables Ci(t) are the measured concentrations of As in fish organs (∀i = 2,6, expressed in μg·L⁻¹) at time t (in days). Variable c_x is the exposure concentration of As in water (expressed in μg·L⁻¹), assumed to be constant over time.

From concentrations within organs, the concentration for the whole organism can be deduced as follows:

Cwhole(t)=1∑i=16wi¯(∑i=16wi¯ci(t))(38)

According to our mathematical formalism, uptake parameters from water are included in the following vector:

U=(0ku2ku3000)T(39)

while elimination parameters towards the water, as well as parameters corresponding to organ-organ connections, are gathered together within the following matrix:

E=(−∑j=26kj,1k1,2k1,3k1,4k1,5k1,6k2,1−ke,2−k1,20000k3,10−ke,3−k1,3000k4,100−k1,400k5,1000−k1,50k6,10000−k1,6)(40)

The matrix ODE system (36) finally leads to the following exact solution:

C(t)=(etE−I)E−1Ucx(41)

where c_x is the exposure concentration in water.

The above matrix solution (41) can then be developed in order to retrieve the six-compartment PBK model as constructed by Zhang et al.⁴⁶ Denoting the exposure concentration in water (variable c_x in our modeling) by C_water, we ultimately get:

{w1¯dC1(t)dt=∑j=27k1,jCj(t)wi¯−(∑j=27kj,1)C1(t)w1¯ (a)w2¯dC2(t)dt=ku,2Cwater+k2,1C1(t)w1¯−(ke,2+k1,2)C2(t)w2¯ (b)w3¯dC3(t)dt=ku,3Cwater+k3,1C1(t)w1¯−(ke,3+k1,3)C3(t)w3¯ (c)w4¯dC4(t)dt=k4,1C1(t)w1¯−k1,4C4(t)w4¯ (d)w5¯dC5(t)dt=k5,1C1(t)w1¯−k1,5C5(t)w5¯ (e)w6¯dC6(t)dt=k6,1C1(t)w1¯−k1,6C6(t)w6¯ (f)(42)

Zhang et al. measured internal concentrations over time in each organ and blood during both accumulation and depuration phases. This data allowed them to estimate all their parameters. Using a Bayesian inference framework also provided them with quantifying the uncertainty of these parameters. Then they deliver parameter estimates as maximum likelihood (ℳℒ) values associated with a 95% credible interval, as partly reported in Table 5.

Based on parameter ℳℒ values in Table 5, using the fresh weight of the different organs as provided in Table S2 from Zhang et al.,⁴⁶ we performed simulations of the internal concentrations within each of the five organs as well as in blood (variables c_i(t), ∀i = 1,6). As shown in Figure 8, our generic solved PBK model (see Equation (41)) again exactly reproduce the median curves provided by Zhang et al. in Figure 1 (solid blue lines). Such an exact match between our curves and the authors’ ones again provides a complete check of the mathematical writing of our generic solution, together with the fact that it produces identical median curves. These results again reinforce the possibility of using this exact generic solution for the further implementation of inference processes with the guarantee of obtaining the parameter estimates much faster, avoiding the time-consuming step of the numerical integration of the ODE system.

Simulation of the total arsenic concentration in fish tissues and blood during the accumulation phase of 16 days (exposure concentration equal to 400 <italic>μg</italic>·<italic>L</italic><sup>−1</sup>) and the subsequent 16 days of depuration in a clean medium.

Fig. 8 Simulation of the total arsenic concentration in fish tissues and blood during the accumulation phase of 16 days (exposure concentration equal to 400 μg·L⁻¹) and the subsequent 16 days of depuration in a clean medium.

The PBK model used for simulations is given as a matrix solution in Equation (41) with parameter values from Table 5.

Discussion and future directions

Benefiting from the generic solution of a PBK model allows us now to envisage the continuation of this work with confidence, in particular the implementation of a Bayesian inference framework to get parameter estimates quickly and efficiently, based on a fully harmonized methodology. The next step will thus be to make freely and easily accessible this generic modeling framework, innovative both in the writing of the PBK model and in its exact resolution and the implementation of the fitting and simulation tools. Most users, who are not necessarily modeling specialists, would be willing to use more complex PBK models, especially if necessity dictates and sufficient input data is available.

Unfortunately, turnkey tools are still rare today. In line with the MOSAIC platform spirit (https://mosaic.univ-lyon1.fr ), we have already started to build a new prototype that will facilitate the use of PBK models, also for beginners, with a step-by-step workflow to first upload data and select the model to fit. By default, the complete PBK model will be automatically proposed from which users can deselect either compartment and/or exchanges between the compartment and/or with the external medium according to prior physiological knowledge. Then the user can try nested models and finally identify the most appropriate model for the question. Users will be accompanied to run the fitting process, get the results (parameter estimates and fitting plots), look at the goodness-of-fit criteria (with guidance on their interpretation), and use the model comparison criteria in case several models would have been tested. Once this calibration step is achieved, users can run simulations, compare with additional data, or plan further experiments. In the end, all expected features of a convenient help-decision service could be offered, with particular attention to further supporting the next generation tiered PBK modeling framework that could become the new paradigm in human and environmental risk assessment. Once our generic PBK modeling framework is firmly anchored in practice, we should be in the right place to consider its coupling with mechanistic models to build an utterly general modeling framework from exposure (pharmaco/toxico-kinetics) to effects (pharmaco/toxico-dynamics) on life history traits, hence defining a unifying PBKD modeling framework.

Conclusions

Our generic solving of any full PBK model comprising as many compartments as physiologically needed, as well as all potential connections between compartments and with the external medium, revealed particularly efficient in simulating diverse situations in terms of species, compounds, and purpose. The four-compartment PBK model for G. fossarum exposed to Cd highlighted the dynamic transfer of Cd among the different organs. The six-compartment PBK model for D. rerio exposed to As showed that intestines were the leading uptake site for waterborne As, instead of gills as the authors expected. Several other examples have complementary been tested (results not shown). Nevertheless, the genericity of the solution we proposed here could still be further extended in order to account for simultaneous but different routes of exposure (via water, food, sediment and/or pore water, for example), as well as several elimination processes, among which the dilution by growth would allow to be more realistic for long-lived species. Moreover, accounting for the metabolization of the parent compound could also be of great interest to better deal with organic compounds.

Supporting information

Supplementary material for this article is available at https://doi.org/10.14218/JERP.2022.00043 .

Supplementary File 1

Generic solving of a multi-compartment physiologically-based kinetic model.

(PDF)

Click here for additional data file.

Abbreviations

As:: arsenic

Cd:: cadmium

ODE:: ordinary differential equation

PBK model:: physiologically-based kinetic model

t_c:: duration of the accumulation phase

Declarations

Acknowledgement

A large part of the work benefited from the French GDR “Aquatic Ecotoxicology” framework, which aims to foster stimulating scientific discussions and collaborations for more integrative approaches. The authors would like to express their sincere gratitude to Julie KLEINE-SCHULTJANN as part of a training course during her 5^th year study at the National Institute of Applied Sciences (INSA) in Lyon (France). At last, the work presented in this paper was performed using the computing facilities of the CC LBBE/PRABI.

Funding

The authors would like to thank the BioEEnViS research federation (Biodiversity, Water, Environment, City, Health) for its financial support for collaboratively achieving this study. Moreover, this work is part of the ANR project APPROve (ANR-18-CE34-0013) for an integrated approach to propose proteomics for biomonitoring: accumulation, fate, and multi-markers (https://anr.fr/Projet-ANR-18-CE34-0013).

Conflict of interest

The authors have no conflicts of interest related to this publication.

Authors’ contributions

SC resolved the full generic PBK model, wrote the manuscript, and coordinated all discussions around this paper. OGes, AC, OGef, TLL, and CL were strongly involved in the G. fossarum case study in close collaboration with OGes’s PhD supervision, who led all the underlying experimental work. JB, DL, and VB paid particular attention to the model equation writing. VB wrote the R source code, which generated simulation curves; this R code was also checked by SC, JB, DL, and CL. All authors contributed to and agreed on the final version of the manuscript.

References

1	Paini A, Leonard JA, Kliment T, Tan YM, Worth A. Investigating the state of physiologically based kinetic modelling practices and challenges associated with gaining regulatory acceptance of model applications. Regul Toxicol Pharmacol 2017;90:104-115 View Article PubMed/NCBI

2	Madden JC, Pawar G, Cronin MTD, Webb S, Tan YM, Painie A. In silico resources to assist in the development and evaluation of physiologically based kinetic models. Computational Toxicology 2019;11:33-49 View Article

3	Thompson CV, Firman JW, Goldsmith MR, Grulke CM, Tan YM, Paini A, et al. A Systematic Review of Published Physiologically-based Kinetic Models and an Assessment of their Chemical Space Coverage. Altern Lab Anim 2021;49(5):197-208 View Article PubMed/NCBI

4	Grech A, Brochot C, Dorne JL, Quignot N, Bois FY, Beaudouin R. Toxicokinetic models and related tools in environmental risk assessment of chemicals. Sci Total Environ 2017;578:1-15 View Article PubMed/NCBI

5	Tebby C, van der Voet H, de Sousa G, Rorije E, Kumar V, de Boer W, et al. A generic PBTK model implemented in the MCRA platform: Predictive performance and uses in risk assessment of chemicals. Food Chem Toxicol 2020;142:111440 View Article PubMed/NCBI

6	Mit C, Bado-Nilles A, Daniele G, Giroud B, Vulliet E, Beaudouin R. The toxicokinetics of bisphenol A and its metabolites in fish elucidated by a PBTK model. Aquat Toxicol 2022;247:106174 View Article PubMed/NCBI

7	Wang WX, Rainbow PS. Comparative approaches to understand metal bioaccumulation in aquatic animals. Comp Biochem Physiol C Toxicol Pharmacol 2008;148(4):315-323 View Article PubMed/NCBI

8	Lautz LS, Oldenkamp R, Dorne JL, Ragas AMJ. Physiologically based kinetic models for farm animals: Critical review of published models and future perspectives for their use in chemical risk assessment. Toxicol In Vitro 2019;60:61-70 View Article PubMed/NCBI

9	Zhao P, Zhang L, Grillo JA, Liu Q, Bullock JM, Moon YJ, et al. Applications of physiologically based pharmacokinetic (PBPK) modeling and simulation during regulatory review. Clin Pharmacol Ther 2011;89(2):259-267 View Article PubMed/NCBI

10	Armitage JM, Hughes L, Sangion A, Arnot JA. Development and intercomparison of single and multicompartment physiologically-based toxicokinetic models: Implications for model selection and tiered modeling frameworks. Environ Int 2021;154:106557 View Article PubMed/NCBI

11	Paini A, Tan YM, Sachana M, Worth A. Gaining acceptance in next generation PBK modelling approaches for regulatory assessments - An OECD international effort. Comput Toxicol 2021;18:100163 View Article PubMed/NCBI

12	Fransson MN, Barregard L, Sallsten G, Akerstrom M, Johanson G. Physiologically-based toxicokinetic model for cadmium using Markov-chain Monte Carlo analysis of concentrations in blood, urine, and kidney cortex from living kidney donors. Toxicol Sci 2014;141(2):365-376 View Article PubMed/NCBI

13	Sweeney LM, MacCalman L, Haber LT, Kuempel ED, Tran CL. Bayesian evaluation of a physiologically-based pharmacokinetic (PBPK) model of long-term kinetics of metal nanoparticles in rats. Regul Toxicol Pharmacol 2015;73(1):151-163 View Article PubMed/NCBI

14	Testai E, Bechaux C, Buratti FM, Darney K, Consiglio E, Kasteel EEJ, et al. Modelling human variability in toxicokinetic and toxicodynamic processes using Bayesian meta-analysis, physiologically based modelling and in vitro systems. EFSA Journal 2021;18(4):EN-6504 View Article

15	Charles S, Ratier A, Lopes C. Generic Solving of One-compartment Toxicokinetic Models. J Explor Res Pharmacol 2021;6(4):158-167 View Article

16	Breen M, Ring CL, Kreutz A, Goldsmith MR, Wambaugh JF. High-throughput PBTK models for in vitro to in vivo extrapolation. Expert Opin Drug Metab Toxicol 2021;17(8):903-921 View Article PubMed/NCBI

17	Zhang Y, Feng J, Gao Y, Liu X, Qu L, Zhu L. Physiologically based toxicokinetic and toxicodynamic (PBTK-TD) modelling of Cd and Pb exposure in adult zebrafish Danio rerio: Accumulation and toxicity. Environ Pollut 2019;249:959-968 View Article PubMed/NCBI

18	Zhang S, Wang Z, Chen J, Xie Q, Zhu M, Han W. Tissue-Specific Accumulation, Biotransformation, and Physiologically Based Toxicokinetic Modeling of Benzotriazole Ultraviolet Stabilizers in Zebrafish (Danio rerio). Environ Sci Technol 2021;55(17):11874-11884 View Article PubMed/NCBI

19	Organisation for Economic Cooperation and Development (OECD). OECD Guidelines for the Testing of Chemicals - Test No. 211: Daphnia magna Reproduction test. Paris: OECD Publishing; 2012 View Article

Ratier A, Lopes C, Labadie P, Budzinski H, Delorme N, Quéau H, et al. A Bayesian framework for estimating parameters of a generic toxicokinetic model for the bioaccumulation of organic chemicals by benthic invertebrates: Proof of concept with PCB153 and two freshwater species. Ecotoxicol Environ Saf 2019;180:33-42 View Article PubMed/NCBI

21	Ratier A, Lopes C, Geffard O, Babut M. The added value of Bayesian inference for estimating biotransformation rates of organic contaminants in aquatic invertebrates. Aquat Toxicol 2021;234:105811 View Article PubMed/NCBI

22	Ratier A, Lopes C, Multari G, Mazerolles V, Carpentier P, Charles S. New perspectives on the calculation of bioaccumulation metrics for active substances in living organisms. Integr Environ Assess Manag 2022;18(1):10-18 View Article PubMed/NCBI

23	Charles S, Ratier A, Baudrot V, Multari G, Siberchicot A, Wu D, et al. Taking full advantage of modelling to better assess environmental risk due to xenobiotics-the all-in-one facility MOSAIC. Environ Sci Pollut Res Int 2022;29(20):29244-29257 View Article PubMed/NCBI

24	Vijver MG, Van Gestel CA, Lanno RP, Van Straalen NM, Peijnenburg WJ. Internal metal sequestration and its ecotoxicological relevance: a review. Environ Sci Technol 2004;38(18):4705-4712 View Article PubMed/NCBI

25	Brinkmann M, Eichbaum K, Kammann U, Hudjetz S, Cofalla C, Buchinger S, et al. Physiologically-based toxicokinetic models help identifying the key factors affecting contaminant uptake during flood events. Aquat Toxicol 2014;152:38-46 View Article PubMed/NCBI

26	Allen GJP, Weihrauch D. Exploring the versatility of the perfused crustacean gill as a model for transbranchial transport processes. Comp Biochem Physiol B Biochem Mol Biol 2021;254:110572 View Article PubMed/NCBI

27	Gestin O, Lacoue-Labarthe T, Coquery M, Delorme N, Garnero L, Dherret L, et al. One and multi-compartments toxico-kinetic modeling to understand metals’ organotropism and fate in Gammarus fossarum. Environ Int 2021;156:106625 View Article PubMed/NCBI

28	Paini A, Leonard JA, Joossens E, Bessems JGM, Desalegn A, Dorne JL, et al. Next generation physiologically based kinetic (NG-PBK) models in support of regulatory decision making. Comput Toxicol 2019;9:61-72 View Article PubMed/NCBI

29	Box GEP. Science and statistics. J Am Stat Assoc 1976;71(356):791-799 View Article

30	Green JW, Springer TA, Holbech H. Statistical Analysis of Ecotoxicity Studies. John Wiley & Sons; 2018 View Article

31	Chojnacka K, Mikulewicz M. Encyclopedia of Toxicology (3rd Edition). Elsevier; 2014, 456-460 View Article

32	Marigomez I. Encyclopedia of Toxicology (3rd Edition). Elsevier; 2014, 398-401 View Article

33	Ratier A, Charles S. Accumulation-depuration data collection in support of toxicokinetic modelling. Sci Data 2022;9(1):130 View Article PubMed/NCBI

34	Gobas FA, Burkhard LP, Doucette WJ, Sappington KG, Verbruggen EM, Hope BK, et al. Review of existing terrestrial bioaccumulation models and terrestrial bioaccumulation modeling needs for organic chemicals. Integr Environ Assess Manag 2016;12(1):123-134 View Article PubMed/NCBI

35	Grech A, Tebby C, Brochot C, Bois FY, Bado-Nilles A, Dorne JL, et al. Generic physiologically-based toxicokinetic modelling for fish: Integration of environmental factors and species variability. Sci Total Environ 2019;651(Pt 1):516-531 View Article PubMed/NCBI

36	Amyot M, Pinel-Alloul B, Campbell PGC, Désy JC. Total metal burdens in the freshwater amphipod Gammarus fasciatus: Contribution of various body parts and influence of gut contents. Freshwater Biology 1996;35:363-373 View Article

37	Ju YR, Chen WY, Singh S, Liao CM. Trade-offs between elimination and detoxification in rainbow trout and common bivalve molluscs exposed to metal stressors. Chemosphere 2011;85(6):1048-1056 View Article PubMed/NCBI

38	Chen W, Hoffmann AD, Liu H, Liu X. Organotropism: new insights into molecular mechanisms of breast cancer metastasis. NPJ Precis Oncol 2018;2(1):4 View Article PubMed/NCBI

39	Rocha TL, Gomes T, Pinheiro JP, Sousa VS, Nunes LM, Teixeira MR, et al. Toxicokinetics and tissue distribution of cadmium-based Quantum Dots in the marine mussel Mytilus galloprovincialis. Environ Pollut 2015;204:207-214 View Article PubMed/NCBI

40	Gust KA, Fleeger JW. Exposure-related effects on Cd bioaccumulation explain toxicity of Cd-phenanthrene mixtures in Hyalella azteca. Environ Toxicol Chem 2005;24(11):2918-2926 View Article PubMed/NCBI

41	Ingersoll CG, Brunson EL, Dwyer FJ, Ankley GT, Benoit DA, Norberg-King TJ, et al. Toxicity and bioaccumulation of sediment-associated contaminants using freshwater invertebrates: a review of methods and applications. Environmental Toxicology & Chemistry 1995;14(11):1885-1894 View Article

42	Van Geest JL, Poirier DG, Sibley PK, Solomon KR. Measuring bioaccumulation of contaminants from field-collected sediment in freshwater organisms: a critical review of laboratory methods. Environ Toxicol Chem 2010;29(11):2391-2401 View Article PubMed/NCBI

43	Kuehr S, Kaegi R, Maletzki D, Schlechtriem C. Testing the bioaccumulation potential of manufactured nanomaterials in the freshwater amphipod Hyalella azteca. Chemosphere 2021;263:127961 View Article PubMed/NCBI

44	Gestin O, Lopes C, Delorme N, Garnero L, Geffard O, Lacoue-Labarthe T. Organ-specific accumulation of cadmium and zinc in Gammarus fossarum exposed to environmentally relevant metal concentrations. Environ Pollut 2022;308:119625 View Article PubMed/NCBI

45	Zhu M, Wang Z, Chen J, Xie H, Zhao H, Yuan X. Bioaccumulation, Biotransformation, and Multicompartmental Toxicokinetic Model of Antibiotics in Sea Cucumber (Apostichopus japonicus). Environ Sci Technol 2020;54(20):13175-13185 View Article PubMed/NCBI

46	Zhang J, Tan QG, Huang L, Ye Z, Wang X, Xiao T, et al. Intestinal uptake and low transformation increase the bioaccumulation of inorganic arsenic in freshwater zebrafish. J Hazard Mater 2022;434:128904 View Article PubMed/NCBI

47	Ratier A, Baudrot V, Kaag M, Siberchicot A, Lopes C, Charles S. rbioacc: An R-package to analyze toxicokinetic data. Ecotoxicol Environ Saf 2022;242:113875 View Article PubMed/NCBI

Copyright © 2023 Authors. This is an Open Access article distributed under the terms of the Creative Commons Attribution-Noncommercial 4.0 License (CC BY-NC 4.0), permitting all non-commercial use, distribution, and reproduction in any medium, provided the original work is properly cited.