## 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^{3} 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.^{11} 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.^{24}

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.^{3} 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.^{29} 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.^{30} 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.

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]^{−1} |

k_{e,i} | elimination rate from compartment i to the external medium | [t]^{−1} |

k_{i,j} | input rate from compartment j to compartment i | [t]^{−1} |

k_{j,i} | output rate from compartment j to compartment j | [t]^{−1} |

#### 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}*:

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:

*t*) gathers all internal concentrations in compartments

*i*at time

*t*,

*i*∊ [1;

*n*]:

Vector U contains all uptake rates from the external medium at exposure concentration *c _{x}*:

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*]:

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:

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:

With Ω_{1} a vector integration constant (∊ℝ* ^{n}*), and the following definition for the matrix exponential:

##### 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:

_{1}(

*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:

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

With Ω_{2} ∊ ℝ^{n} a constant to be determined.

From an initial condition C_{wsm}(*t* = 0) = C_{0}, we finally get Ω_{2} = C_{0}, which leads to the following final particular solution of Equation (2):

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:

*i.e.*the (

*n*×

*n*) square matrix with ones on the main diagonal and zeros elsewhere, and

**E**

^{−1}the inverse matrix of

**E**.

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

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.^{31} 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).

Bioaccumulation tests are usually mid-to-long term laboratory experiments designed to identify all the potential uptake pathways, including food and waterborne exposure routes.^{32} 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*fixed by the experimental design, organisms are transferred to a clean medium for the depuration phase (for

_{c}*t*>

*t*). 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.

_{c}^{22}

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.^{36} 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:

Equation (15) is identical to Equation (2), denoting C_{A}(*t*) the internal concentration at time *t* during the accumulation phase.

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*_{0}.

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

Regarding the depuration phase, Equation (16) is similar to Equation (6), with the corresponding generic solution given by Equation (7). The constant vector Ω_{1} 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:

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

Hence, the constant vector Ω_{1} derives from the following equation:

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

## 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;^{45} the fourth one concerns species *Danio rerio* exposed to arsenic (As).^{46} 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).

Process | Organ | Parameter | Median value (in [t]^{−1}) |
---|---|---|---|

Uptake | Intestines | k_{u,}_{1} | 1917 |

Caeca | k_{u,}_{2} | 1571 | |

Cephalons | k_{u,}_{3} | 91.1 | |

Remaining tissues | k_{u,}_{4} | 135 | |

Elimination | Intestines | k_{e,}_{1} | 0.506 |

Caeca | k_{e,}_{2} | 0.053 | |

Cephalons | k_{e,}_{3} | 0.060 | |

Remaining tissues | k_{e,}_{4} | 0.026 |

### Case study with one compartment

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

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:

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*) rates (Table 2). Model parameters were estimated under a unified Bayesian framework.

_{e,i}^{20}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*

^{−1}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).

### 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:

**E**defined by:

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

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}*):

For the depuration phase (*t* > *t _{c}*):

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.

Organ-Connection | Parameter | Median | Q_{2.5%} | Q_{97.5%} |
---|---|---|---|---|

Intestines (uptake) | k_{u,}_{1} | 3342 | 2720 | 3707 |

Intestines (elimination) | k_{e,}_{1} | 0.54 | 0.415 | 1.402 |

Intestines-Caeca | k_{21} | 0.873 | 0.603 | 1.739 |

Caeca-Intestines | k_{12} | 0.218 | 0.132 | 0.376 |

Intestines-Cephalons | k_{31} | 0.059 | 0.034 | 0.124 |

Cephalons-Intestines | k_{13} | 0.262 | 0.124 | 0.871 |

Intestines-remaining tissues | k_{41} | 0.069 | 0.049 | 0.126 |

Remaining tissues-Intestines | k_{14} | 0.14 | 0.086 | 0.238 |

Intestines | σ_{1} | 8.974 | 6.469 | 15.28 |

Caeca | σ_{2} | 17.94 | 13.07 | 26.84 |

Cephalons | σ_{3} | 1.223 | 0.863 | 1.818 |

Remaining tissues | σ_{4} | 1.468 | 1.06 | 2.242 |

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*.^{27} 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.^{47} 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.

### Case study with five compartments

Zhu and colleagues^{45} 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,^{45} biological parameter values were extracted to calculate the model’s coefficients used for the simulations. All these parameters are provided in Table 4.

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 |

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).

**E**defined by:

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

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

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

*c*= 10

_{w}*μg*·

*L*

^{−1}the concentration in the seawater to which they are exposed.

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

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 ^{45} in Figure S10 of their supplementary information.

### 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^{46} concerning exposure of zebrafish (*Danio rerio*) to arsenic (As). To this end, Zhang *et al*.^{46} 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).

Connection | ℳℒ | Value | Unit | Connection | ℳℒ | Value | Unit |
---|---|---|---|---|---|---|---|

Water to gills | k_{u}_{,2} | 5.28 10^{−5} | L.d^{−1} | Gill to water | k_{e}_{,2} | 0.152 | d^{−1} |

Water to intestine | k_{u}_{,2} | 1.52 10^{−4} | L.d^{−1} | Intestine to water | k_{e}_{,2} | 0.672 | d^{−1} |

Blood to gills | k_{2,1} | 76.0 | d^{−1} | Gill to blood | k_{1,2} | 58.2 | d^{−1} |

Blood to intestine | k_{3,1} | 27.8 | d^{−1} | Intestine to blood | k_{1,3} | 3.94 | d^{−1} |

Blood to liver | k_{4,1} | 38.5 | d^{−1} | Liver to blood | k_{1,4} | 23.6 | d^{−1} |

Blood to head | k_{5,1} | 97.5 | d^{−1} | head to blood | k_{1,5} | 8.48 | d^{−1} |

Blood to carcass | k_{6,1} | 7.20 | d^{−1} | Carcass to blood | k_{1,6} | 0.317 | d^{−1} |

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

*c*(

_{i}*t*), ∀

*i*= 1,6, are the variables corresponding to internal concentrations to be simulated. Variables

*c*(

_{i}*t*) are equal to

*(*w

_{i}C_{i}*t*), ∀

*i*= 1,6, with

*the mean wet weight of blood volume (when*w

_{i}*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*

^{−1}) at time

*t*(in days). Variable

*c*is the exposure concentration of As in water (expressed in

_{x}*μg*·

*L*

^{−1}), assumed to be constant over time.

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

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

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

*c*is the exposure concentration in water.

_{x}The above matrix solution (41) can then be developed in order to retrieve the six-compartment PBK model as constructed by Zhang *et al*.^{46} Denoting the exposure concentration in water (variable *c _{x}* in our modeling) by

*C*, we ultimately get:

_{water}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*.,^{46} 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.

## 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)

## 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 (

### 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.