# Derivation and validation of compartment models: implications for dynamic imaging

Jérôme Kowalski, Lorenzo Sala, Dirk Drasdo, Irene Vignon-Clementel

► To cite this version:

Jérôme Kowalski, Lorenzo Sala, Dirk Drasdo, Irene Vignon-Clementel. Derivation and validation of compartment models: implications for dynamic imaging. IEEE Transactions on Medical Imaging, In press, 10.1109/TMI.2025.3630705. [hal-05350568](https://HAL Id: hal-05350568)

Submitted on 6 Nov 2025

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
# Derivation and validation of compartment models: implications for dynamic imaging

Jérôme Kowalski, Lorenzo Sala, Dirk Drasdo and Irene Vignon-Clementel

**Abstract**— Compartment models are mathematical representations of the transport of a chemical substance in the human body, assuming uniform concentration within each compartment, which corresponds to a distinct body part. Widely used in pharmaceutical and medical imaging through pharmacokinetic and tracer kinetic (TK) models, they are formulated as systems of time-dependent ordinary differential equations. However, these models do not account for the spatial dependence of physiological processes within individual compartments. This work aims to explore the physical interpretation of these models through mathematical derivation and to analyze their limitations. Three TK models are derived from more complex models regarding the physiological processes represented. The derivation introduces the hypotheses relevant to these processes. The most notable is the *well-mixed hypothesis*, resembling concentration homogeneity in the region of interest. The hypotheses are numerically tested by simulating the more complex model and evaluating the reduced models' residuals. For a given tracer molecule, the validity of TK models is found to decrease as voxel size increases. This work proposes an algorithm to compute the maximum voxel size, above which the reduced models' residuals exceed a chosen error threshold. For a 3% error, voxel sizes larger than 250 µm are found critical for highly permeable vessels, while sizes smaller than 4 mm induce less than 3 % error in non-permeable vessels. Tested on literature-derived datasets, these findings identify diffusion as the underlying phenomenon leading to well-mixed compartments, and highlight the need for high spatial resolution imaging to improve TK parameters ground truth estimation.

**Index Terms**— Applicability range, Mathematical derivation, Numerical simulation, Transport models, Dynamic medical imaging, Contrast agent

## I. INTRODUCTION

DYNAMIC imaging plays a crucial role in measuring perfusion, vessel permeability, or microvascular density, especially in oncology for tumor growth and treatment assessment. This technique involves injecting a contrast agent into the patient's bloodstream and capturing a series of images of the target organ. One then observes the temporal changes

in image intensity caused by the contrast agent. Image intensity is interpreted as being proportional to the contrast agent concentration. To track its dynamics over time, a single mathematical concentration transport model is employed for the entire image. The model's parameters are then estimated individually for each region of interest (ROI) – whether a single voxel or a cluster of voxels – by fitting the model solution to the time dynamics observed in the selected ROI. The mathematical models in question generally fall into one of the following frameworks:

1) *Pharmacokinetic (PK) models:* They investigate the transport and dynamics of chemical species within the body, from their point of administration to their excretion, including distribution to organs and cellular metabolization. PK models are popular in biomedical research for their precision and simplicity, in pharmaceutical and toxicity studies [1]–[4], and dynamic imaging [5]–[10]. Originally developed to track a molecule's journey through the human body and its uptake by tissues [11], PK models aim to measure drug characteristics or predict the body's response when combined with pharmacodynamics. They rely on *compartments*, which are representations of a region where only drug average concentration is considered, under the *well-stirred* or *well-mixed* hypothesis [12]. PK model parameters are therefore aggregated representations of the physiological quantities at play, usually involving basic mathematical forms. PK studies are typically conducted over a few hours to days.

2) *Tracer Kinetic (TK) models:* They specifically focus on the transport of contrast agents (or *tracers*) in the body. Advances in both high-spatial-resolution medical imaging and contrast agent development have increased the importance of dynamic contrast-enhanced imaging for diagnosis (DCE-MRI, DCE-CT, DCE-US)[13], [14]. Employing PK equations at shorter time and space scales (minutes vs. hours, single voxel vs. entire organ), TK studies are indeed of the order of a radiological examination: a few seconds to several minutes [12], [15].

3) *Physical transport models:* Derived from conservation principles, these models provide a spatio-temporal representation of mass transport dynamics. Though not limited to biofluids – with applications in geophysics [16] and plasma physics [17] – they are particularly valuable in biomedical contexts for targeting precise spatial resolution, such as studying spatial anisotropies in the region of interest. Examples include platelet or fat deposition in blood vessels [18], aerosol deposition in airways [19], transport under complex flow

Submitted November 6, 2025. We acknowledge funding from the European Research Council (ERC) under the European Union's Horizon 2020 research and innovation program (Grant agreement No. 864313).

J. Kowalski, D. Drasdo, and I. Vignon-Clementel are with Inria Saclay, Palaiseau, France (e-mail: jerome.kowalski@inria.fr; dirk.drasdo@inria.fr; irene.vignon-clementel@inria.fr). L. Sala is with Université Paris-Saclay, INRAE, Jouy-en-Josas, France (e-mail: lorenzo.sala@inrae.fr)
regimes [20], and heterogeneous oxygen distribution in tumors [21]. However, developing accurate patient-specific models is challenging due to the intricate spatial structure of the region of interest, and the complex physiological and chemical interactions within the body. Dynamic imaging data usually cover large regions, making it difficult to resolve vascular structures without assuming a particular shape [22], [23]. Some studies have combined PK (or TK) equations with physical approaches. For example, TK Field Theory [24] proposes a spatially resolved tracer-kinetic model with interconnected compartments. Another model, the axially distributed model [25], aims to predict transport dynamics more realistically than conventional PK models. However, these models have yet to be applied to imaging data for parameter estimation tasks.

While PK and TK models offer fast and valuable insights into drug transport into blood and the surrounding tissues, their simplification compared to three-dimensional (3D) physical models can introduce significant limitations. These models rely on a strong hypothesis proposed by [11], which assumes well-mixed compartments: a hypothesis precisely revisited in this work. Our goal is to critically analyze the physical meaning and validity of this hypothesis (and others) through analytical derivation and numerical simulation. For the first time, we mathematically derive the simpler models from the physical equations of transport and analyze them under various parameter conditions, thereby exploring the validity of the derivation hypotheses with respect to the mass conservation equations. We consider this work as an extension to prior research focused on further reducing the TK models [14]. We validate our findings with numerical simulations and discuss the implications of these limitations for the interpretation of dynamic imaging data. Finally, we present an innovative algorithm to determine the maximum compartment or ROI size for each mathematical model to ensure minimal interpretation errors, and compare the presented models with imaging data from the literature.

## II. FROM PHYSICAL TRANSPORT MODELS TO TRACER KINETICS

The studied system is represented in figure 1. It consists in a tissue volume with two phases: the blood vessels, which contain plasma (indexed with $p$), and the extravascular space (which includes cells, interstitial space, etc; indexed with $e$). The total volume $V_{tot}$ is the sum of the two phases. In practice, the entire volume may not be accessible to the tracer (see remark II.3). Blood and the tracer (diluted in blood) enter through the plasma phase at $x = 0$ and exit the volume at $x = L$. Transport is driven by an arterial input function $c_{in}(t)$ entering the plasma phase at $x = 0$. We assume symmetric boundary conditions (meaning that the concentration is locally uniform across the border) on all other exterior boundaries. From 3D equations, we first present a derivation of the one-dimensional (1D) models that serve as a basis for all TK models presented later.

### A. One-dimensional physical models

We start by defining the total cross-section area of phase $i$ $A_i(x)$ at position $x$. To reduce the number of dimensions of the

<figure><img src="image_2.png" /><figcaption>Fig. 1: Left: Schematic representation of the system, composed of blood vessels (in red) and the extravascular space (in light blue). The blood velocity $\mathbf{u}$ is a vector field defined at each point $\mathbf{x} = (x, y, z)$ inside the blood vessels. The $x$ axis is defined along the main direction of blood flow $\mathbf{e}_x$. Right: planar $(\mathbf{e}_y, \mathbf{e}_z)$ representation of the volume $V_{tot}$ cross-section, normal to vector $\mathbf{e}_x$ at position $x$. The intersection of this plane with the blood vessels defines the cross-section of the blood vessels $\Sigma_p(x)$ (2D domain, union of all the red surfaces) and their interface with the extravascular 2D space $\partial\Sigma(x)$ (solid black lines). $\Sigma_e(x)$ is the cross-section of the extravascular phase. The concentration $c(\mathbf{x})$ of a chemical species of isotropic diffusivity $D(\mathbf{x})$ is followed inside the system (box).</figcaption></figure>

problem, we study the evolution over time $t$ of the quantities $\mathbf{f}(\mathbf{x}, t)$ (vector or scalar field), averaged over the cross-section:

$$A_i(x, t) = \int_{\Sigma_i(x, t)} d\sigma; \quad \bar{\mathbf{f}}_i(x, t) = \frac{1}{A_i(x, t)} \int_{\Sigma_i(x, t)} \mathbf{f}(\mathbf{x}, t) d\sigma. \quad (1)$$

$d\sigma = dydz$. Let $c(\mathbf{x}, t)$ be the tracer concentration. For the sake of conciseness, we now focus on a single phase, removing the index $i$. The quantity of molecules per unit length is $n_l = \int_{\Sigma} c d\sigma$. Applying the Reynolds transport theorem on the surface $\Sigma$, the time derivative of $n_l$ obeys the following equation:

$$\frac{\partial n_l}{\partial t} = \frac{\partial}{\partial t} \left( \int_{\Sigma} c d\sigma \right) = \int_{\Sigma} \frac{\partial c}{\partial t} d\sigma + \int_{\partial \Sigma} c(\mathbf{u}_b \cdot \mathbf{n}) dl, \quad (2)$$

where $dl$ is the line element on the boundary, $\mathbf{u}_b$ is the velocity of the boundary, and $\mathbf{n}$ is the local vector normal to the boundary $\partial\Sigma$. In addition, inside a solvent, if the solute concentration is small enough, the concentration follows the transport equation in 3D:

$$\frac{\partial c}{\partial t} = -\nabla \cdot (c\mathbf{u}) + \nabla \cdot (D\nabla c) + \mathcal{R}(c), \quad (3)$$

where $\mathbf{u}$ is the solvent velocity (blood or interstitial fluid), $D$ the diffusion coefficient of the chemical (tracer), and $\mathcal{R}(c)$ gathers the phenomena that act as a source or a sink of matter (such as chemical reactions, for instance). Since $n_l = A\bar{c}$ thanks to (1), and replacing the partial derivative in time $\frac{\partial c}{\partial t}$ in the right-hand side of (2) by its expression in (3), (2) reads:

$$\frac{\partial A\bar{c}}{\partial t} = \int_{\Sigma} \nabla \cdot (D\nabla c - c\mathbf{u}) d\sigma + A\overline{\mathcal{R}(c)} + \int_{\partial\Sigma} c(\mathbf{u}_b \cdot \mathbf{n}) dl. \quad (4)$$
Recall that for an arbitrary vector function $\mathbf{f} = (f_x, f_y, f_z)$, the divergence theorem on a surface $\Sigma$ normal to $\mathbf{e}_x$ writes

$$ \int_{\Sigma} \nabla \cdot \mathbf{f} d\sigma = \frac{\partial A \overline{f_x}}{\partial x} + \int_{\partial \Sigma} \mathbf{f} \cdot \mathbf{n} dl. \quad (5) $$

Introducing the divergence theorem (5) with $\mathbf{f} = D\nabla c - c\mathbf{u}$ into (4) leads to the 1D transport equation:

$$ \frac{\partial A \bar{c}}{\partial t} = \frac{\partial}{\partial x} \left( A D \frac{\partial c}{\partial x} \right) - \frac{\partial A c u_x}{\partial x} + A \bar{R} \\ + \int_{\partial \Sigma} \left( D \frac{\partial c}{\partial \mathbf{n}} + c(\mathbf{u}_b - \mathbf{u}) \right) \cdot \mathbf{n} dl. \quad (6) $$

$u_x$ is the fluid velocity component along $\mathbf{e}_x$. Equation (6) is the leading reduced equation for transport in 1D. The left-hand side (LHS) of the equation specifies the variation of linear density at a given cross-section normal to a considered flow direction $\mathbf{e}_x$ per unit time. The first two terms in the right-hand side (RHS) are the diffusion in the direction of $\mathbf{e}_x$ and the advective transport of mass in the direction of $\mathbf{e}_x$. Hence, they describe the transport fluxes in the direction normal to the cross-section. The third term in the RHS represents sources and sinks; the integral represents transport (diffusion and advection) within the plane of the cross section, relative to the velocity of the vessel wall $\mathbf{u}_b$.

The next subsections present the 1D models that are derived as a basis for the further-reduced models, and the derivation of most TK parameters. Below we sometimes write out the arguments of $c$ for clarity.

1) **List of hypotheses to study specific 1D cases:** First a full list of hypotheses ("H.X", X=1,2,... refers to hypothesis X) is given together with their mathematical expression. Then specific cases are considered: equation (6) is modified by applying hypotheses reflecting each specific case.

### H.1: No reaction

Reactions effects are negligible in phase i, $R_i = 0$

### H.2: Non-deformable walls

The vessels walls are fixed $\mathbf{u}_b = 0$

### H.3: Uniform diffusion

In phase i, $\partial_y D_i = \partial_z D_i = 0$

### H.4: Parabolic profile

In a cylindrical pipe of radius $a$, in locally cylindrical coordinates $\mathbf{x} = (x, r, \theta)$, functions of $(x, t)$ $u_{x,0}$, $u_{x,1}$, $c_0$, $c_1$ exist such that $u_x(x, r, \theta, t) = u_{x,0}(x, t) + u_{x,1}(x, t)(r^2/a^2 - 1)$ and $c(x, r, \theta, t) = c_0(x, t) + c_1(x, t)(r^2/a^2 - 1)$.

### H.5: Impermeable walls

In phase i, the flux of $c_i$ writes $(\mathbf{u}c_i - D\frac{\partial c_i}{\partial \mathbf{n}}) \cdot \mathbf{n} = 0$

### H.6: Uniform cross-section

If applied to phase i, $\partial_x A_i = 0$

### H.7: Low velocity

If applied to phase i, $\mathbf{u}_i = 0$

### H.8: Permeable walls

At the interface between i and j the flux of $c_i$ writes $(\mathbf{u}_i c_i - D_i \frac{\partial c_i}{\partial \mathbf{n}_{i,j}}) \cdot \mathbf{n}_{i,j} = P(c_i - c_j)$ With a linear permeability $P$ (P has the unit of a velocity).

### H.9: Negligible diffusion

In phase i, $D_i = 0$

### H.10: Uniform velocity

In phase i, $\partial_x u_i = 0$

### H.11: Well-mixed compartment

In phase i, there is a constant $\lambda \in [0, 1]$, such that $\bar{c}_i(t) = \frac{1}{L} \int_0^L c_i(x, t) dx = c_i(0, t) + \lambda(c_i(L, t) - c_i(0, t))$

### H.12: Concentration homogeneity

In phase i, $\partial_x c_i = 0$

### Remark II.1

Mathematically, H.11 is valid whenever $c_i$ is independent of time and is a monotonous function of position. For a phase i, it is mathematically weaker than H.12.

### H.13: Adiabatic approximation

At the interface between i and j the flux of $c_i$ writes $(\mathbf{u}_i c_i - D_i \frac{\partial c_i}{\partial \mathbf{n}_{i,j}}) \cdot \mathbf{n}_{i,j} = L P(c_i - c_j) \delta(L - x)$

2) **Single-phase advection-diffusion:** Here, the tracer is transported in a representative blood vessel and cannot escape (or extravasate). This is usually the case in larger vessels (arteries, veins) or smaller vessels (arterioles, venules, capillaries) with very low or no exchange with the extravascular space (negligible permeability). In absence of reaction (H.1), and for non-deformable (H.2), impermeable (H.5) walls, the two right-most terms in (6) vanish. In case of a parabolic velocity profile (H.4), corresponding to steady unidirectional flow of the solvent, in a long tube of circular section of radius $a$ [26], this leads to the advection-diffusion equation in 1D:

$$ \frac{\partial A \bar{c}}{\partial t} + \frac{\partial A u \bar{c}}{\partial x} = \frac{\partial}{\partial x} \left( D A \frac{\partial \bar{c}}{\partial x} \right). \quad (7) $$

$u$ is an affine function of $\bar{u}_x$ (see remark II.2). It is equal to $\bar{u}_x$ for a uniform concentration cross-section profile (degenerate parabola; $\alpha = 1$). $\bar{c}(x, t)$ is the cross-sectional average of the concentration and will from now on be referred to as $c(x, t)$.

### Remark II.2

If $c_0$ (as defined in H.4) is written $c_0 = \alpha \bar{c}$, then H.3 and H.4 give $A \bar{u}_x \bar{c} - A D \partial_x \bar{c} = A u \bar{c} - A D \partial_x \bar{c}$, with $Au = (A \bar{u}_x \frac{4-\alpha}{3} - (1-\alpha)A'D)$.

3) **Two-phase exchange model:** Blood vessels are now assumed to be permeable for the tracer that may leave the plasma and enter another phase. The two phases – plasma with
accessible volume fraction $\phi_p$, extravascular with accessible volume fraction $\phi_e$ – share the interface $\partial\Sigma$ along the $x$ axis.

### Remark II.3

In practice, the total volume $V_{tot}$ may not be entirely accessible to the tracer. When the tracer cannot penetrate the cells' membrane, for example. $\phi_p$ and $\phi_e$ thus represent the accessible volume fractions. In this case, $\phi_p + \phi_e \le 1$; otherwise, $\phi_p + \phi_e = 1$.

We assume linear passive exchange (H.8), with a permeability $P$. By then assuming H.1-H.3, uniform cross-section (H.6) in both phases, H.4 in the plasma and low fluid velocity (H.7) in the extravascular phase, (6) becomes for the two phases

$$A_p \frac{\partial c_p}{\partial t} + u_p A_p \frac{\partial c_p}{\partial x} = D_p A_p \frac{\partial^2 c_p}{\partial x^2} + \mathcal{L}_{\partial\Sigma} \mathcal{P}(c_e - c_p), \\ A_e \frac{\partial c_e}{\partial t} = D_e A_e \frac{\partial^2 c_e}{\partial x^2} + \mathcal{L}_{\partial\Sigma} \mathcal{P}(c_p - c_e).$$

$\mathcal{L}_{\partial\Sigma} = \int_{\partial\Sigma} dl$ is the wall perimeter at the interface. To compare the current mathematical framework with existing works, we write $L_p = \phi_p \frac{V_{tot}}{A_p} = L$ the length of the voxel or ROI in $x$-direction, $F_p = \frac{u_p A_p}{V_{tot}} = \phi_p \frac{u_p}{L_p}$ the plasma flow rate, $P = \frac{P}{V_{tot}}$ the permeability coefficient, and $S = L_p \mathcal{L}_{\partial\Sigma}$ the surface area of the interface. This leads to the model of transport in two phases, with passive exchange:

$$\phi_p \frac{\partial c_p}{\partial t} = PS(c_e - c_p) - F_p L_p \frac{\partial c_p}{\partial x} + D_p \phi_p \frac{\partial^2 c_p}{\partial x^2} \quad (8) \\ \phi_e \frac{\partial c_e}{\partial t} = PS(c_p - c_e) + D_e \phi_e \frac{\partial^2 c_e}{\partial x^2},$$

$c_i = c_i(x, t)$, $i \in \{p, e\}$. The product $PS$ is called the permeability-surface area and considered as one parameter.

Equation (8) will be referred to as the full order model (FOM) as it introduces fewer hypotheses compared to the models presented in the following. Compared to other works, the current analysis explicitly states the derivation hypotheses to relate the TK parameters with the physical quantities. Figure 2 provides a schematic outline of the derivations, and table I summarizes the found relation between TK parameters and physical quantities. Section II-B analytically analyzes the asymptotic regimes in the single-phase case, and section II-C offers a new derivation for some TK models, which ensures the respect of the physical processes involved.

## B. The (one) compartment model

This section focuses on equation (7). However, in applications such as PK or TK, when the full architecture of the observed vascular network is unknown, the compartment model (1C) derived below is preferred to represent transport dynamics in a single phase:

$$\frac{dc}{dt}(t) = F(c_{in}(t) - c(t)), \quad (9)$$

In this section, we present the hypotheses required to mathematically derive this compartment model, before exploring the extreme regimes of the advection-diffusion equation.

<figure><img src="image_3.png" /><figcaption>Fig. 2: Schematic representation of the reduced models and clickable links to the required hypotheses. The shade of green qualitatively represents the contrast agent concentration.</figcaption></figure>

1) About the transfer function: Historically, the compartment model is derived from (7) thanks to H.9-H.11 [11]. We integrate (7) over the length $L$ of the ROI. We assume no diffusion (H.9), uniform cross-section area and velocity (H.6 and H.10). These allow us to factorize the difference in concentration between the inlet and outlet. Finally, assuming H.11 sets the average concentration $\bar{c}$ as the state variable and introduces the constant $\lambda$ to obtain the following

$$\frac{d\bar{c}}{dt} = \frac{u}{\lambda L} (c_{in} - \bar{c}).$$

Here, $c_{in}(t) = c(0, t)$ and $\bar{c}(t) = \frac{1}{L} \int_{0}^{L} c(x, t)dx$.

### Remark II.4

Hypothesis H.10 may be suppressed by adding a correction term for the input concentration:

$$\tilde{c}_{in} = \left( (1 - \lambda) + \lambda \frac{u(0)}{u(L)} \right) c_{in}, \quad \text{leading to} \\ \frac{d\bar{c}}{dt} = \frac{u(L)}{\lambda L} (\tilde{c}_{in} - \bar{c}).$$

Now we aim to verify the validity of hypotheses H.9 & H.11 and suggest a value for parameter $\lambda$. We assume uniform velocity (H.10) and non-zero diffusion. We impose a temporal concentration profile at the inlet and a Neumann boundary condition at the outlet. Equation (7) translates into:

$$\left\{ \begin{array}{l} \frac{\partial c}{\partial t} + u \frac{\partial c}{\partial x} - D \frac{\partial^2 c}{\partial x^2} = 0, \\ c(x = 0, t) = c_{in}(t), \\ \partial_x c(x = L, t) = 0. \end{array} \right. \quad (10)$$
An analytical solution of system (10) for in-depth analysis of different regimes can be obtained by a Fourier transform. For this purpose, we introduce the characteristic time $\tau$, the reduced time $\theta = \frac{t}{\tau}$, the reduced frequency $\nu$, and the reduced length $\chi = \frac{x}{L}$. We nondimensionalize (10) and study its behavior in the frequency domain:

$$ egin{cases} i\nu \frac{\tau_d}{\tau} \hat{h} + P_e \partial_\chi \hat{h} - \partial_\chi^2 \hat{h} = 0, \\ \hat{h}(0, \nu) = 1, \\ \partial_\chi \hat{h}(1, \nu) = 0, \end{cases} \quad (11) $$

### Remark II.5

For a function $f: \theta \mapsto f(\theta)$, we define $\hat{f}(\nu) = \mathcal{F}(f)(\nu) = \int_{-\infty}^{+\infty} f(\theta)e^{i\nu\theta}d\theta$ the Fourier transform of $f$. A consequence is $\int_{-\infty}^{+\infty} \theta f(\theta)d\theta = i\hat{f}'(0)$

which introduces the transfer function $h(\chi, \theta)$ such that

$$ c(\chi L, \theta \tau) = \int_{\theta^*} c_{\text{in}}(\theta^* \tau) h(\chi, \theta - \theta^*) d\theta^*, \quad (12) $$

its Fourier transform in time $\hat{h}(\chi, \nu)$ (such that $\hat{c}(\chi L, \nu) = \hat{h}(\chi, \nu)\hat{c}_{\text{in}}(\nu)$), the advection time $\tau_a = \frac{L}{u}$, the diffusion time $\tau_d = \frac{L^2}{D}$, and the Péclet number $P_e = \frac{\tau_d}{\tau_a} = \frac{uL}{D}$. For high $P_e$, molecular transport is advection-driven; for low $P_e$, it is diffusion-driven.

Solving (11) as a 2nd-order ODE in $\chi$ for each $\nu$, leads to the following expression of the transfer function:

$$ \hat{h}(\chi, \nu) = \frac{(P_e + 2k)e^{k(1-\chi)} - (P_e - 2k)e^{-k(1-\chi)}}{(P_e + 2k)e^k - (P_e - 2k)e^{-k}} e^{\frac{P_e}{2}\chi} \quad (13) $$

$k^2 = (\frac{P_e}{2})^2 + i\nu \frac{\tau_d}{\tau}$ bears the signature of the advection and diffusion regimes.

2) **Biomarkers and measurables:** As $\tau_a$ and $\tau_d$ are difficult to measure from the temporal evolution of the concentration, it is more common to study the mean transit time, defined as the average time it takes for a particle to exit the system after entering it. This is mathematically defined as [27]

$$ \tau_{t,1} = \tau \int_{-\infty}^{+\infty} \theta h(1, \theta) d\theta. \quad (14) $$

We refer to (12) for the definition of $h$. Applying remark II.5 to the transfer function (13), the mean transit time writes

$$ \tau_{t,1} = i\tau \partial_{\nu} \hat{h}(1, 0) = \tau_d \frac{e^{-P_e} - (1 - P_e)}{P_e^2} \quad (15) $$

Equation (15) translates into the advection time at high, and the diffusion time at low Péclet numbers, respectively, while remaining always smaller than the two transport characteristic times $\tau_a$ and $\tau_d$. Practically, if one has access to the inlet and outlet concentration signals – respectively $c_{\text{in}}(t)$ and $c_{\text{out}}(t) = c(L, t)$ – one can extend (14) and compute the mean transit time as the difference between the average exit time and the average entry time (16) (see appendix A for derivation). This

formulation eliminates the need for deconvolution tools in a practical setting.

$$ \tau_{t,1} = \frac{\int_t c_{\text{out}}(t) dt}{\int_t c_{\text{out}}(t) dt} - \frac{\int_t c_{\text{in}}(t) dt}{\int_t c_{\text{in}}(t) dt} \quad (16) $$

3) **Asymptotic behaviors:** As presented earlier, the Péclet number gives rise to two asymptotic regimes for (7):

a) *Time translation at high Péclet number (H.9 leads to this setting):* the transfer function writes

$$ \hat{h}(\chi, \nu) = e^{-i\nu\chi\frac{\tau_d}{\tau}} $$

### Remark II.6

The Neumann boundary condition is still preserved when considering higher-order terms that vanish as $P_e \to +\infty$.

The approximation results in $c(x, t) = c_{\text{in}}(t - \frac{x}{L}\tau_a)$. The mean transit time is $\tau_{t,1} = \tau_a$. Besides, since $\bar{c}(t) = \frac{1}{\tau_a} \int_{t-\tau_a}^{t} c_{\text{in}}(\tau) d\tau$ and $c_{\text{out}}(t) = c_{\text{in}}(t - \tau_a)$, $\bar{c}$ cannot equal $(1 - \lambda) c_{\text{in}} + \lambda c_{\text{out}}$ for any given $\lambda \in [0, 1]$ in the general case, showing that hypotheses H.9 and H.11 are incompatible.

b) *Compartment-like behavior at low Péclet number:* In this case $D \gg uL$. Under some frequency conditions (see appendix B for more details), this leads to

$$ \frac{d\bar{c}}{dt} = \frac{3}{\tau_d} (c_{\text{in}} - \bar{c}) \quad (17) $$

Comparing (17) with (9), we observe that for low Péclet number and slowly-varying input signals, the advection-diffusion equation can be approximated by the compartment model with an extraction rate $F = \frac{3}{\tau_d}$. Besides, we can show that $\lambda$ exists and equals $\frac{2}{3}$ in this case (see appendix B). This shows that a low Péclet number is a sufficient condition for hypothesis H.11 to be valid. However, historically, when assuming H.9 and H.11, the extraction rate in equation (9) is $F = \frac{u}{\lambda L} = \frac{1}{\lambda \tau_a}$: this duality may lead to a false interpretation of this constant as advection-driven rather than diffusion-driven.

## C. Two-phase models

As references [12] and [14] already present the context and applications of these models, this section focuses on deriving the classical TK two-phase models, (8) being the starting point. Then, we establish the framework for the numerical analysis of the models (section III).

1) **Derivation:** The first three derivations lead to well-known TK models. To our knowledge, the derivation of the 2CXM, TH and AATH models as well as the K'trans parameter (as presented in the following section) has not yet been established.

a) *Distributed Parameter (DP) model:* The DP model is the most general formula for tracer-kinetics models. From equation (8), negligible diffusion (H.9) is assumed in both phases. Although already discussed in section II-B.3, this hypothesis will be examined in more detail in section III. Equation (8) then reads:

$$ \begin{aligned} \phi_p \frac{\partial c_p}{\partial t}(x, t) &= PS(c_e - c_p) - F_p L_p \frac{\partial c_p}{\partial x} \\ \phi_e \frac{\partial c_e}{\partial t}(x, t) &= PS(c_p - c_e) \end{aligned} \quad (18) $$
<table><caption>TABLE I: Summary of the main TK parameters as a function of the system's physical parameters.</caption><thead><tr><th>Parameter</th><th>PS</th><th>F<sub>p</sub></th><th>K<sup>trans</sup></th><th>P<sub>e</sub></th></tr></thead><tbody><tr><td>Formula</td><td><sup>PS</sup>/<sub>P<sub>p</sub>L<sub>p</sub>L<sub>tot</sub></sub></td><td><sup>φ<sub>p</sub>u<sub>p</sub></sup>/<sub>L<sub>p</sub></sub></td><td><sup>PS · F<sub>p</sub></sup>/<sub>PS + F<sub>p</sub></sub></td><td><sup>u<sub>p</sub>L<sub>p</sub></sup>/<sub>D<sub>p</sub></sub></td></tr></tbody></table>

b) Tissue Homogeneity (TH): As the name suggests, the TH model assumes a homogeneous extravascular (or tissue) phase (H.12). Integrating the second line over the extravascular space then applying H.12, (18) becomes

$$ \phi_p \frac{\partial c_p}{\partial t} = PS(c_e(t) - c_p(x, t)) - F_p L_p \frac{\partial c_p}{\partial x}. \quad (19) $$

$$ \phi_e \frac{dc_e}{dt} = PS(\bar{c}_p(t) - c_e(t)) $$

This is the TH model, where $\bar{c}_p(t) = \frac{1}{L} \int_0^L c_p(x, t) dx$ ($c_e$ no longer depends on $x$).

c) Two-compartments exchange model (2CXM): Historically ([11], [15]) and similarly to the one-phase case, the main assumption for compartment models is the well-mixed assumption H.11. First, we integrate (18) over the length of the compartment L:

$$ \phi_p \frac{d\bar{c}_p}{dt} = PS(\bar{c}_e(t) - \bar{c}_p(t)) + F_p (c_p(0, t) - c_p(L, t)). $$

Assuming a well-mixed plasma compartment (H.11) then leads to the two-compartment exchange model (2CXM):

$$ \begin{aligned} \phi_p \frac{d\bar{c}_p}{dt} &= PS(\bar{c}_e - \bar{c}_p) + \frac{F_p}{\lambda} (c_{in} - \bar{c}_p) \\ \phi_e \frac{d\bar{c}_e}{dt} &= PS(\bar{c}_p - \bar{c}_e) \end{aligned} \quad (20) $$

d) Adiabatic approximation to the Tissue Homogeneity (AATH): The AATH model assumes that all exchanges with the extravascular volume are made at the end of the plasma volume (mathematically defined within H.13). It is usually perceived as a replacement for the TH model, as it offers an analytical solution [12]. It is defined as:

$$ \begin{aligned} \phi_p \frac{\partial c_p}{\partial t}(x, t) &= -F_p L_p \frac{\partial c_p}{\partial x}(x, t) \\ \phi_e \frac{\partial c_e}{\partial t}(t) &= K^{\text{trans}}(c_p(L, t) - c_e(t)) \end{aligned} \quad (21) $$

Here we derive the AATH model from the initial transport model (6), assuming the adiabatic approximation (H.13). To do so, one follows the same steps as for (8) but replacing H.8 by H.13. These hypotheses allow rewriting (8) as:

$$ \phi_p \frac{\partial c_p}{\partial t}(x, t) = PSL(c_e - c_p)\delta(L - x) - F_p L \frac{\partial c_p}{\partial x}, \quad (22) $$

$$ \phi_e \frac{\partial c_e}{\partial t}(x, t) = PSL(c_p - c_e)\delta(L - x). \quad (23) $$

With $\delta$ the Dirac delta function. Defining $K^{\text{trans}} = \frac{PS \cdot F_p}{PS + F_p}$, one can translate (22, 23) into (21), recovering the well-known AATH model. The derivation is given in appendix C.

2) Biomarkers and measurables: The presented derivations introduce a variety of known TK parameters, for which we gave the mathematical relation to the physical quantities (see table I for a summary). To compare them with one another, we chose to represent them with the following characteristic times: $\tau_a = \frac{L}{u_p}$ which characterizes advection, $\tau_{d,p} = \frac{L^2}{D_p}$ and $\tau_{d,e} = \frac{L^2}{D_e}$ which characterize diffusion respectively in the plasma and the extravascular space, and $\tau_{PS} = (PS)^{-1}$ which characterizes exchanges between the two phases. In section II-B.2 we defined the mean transit time for one phase (15). Based on this and [5], we define $\tau_{t,2} = \tau_{t,1}(\phi_p + \phi_e)$ as the mean transit time for two phases, which is set as the characteristic time for transport since it takes into account all transport processes here represented. The mean transit time, as defined in (14), cannot be recovered analytically in the general case. Non-dimensionizing (8) by defining $\tau = \tau_{t,2}$ and keeping the notations of II-B.1 leads to the system:

$$ \left\{ \begin{aligned} \phi_p \frac{\partial c_p}{\partial \theta} &= \frac{\tau}{\tau_{PS}} (c_e - c_p) - \frac{\tau}{\tau_a} \phi_p \frac{\partial c_p}{\partial \chi} + \frac{\tau}{\tau_{d,p}} \phi_p \frac{\partial^2 c_p}{\partial \chi^2} \\ \phi_e \frac{\partial c_e}{\partial \theta} &= \frac{\tau}{\tau_{PS}} (c_p - c_e) + \frac{\tau}{\tau_{d,e}} \phi_e \frac{\partial^2 c_e}{\partial \chi^2} \\ c_p(0, \theta) &= c_{in}(\tau \theta); & \partial_{\chi} c_p(1, \theta) &= 0 \\ \partial_{\chi} c_e(0, \theta) &= 0; & \partial_{\chi} c_e(1, \theta) &= 0 \end{aligned} \right. . \quad (24) $$

To assess the physical validity of the reduced models – i.e. to what extent the models’ hypotheses are coherent with the physical setting, – we choose a metric that contains information on the derivation hypotheses. For a reduced model $i$ ($i \in \{1C, TH, \dots\}$), we evaluate its residual $\mathcal{R}_i$ on the solution of the full-order model $c_f$. By definition, the residual is the LHS of the model equation minus the RHS. A demonstration example is provided in appendix D, giving

$$ \mathcal{R}_{1C}(c_f) = \frac{\tau}{\lambda \tau_a} ((\lambda - 1)c_{in} - \lambda c_f(L) + \bar{c}_f) + \frac{\tau}{\tau_d} \partial_{\chi} c_f(0) \quad (25) $$

in the one-compartment case. The residual highlights the key hypotheses for (9), namely H.9 (highlighted in the right-most term) and H.11 (highlighted in the left-most term), leading to $\mathcal{R}_{1C}(c_f) = 0$ whenever the physical system satisfies H.9 and H.11. We then apply this method on the residuals of (19), (20), and (21). In the general case, as the residual of the reduced model includes the terms of the partial differential system, and is evaluated on the solution of the full model, it directly embeds the derivation hypotheses. This allows for a more targeted interpretation, compared to a more naive $\mathcal{L}^2$ comparison of the solutions ($\|c_f - c_{reduced}\|_2^2$).

Finally, reference [5] suggests a value of $\lambda = 0.5$ for the 2CXM. We challenge this value and compute an optimal $\lambda_{opt}$ that reduces the difference between $\bar{c}$ and $\lambda c_{out} + (1 - \lambda)c_{in}$. This is equivalent to minimizing the integral over time of the square difference between the two terms ($\mathcal{E}_{\lambda}$). This leads to a time-independent optimal $\lambda_{opt}$,

$$ \mathcal{E}_{\lambda}^2 = \frac{\int_{0}^{+\infty} (\bar{c} - \lambda c_{out} - (1 - \lambda)c_{in})^2 dt}{\int_{0}^{+\infty} \bar{c}^2 dt}; \lambda_{opt} = \arg\min_{\lambda} \mathcal{E}_{\lambda} \quad (26) $$
The 2CXM with $\lambda = \lambda_{opt}$ is noted 2CXM* in the following.

3) Asymptotic behaviors: As the governing equations cannot be solved analytically in the general case, we study the different regimes through simulation. System (24) is solved through an in-house finite differences scheme with an implicit Euler scheme over time. As the equation is linear, this results in the solution of an equivalent linear system. The Péclet number and the ratio $\tau_{d,e}/\tau$ are first set to a high value: negligible diffusion. They are gradually reduced to observe the impact of diffusion in the plasma and the extravascular space on the performance of the reduced models (in terms of the residual error). The residuals of the reduced models (19), (20), and (21) are evaluated with the solution of the FOM (24). To get a relative error value, the $L^2$ norm of the residuals is then normalized by dividing by the integral over time of the squared input signal. This enables us to quantify to what extent the hypotheses are valid, and how good the reduced models are at approximating the full model, depending on the time characteristics $\tau_a, \tau_{d,p}, \tau_{d,e}$ and $\tau_{PS}$.

### III. SIMULATION RESULTS

In this section, we first present the parameter study of the TH, 2CXM, and AATH models (section III-A). Based on this, we have developed an algorithm that determines the maximum length of a compartment that ensures validity of the models' hypotheses. The algorithm is then expressed in pseudocode (section III-B) before presenting a few examples applied to molecules widely adopted in the medical imaging community (section III-C).

#### A. The relevance of diffusion

We here study the effect of having a diffusion term in the plasma and extravascular compartments ($\tau_{d,p}, \tau_{d,e} < +\infty$) (figure 4). The residual error is calculated on a grid of 50 × 50 parameter values for the couple ($\tau_{d,e}/\tau, P_e$), both ranging logarithmically between $10^{-4}$ and $10^4$, and for several values of PS. The compartment length is set at $L = 1$ mm, which is the typical spatial resolution for imaging the heart and other organs [28]. Blood velocity is chosen as the velocity in capillaries, as it is the place where diffusion is supposedly most impactful. A standard value is $u = 500$ µm/s [29]. A realistic input function is taken as the inlet concentration profile. Its mathematical formula is given by (27), and its shape and parameters are given in figure 3.

$$c_{in}(t) = \alpha e^{-\frac{(t-t_0)^2}{2\sigma_0^2}} + \beta \frac{e^{-\frac{t}{t_1}}}{1 + e^{-\frac{t-t_2}{\sigma_1}}} \quad (27)$$

Time and space steps $n = 5000$ were required for appropriate convergence of the numerical solver. The blood volume fraction is set the same in all simulations and chosen at a typical value $\phi_p = 0.04$ [22]. The maps for residual errors, taking PS values ranging from 0.01/s to 100/s, are shown in figure 4.

For small $P_e \lesssim 1$ (corresponding to a diffusion-dominated transport) all three models approximate the FOM well, as diffusion in the FOM homogenizes the concentration of transported molecules in the ROI.

<figure><img src="image_4.png" /><figcaption>Fig. 3: Shape of the signal given as input to the system. Inset: values of the parameters in (27).</figcaption></figure>

With increasing $\tau_{d,e}/\tau$ (corresponding to a smaller $D_e$ for a given $L$), extravascular inhomogeneities increase, resulting in an increasing residual error, especially for the TH models. Increasing $P_e$ (decreasing $D_p$ for a given $L$ and $u_p$) has a similar effect on the plasma phase, explaining the reduced performance of the 2CXM*. If $\tau_{d,e}/\tau$ and $P_e$ are increased, the effect is amplified. The larger the PS, the greater the effect, as more molecules leave the plasma and enter the extravascular space, increasing heterogeneity in both phases (ultimately enlarging the white regions of higher errors).

The AATH performs worst here as the exchange with the extravascular space is assumed to occur only at the end of the plasma volume, so significant residual error accumulates over the length $[0, L]$ in the ROI. Although it is a less reduced model in the sense of still considering the spatial dimension, TH performs worse than 2CXM, especially at high $PS \gg 1/s$, since the latter is compatible with certain inhomogeneities. At low $PS \ll 1/s$ (figure 4a), the residuals of the three reduced models do not exceed 5 % of the input concentration magnitude for all values of the other parameters.

These observations are summarized in table II, where only the extreme regimes are reported.

#### B. Definition of the maximum ROI size

Figure 4 shows an increase in the residual error with $P_e$ and $\tau_{d,e}/\tau$. As $P_e$ and $\tau_{d,e}/\tau$ both depend on $L$, this justifies the existence of a maximum ROI size $L_{max}$ above which the reduced models' residuals introduce more error than the desired value. Above this threshold, the physical hypotheses are assumed to be violated (refer to the paragraph about the residual calculation around (25)).

For a given molecule (transported inside a blood vessel with specified blood velocity and permeability PS) the maximum ROI size is estimated according to the algorithm III.1.

**Algorithm III.1: Maximum ROI size**

INPUT: Molecule, Blood velocity, Permeability, Blood volume fraction, Temperature

OUTPUT: $L_{max}$

START:
<figure><img src="image_5.png" /><figcaption>Fig. 4: Residual maps of the three reduced two-phase models, namely 2CXM*, TH, and AATH, under multiple values of PS. The $\tau_{d,e}/\tau$ ratio is represented on the horizontal axis; $P_e$ is represented on the vertical axis. (a) PS = 0.01/s. (b) PS = 1/s. (c) PS = 100/s. The dashed line represents the level line at 5% relative residual error. The compartment length is fixed at L = 1 mm and velocity at u = 0.5 mm/s.</figcaption></figure>

```
INITIALIZE L ← 1
INITIALIZE R ← 0
DO Monte Carlo Volume estimation
COMPUTE Molecular radius
COMPUTE Molecular sphericity
COMPUTE Stokes-Einstein Diffusivity
WHILE R < ε DO
    SET c_f ← solution of (24)
    SET R ← R_i(c_f)
    SET L ← L · 10^0.02
    SET L_max ← L · 10^-0.02
END.
```

The molecular diffusivity coefficient is computed via the Stokes-Einstein-Sutherland (SES) formula $D = \frac{k_B T}{6\pi\mu r}$, where $k_B$ is the Boltzmann constant, $T$ is the temperature, $\mu$ is the viscosity of the considered solvent and $r$ is the equivalent hydrodynamic radius of the tracing molecule. The latter refers to the radius of the molecule if it were spherical. In reality, no molecule is purely spherical, so an estimate of the sphericity of the molecule indicates to what extent one can trust the calculated diffusivity coefficient. The molecular radius and

<table><caption>TABLE II: Model validity (color-coded by green when the error is below 0.05 – meaning 5 % – and red otherwise) and relative residual error (numerical values) for low ($P_e, \tau_{d,e}/\tau = 10^{-4}$; PS = $10^{-2}$/s) and high ($P_e, \tau_{d,e}/\tau = 10^4$; PS = $10^2$/s) values of the parameters.</caption><thead><tr><th rowspan="2">Model</th><th rowspan="2"></th><th colspan="2">Low PS</th><th colspan="2">High PS</th></tr><tr><th>Low τ<sub>d,e</sub></th><th>High τ<sub>d,e</sub></th><th>Low τ<sub>d,e</sub></th><th>High τ<sub>d,e</sub></th></tr></thead><tbody><tr><td rowspan="2">2CXM*</td><td>Low P<sub>e</sub></td><td>2 · 10<sup>-5</sup></td><td>5 · 10<sup>-7</sup></td><td>9 · 10<sup>-6</sup></td><td>2 · 10<sup>-6</sup></td></tr><tr><td>High P<sub>e</sub></td><td>2 · 10<sup>-3</sup></td><td>5 · 10<sup>-3</sup></td><td>4 · 10<sup>-3</sup></td><td>2 · 10<sup>-1</sup></td></tr><tr><td rowspan="2">TH</td><td>Low P<sub>e</sub></td><td>9 · 10<sup>-16</sup></td><td>4 · 10<sup>-15</sup></td><td>8 · 10<sup>-13</sup></td><td>2 · 10<sup>-7</sup></td></tr><tr><td>High P<sub>e</sub></td><td>3 · 10<sup>-10</sup></td><td>8 · 10<sup>-4</sup></td><td>5 · 10<sup>-5</sup></td><td>5 · 10<sup>1</sup></td></tr><tr><td rowspan="2">AATH</td><td>Low P<sub>e</sub></td><td>2 · 10<sup>-5</sup></td><td>5 · 10<sup>-7</sup></td><td>1 · 10<sup>-5</sup></td><td>3 · 10<sup>-6</sup></td></tr><tr><td>High P<sub>e</sub></td><td>9 · 10<sup>-3</sup></td><td>9 · 10<sup>-3</sup></td><td>1 · 10<sup>0</sup></td><td>5 · 10<sup>1</sup></td></tr></tbody></table>

<table><caption>TABLE III: Estimation of diffusion coefficients with a Monte-Carlo (MC) based approach, compared to the TDA technique [30]. Values are given in µm²/s. The temperature and viscosity are taken as those of the experimental setting.</caption><thead><tr><th>Molecule</th><th>Adapted MC method</th><th>Literature values</th></tr></thead><tbody><tr><td>Omniscan®</td><td>537</td><td>501</td></tr><tr><td>Dotarem®</td><td>539</td><td>546</td></tr><tr><td>Elucirem®</td><td>440</td><td>332</td></tr></tbody></table>

sphericity are estimated thanks to a Monte Carlo volume estimation method. This method provides similar results compared to experimental measurements on Gadolinium-based contrast agents [30] (see table III for several examples).

To eliminate the dependency on extravascular diffusion ($\tau_{d,e}$), we assume that the tracer molecules diffuse normally in the extravascular space. According to the SES formula, $P_e$ and $\tau_{d,e}/\tau$ are thus related through:

$$ \frac{\tau_{d,e}}{\tau} = f(P_e) = \frac{P_e^2}{e - P_e - (1 - P_e) \frac{\mu_e}{\mu_p} \phi_p + \phi_e} \quad (28) $$

$\mu_p$ (resp. $\mu_e$) is the blood (resp. extravascular space) viscosity. Therefore, we reduce the exploration of parameters to the graph of the function $\tau_{d,e}/\tau = f(P_e)$ and more specifically to the variation of $L$ only (see definitions in II-B.2). The length is initialized at a minimum value $L = 1$ µm, then increased iteratively – multiplying the previous value by $10^{\frac{1}{50}} \approx 1.047$, multiplying by 10 after 50 steps – until the error exceeds the desired precision. Algorithm III.1 has been implemented in PYTHON and is available online<fnref n="1" />.

### C. Orders of magnitude

In this section, the threshold for the residual error is chosen at 0.03 (i.e., 3 %, placing the threshold line in the dark blue regions in figure 4) as an example. For $e$ = interstitium space (essentially water) and $p$ = blood plasma, the ratio $\frac{\mu_e}{\mu_p} \approx 0.16$ at $T = 37$ °C [31].

<footnote n="1">¹https://gitlab.com/imaging-flow-transport-inria/diffusion-length-scale</footnote>
<figure><img src="image_6.png" /><figcaption>Fig. 5: Maximum ROI length that ensures the residual error for the reduced models is < 0.03, as a function of the PS product (the region below the curve gives errors lower than 3 %, the region above, higher than 3 %). The 2CXM* corresponds to the 2CXM with the optimized λ<sub>opt</sub> factor (see (26) and section III-D). For 2CXM, λ is fixed at λ = 0.5.</figcaption></figure>

1) Evolution with PS: The maximum length is computed for increasing values of PS. Diffusivity in blood is set at the average value of D<sub>p</sub> = 110 µm<sup>2</sup>/s (table IV). The results are displayed in figure 5 and discussed here.

a) 2CXM: The variation of the compartment length with PS is not monotonous: it starts around 4.4 mm for low PS values, reaches a minimum value of 48.7 µm for PS around 10/s, and saturates at 200 µm at high PS.

b) 2CXM*: Unlike the pure 2CXM, the maximum ROI length here always decreases as PS increases. It is always higher here than for the simple 2CXM demonstrating that λ = 0.5 is not always the optimal choice, especially at relatively high PS ≳ 0.1/s (see section III-D). For the studied range of PS, the maximum length drops by a factor of around 20.

c) TH: The maximum validity length strictly decreases as the vessel permeability increases, reaching higher values than for the 2CXM at low permeability and lower values at higher permeability, which is coherent with the results in III-A.

d) AATH: The variation of the maximum length resembles the one of the TH model but remains always lower at a given PS (coherent with section III-A). Nevertheless, at high PS, the maximum lengths are similar for the AATH and TH model, highlighting a physical similarity in the two models despite their different derivation hypotheses (see section II-B.1).

2) Application to imaging molecules: The maximum ROI size for each model is computed for four tracer molecules, spanning different imaging techniques. Their hydrodynamic radius is estimated based on their three-dimensional configuration available online<fnref n="2" />. The proposed algorithm is applied to each molecule, and the results are gathered in table IV.

The maximum length of a 2CXM is around 4 mm for all molecules and is larger by up to 1 mm in the 2CXM* case. It is larger than 7 cm for a TH model and varies from 1 cm to higher than 10 cm for the AATH. Overall, it remains larger than the millimeter (corresponding to a standard imaging ROI

<figure><img src="image_7.png" /><figcaption>Fig. 6: (a) Evolution of the λ<sub>opt</sub> factor, (b) ε<sub>λ<sub>opt</sub></sub>, the root relative mean square difference between c̄ and (1 − λ<sub>opt</sub>)c<sub>in</sub> + λ<sub>opt</sub>c<sub>out</sub>, (d) ε<sub>0.5</sub>, relative error introduced by setting λ = 0.5 as in [5], and (c) ε<sub>add</sub> = ε<sub>0.5</sub> − ε<sub>λ<sub>opt</sub></sub>.</figcaption></figure>

size) for all molecules and the three models studied, even reaching values larger than 10 cm for the TH model. The maximum length also appears to increase with increasing molecular diffusivity, while it decreases with increasing vessel permeability, as it was hinted at in sections III-A and III-C.1.

D. About the well-mixed hypothesis

Finally, we here study the origin of the difference between the 2CXM and 2CXM* performances at high PS seen on figure 5. For that, let us compute the factor λ<sub>opt</sub>, introduced in H.11 and recalled in (29), for different values of the Péclet number and the permeability-surface area product.

$$
\bar{c} = (1 - \lambda)c_{\text{in}} + \lambda c_{\text{out}} \quad (29)
$$

The optimal interpolation factor is computed through (26) – based on the solution of the FOM (24) – and presented in figure 6a along with the error made by doing this approximation, relative to the effective average concentration $\bar{c}$ (on figure 6b). For low PS and $P_e$, the computed value is $\lambda_{opt} \approx 0.66$, which aligns with the $\frac{2}{3}$ value derived in section II-B.3.b. Furthermore, we observe a variation from 0.5 to 1 for increasing PS/F, as it is argued in [5]. Finally, the constant value $\lambda = 0.5$ introduces an additional relative error (up to 16% at medium permeability values) compared to the case when $\lambda = \lambda_{opt}$ (figure 6c and 6d). This may explain the non-monotonous behavior of $P_e$ vs. PS for the 2CXM in figure 5.

IV. DISCUSSIONS

The analytical derivation enlightened several physical hy-
potheses, which validity has been challenged analytically and
numerically. Section IV-A discusses the main results while
putting them into perspective with literature-derived imaging
data. Section IV-B then recapitulates the main contributions of
this work while defining their limits. Finally, we summarize

<footnote n="2">²https://pubchem.ncbi.nlm.nih.gov, https://www.rcsb.org/</footnote>
<table><caption>TABLE IV: Physical properties of several molecules injected as tracers for different imaging techniques. The last four columns represent the maximum size of an ROI that ensures validity of the models' hypotheses, depending on the molecule. The PS column includes the reference for the values of the permeability-surface area.</caption><thead><tr><th>Molecule</th><th>Imaging modality</th><th>Molar mass g/mol</th><th>Radius Å</th><th>D<sub>blood</sub> µm<sup>2</sup>/s</th><th>PS 10<sup>-3</sup>/s</th><th>2CXM mm</th><th>2CXM* mm</th><th>TH mm</th><th>AATH mm</th></tr></thead><tbody><tr><td>Gadopentetic acid</td><td>Magnetic resonance</td><td>547.57</td><td>4.46</td><td>114.41</td><td>0.20 [32]</td><td>4.34</td><td>4.34</td><td>&gt; 100</td><td>&gt; 100</td></tr><tr><td>Iopamidol</td><td>Computer Tomography</td><td>777.09</td><td>4.63</td><td>110.31</td><td>2.32 [33]</td><td>3.87</td><td>4.60</td><td>&gt; 100</td><td>11.22</td></tr><tr><td>Indocyanine Green</td><td>ICG Angiography</td><td>774.96</td><td>5.48</td><td>93.12</td><td>0.37 [34]</td><td>4.34</td><td>4.47</td><td>&gt; 100</td><td>59.57</td></tr><tr><td>Cholyl-lysyl-fluorescein</td><td>Intravital microscopy</td><td>926.12</td><td>5.76</td><td>88.63</td><td>2.64 [35]</td><td>3.65</td><td>4.60</td><td>77.18</td><td>10.00</td></tr></tbody></table>

the main conclusions of this work and introduce possible future explorations (section IV-C).

## A. Physical validation in different regimes

1) **Diffusion:** The analytical study in section II-B highlights the two extreme regimes of the advection-diffusion equation: the time-translation at low diffusion (paragraph II-B.3.a) and the compartment model at high diffusion (paragraph II-B.3.b). This result is further generalized through the numerical analysis of section III-A: diffusion processes tend to reduce concentration gradients, leading to well-mixed or homogeneous phases. The TK models appear to be coherent with the more physical equation (24) at highly diffusive regimes.

2) **Permeability:** Conversely, high vascular permeability weakens the physical coherence of reduced TK models, as the maximum ROI size to reach model validity decreases with increasing permeability. High vessel permeability accelerates tracer exchange between plasma and extravascular phases compared to low permeability. This increases concentration heterogeneity – especially in the extravascular phase. This explains why the AATH and TH models’ assumptions (tissue homogeneity) perform poorly at high permeability. However, AATH requires smaller ROIs than TH at low permeability, as its derivation relies on stronger physical assumptions.

3) **Literature-derived imaging data:** Thanks to the presented framework, it is possible to examine the choice of model that has been made to study imaging data in the literature. Let us take as an example the data analyzed by [6]: the voxel (ROI) size is $L = 1.25$ mm with a measured $PS = 2.17 \cdot 10^{-3}$ s. According to our present work (figure 5), the ROI size is small enough for the 2CXM to be coherent with the physical model. Indeed, in the article, the deconvolution displays a bi-exponential kernel, which is coherent with a 2CXM model (see [6], [14] for more explanation on the bi-exponential kernel).

However, as the ground truth value for the TK parameters cannot be measured directly (as stressed out in [6], [7]), we here validate the authors model choice based on the parameters estimated by that same model: it is a vicious circle.

To go beyond this circular validation process, we analyze the parameter estimation capacity of the 2CXM on two synthetic cases S.1 and S.2, and two literature-derived datasets [6], [7]. The synthetic cases establish ground truth by simulating the full model (24) with chosen parameters reported in the first two lines of table V, $\tau_{d,p} = 13.7 \cdot 10^3$ s (representing Gd-DTPA molecular diffusion in a voxel of size $L = 1.25$ mm), $\lambda = 0.5$ and $c_{in}$ interpolated from the data in figure 4a of [6]. The literature data is represented as black dots in

<figure><img src="image_8.png" /><figcaption>Fig. 7: Comparison of the estimated parameters with the reference. The value estimated through a 2CXM inverse problem is divided by the literature or synthetic reference value (see table V). For example, a value of 1 means that the estimated parameter is equal to the reference one.</figcaption></figure>

<table><caption>TABLE V: Chosen reference parameters for the two synthetic cases S.1 and S.2. Reported literature parameters for the two examples: [6], [7].</caption><thead><tr><th>Quantity (unit)</th><th>F<sub>p</sub>/λ (s<sup>-1</sup>)</th><th>PS (s<sup>-1</sup>)</th><th>φ<sub>p</sub> (-)</th><th>φ<sub>e</sub> (-)</th></tr></thead><tbody><tr><td>S.1</td><td>0.4</td><td>0.56</td><td>0.04</td><td>0.65</td></tr><tr><td>S.2</td><td>0.4</td><td>0.0056</td><td>0.04</td><td>0.65</td></tr><tr><td>[6]</td><td>0.0075</td><td>0.0022</td><td>0.20</td><td>0.20</td></tr><tr><td>[7]</td><td>0.020</td><td>0.0056</td><td>0.08</td><td>0.40</td></tr></tbody></table>

figure 8 and gathered from figure 4b of [6] and figure 1b of [7]. The 2CXM parameters ($F_p/\lambda$, $PS$, $\phi_p$, and $\phi_e$) estimated from these curves and reported by the respective authors are gathered in the last two lines of table V. For all four datasets, we estimated 2CXM parameters via the CMA-ES optimization method [36], [37]. Figure 7 shows the relative deviations of estimated parameters from reference values across all test cases. In figure 8, we compared the raw literature data (black dots) with several simulations: 1) (20) with reference parameters (dashed black lines), 2) (20) with 2CXM-estimated parameters (solid red lines), and 3) (24) with 2CXM-estimated parameters, $\lambda = 0.5$, and $\tau_{d,p} = 13.7 \cdot 10^3$ s (solid blue lines).

In the synthetic case S.1, $F_p/\lambda$ and $PS$ are well estimated, $\phi_p$ is overestimated and $\phi_e$ is underestimated. In S.2, $F_p/\lambda$, $PS$, $\phi_p$, and $\phi_e$ are better estimated, even though $\phi_p$ is still underestimated (this is solved by fixing the $\phi_e$ to $\phi_p$ relation, which suggests identifiability issues). With data from [6], $F_p/\lambda$ is well estimated, yet all other parameters differ drastically from the reference value. Yet, as displayed on figure 8a in dashed lines, the 2CXM model evaluated with the reference
<figure><img src="image_9.png" /><figcaption>Fig. 8: Concentration-time curves obtained from the DCE-MRI literature. $C_T = \phi_p c_p + \phi_e c_e$ is the total tissue concentration. The black dashed lines are evaluations of the 2CXM (20) with values presented in the reference. The red solid lines are evaluations of the 2CXM with fitted parameters. The blue solid lines are evaluations of the full model (24) with the 2CXM-fitted parameters (with molecular diffusion, i.e., $\tau_{d,p} = 13.7 \cdot 10^3$s).</figcaption></figure>

parameters does not match the data (more than 44 % relative error). This suggests a possible mistake in the reporting of the reference parameters. Moreover, the 2CXM fit estimates PS = 0.34/s, which places the physical system outside the model validity range (this PS would require $L \approx 200$µm, see figure 5). The parameters we fitted appear to be fairly close to the estimates in [7], although our fit seems to match the data more closely (a difference of 6.9 % when evaluating the 2CXM model with the estimates in [7], compared to 2.2 % for our fit). Moreover, the evaluation of the full model with the 2CXM-estimated parameters (solid blue lines in figure 8) fits the data closer for [7] than for [6]. This endorses the physical invalidity of the reduced model in the latter case, consistent with that a fit with the 2CXM model yields values outside of the model's validity range. As mentioned earlier in this article, the full model may nevertheless fit the data when a higher diffusion coefficient compared to molecular diffusion alone is taken into account.

This study highlights a few points: 1) the reduced models (here 2CXM) may lead to misinterpretation of the observed phenomena whenever the "ground truth" parameters are outside the validity range, defined by the error on the models' residual, even in case of a good fit; 2) the *no diffusion* hypothesis (H.9) may not be assumed for compartment modeling whenever the model-inferred parameters are outside of the range within which the model residuals are below a defined threshold; 3) the 2CXM model carries identifiability issues which may lead to different results depending on the optimization method to fit the model to the data. This has already been pointed out in [9], [22], and can be analyzed by the framework presented in [38]. These observations challenge the interpretation of classic TK parameters – estimated through reduced models – and encourage further research in the estimation of ground truth values for model validation.

## B. Main impacts

1) **Contributions:** This work presents, to our knowledge for the first time, the derivation of the tracer-kinetics models from

the physical principles of mass conservation. Figure 2 summarizes the hypotheses and derivations of this study. Derived from PK compartment models, the presented TK models have a priori not been validated against physical transport equations [24], [39]. For example, [14] presents further reductions of the three two-phase models: 2CXM, TH, and AATH. This study offers a new perspective on the physical origins of key TK parameters, presenting mathematical formulas for mean transit time, $F_p$, PS, and $K^{\text{trans}}$ (Table I). Our derivation also retains the diffusion term, potentially revealing additional insights into transport dynamics in literature-derived imaging data, as suggested in section IV-A.3. A refinement of the historical well-mixed hypothesis is proposed for the derivation of compartment models facilitated by numerical simulations. The approach is illustrated through examples involving commonly prescribed contrast agents.

2) **Code availability:** In conjunction with the present manuscript, we introduce a Python software tool¹ that computes the upper limit of a compartment's volume for it to be considered well-mixed, as well as the best interpolation factor $\lambda_{\text{opt}}$ for given concentration profiles $c_{\text{in}}$, $c_{\text{out}}$ and $\bar{c}$.

3) **Limitations:** Table IV displays literature-derived values of the PS for four tracer molecules. These values were obtained through various measurement techniques across different organs, for multiple disease states, including colorectal cancer and healthy liver tissues. Given that the physical validity of TK models' assumptions significantly deteriorates as this parameter increases (see section III-C.1), the specific outcomes for each tracer molecule may vary depending on the organ of interest and the patient's disease state. For instance, hepatic sinusoids, with their multiple fenestrations, exhibit far greater permeability than regular capillaries, significantly reducing the maximum ROI size in liver images compared to other organs. In the brain, blood-brain barrier disruption, as seen in stroke, increases permeability 6 to 8 times [40], dividing the maximum ROI size by approximately 3 for the TH and AATH models. Moreover, as presented in table I, PS appears to depend on $L_p L_{\partial\Sigma}/V_{\text{tot}}$ and thus on ROI size, potentially explaining the voxel size-dependent PS variations reported in [41]. This aspect deserves further investigation.

The current voxel (or ROI) size thresholding approach should be contextualized within the framework of dynamic imaging protocols: the choice of spatial resolution (usually related to the matrix size and field of view) depends on the desired level of details, the organ and region of interest, and on the desired temporal resolution of the dynamic signal [42]. The presented numbers provide a theoretical threshold, which is not intended to be used as-is, but rather serves as a warning regarding the model's interpretability (see section IV-A.3).

In this study, even though the impact of input signal duration is analyzed analytically in the simpler case (section II-B.3), it is not analyzed numerically in the general case; however, a shorter input signal time might reduce model validity. This may be further analyzed within the proposed framework.

## C. Conclusion and future work

The 2CXM presents an optimal balance between model complexity and validity, outperforming the TH and AATH
models regarding the physical interpretability. In particular, its maximum validity size exceeds 250 µm, even under conditions of high permeability, despite no longer having any space dependency. However, conventional dynamic imaging techniques operate at voxel sizes above 500 µm [28], [42]. Our findings suggest that keeping the same accuracy at a lower spatial resolution would require an increased diffusivity coefficient. At high spatial resolutions (small voxel size, $L < 10$ µm), molecular diffusion ensures well-mixed compartments. For larger voxel sizes (lower resolution), our study demonstrates that the well-mixed hypothesis becomes less valid. Assuming well-mixed compartments in such cases increases model residuals, leading to inaccurate estimates of tissue properties from dynamic imaging. Our framework reveals that the current literature data appear to include apparent diffusion, previously undetected by other models. This may offer new interpretations of the data.

With further validation, our algorithm could define a validity range for voxel size based on expected or known ranges of permeability and blood flow. This would help guide imaging protocols or inform the development of next-generation devices. This model-guided approach to voxel size estimation may encourage further research to either recommend larger voxel sizes – saving time and cost – or set targets for reaching smaller voxel sizes through improved protocols or devices.

Further research is also required to understand the apparent diffusion coefficient. Indeed, in a network of interconnected blood vessels, the variability of transport velocities leads to a well-mixed distribution of the substance in the larger encompassing region of interest. This seems to qualitatively originate from a diffusion behavior associated with the network's geometry and hemodynamic properties. The study in [43] offers mathematical insights into the emergence of diffusive behavior from network structures. Some research efforts have been dedicated to linking tree structures to the manifestation of apparent diffusion [44], [45]. An initial attempt to relate this phenomenon to the geometrical characteristics of the tree was undertaken in [46], which, although incomplete, establishes a foundation for extending the interpretation of TK models.

## REFERENCES

[1] L. Kuepfer, C. Niederalt, T. Wendl, et al., "Applied Concepts in PBPK Modeling: How to Build a PBPK/PD Model," *CPT Pharmacometrics Syst Pharmacol*, vol. 5, no. 10, pp. 516-531, Oct. 2016, ISSN: 2163-8306. DOI: 10.1002/psp4.12134.

[2] L. O. Schwen, L. Kuepfer, and T. Preusser, "Modeling approaches for hepatic spatial heterogeneity in pharmacokinetic simulations," en, *Drug Discovery Today: Disease Models, Computational Models of Liver Disease* 2016, vol. 22, pp. 35-43, Dec. 2016, ISSN: 1740-6757. DOI: 10.1016/j.ddmod.2017.09.002.

[3] W. Albrecht, F. Kappenberg, T. Brecklinghaus, et al., "Prediction of human drug-induced liver injury (DILI) in relation to oral doses and blood concentrations," en, *Arch Toxicol*, vol. 93, no. 6, pp. 1609-1637, Jun. 2019, ISSN: 0340-5761, 1432-0738. DOI: 10.1007/s00204-019-02492-9.

[4] J. Dichamp, G. Cellière, A. Ghallab, et al., "In vitro to in vivo acetaminophen hepatotoxicity extrapolation using classical schemes, pharmacodynamic models and a multiscale spatial-temporal liver twin," *Frontiers in Bioengineering and Biotechnology*, vol. 11, 2023, ISSN: 2296-4185.

[5] G. Brix, S. Zwick, F. Kiessling, and J. Griebel, "Pharmacokinetic analysis of tissue microcirculation using nested models: Multimodel inference and parameter identifiability: Analysis of tissue microcirculation using nested models," en, *Med. Phys.*, vol. 36, no. 7, pp. 2923-2933, Jun. 2009, ISSN: 00942405. DOI: 10.1118/1.3147145.

[6] G. Brix, J. Griebel, F. Kiessling, and F. Wenz, "Tracer kinetic modelling of tumour angiogenesis based on dynamic contrast-enhanced CT and MRI measurements," en, *Eur J Nucl Med Mol Imaging*, vol. 37, no. S1, pp. 30-51, Aug. 2010, ISSN: 1619-7070, 1619-7089. DOI: 10.1007/s00259-010-1448-7.

[7] D. L. Buckley, "Uncertainty in the analysis of tracer kinetics using dynamic contrast-enhanced T1-weighted MRI," en, *Magn. Reson. Med.*, vol. 47, no. 3, pp. 601-606, Mar. 2002, ISSN: 0740-3194, 1522-2594. DOI: 10.1002/mrm.10080.

[8] W. H. Sommer, S. Sourbron, A. Huppertz, M. Ingrisch, M. F. Reiser, and C. J. Zech, "Contrast agents as a biological marker in magnetic resonance imaging of the liver: Conventional and new approaches," en, *Abdom Imaging*, vol. 37, no. 2, pp. 164-179, Apr. 2012, ISSN: 0942-8925, 1432-0509. DOI: 10.1007/s00261-011-9734-9.

[9] C. Audebert and I. E. Vignon-Clementel, "Model and methods to assess hepatic function from indocyanine green fluorescence dynamical measurements of liver tissue," en, *European Journal of Pharmaceutical Sciences*, vol. 115, pp. 304-319, Mar. 2018, ISSN: 09280987. DOI: 10.1016/j.ejps.2018.01.008.

[10] H.-S. Lim, "Evolving role of modeling and simulation in drug development," *Transl Clin Pharmacol*, vol. 27, no. 1, pp. 19-23, Mar. 2019, Publisher: Korean Society for Clinical Pharmacology and Therapeutics. DOI: 10.12793/tcp.2019.27.1.19.

[11] M. F. Morales and R. E. Smith, "On the theory of blood-tissue exchange of inert gases: VI. Validity of approximate uptake expressions," en, *Bulletin of Mathematical Biophysics*, vol. 10, no. 3, pp. 191-200, Sep. 1948, ISSN: 0007-4985, 1522-9602. DOI: 10.1007/BF02477493.

[12] M. Ingrisch and S. Sourbron, "Tracer-kinetic modeling of dynamic contrast-enhanced MRI and CT: A primer," en, *J Pharmacokinet Pharmacodyn*, vol. 40, no. 3, pp. 281-300, Jun. 2013, ISSN: 1567-567X, 1573-8744. DOI: 10.1007/s10928-013-9315-3.

[13] S. P. Sourbron and D. L. Buckley, "On the scope and interpretation of the Tofts models for DCE-MRI: On the Scope and Interpretation of the Tofts models," en, *Magn. Reson. Med.*, vol. 66, no. 3, pp. 735-745, Sep. 2011, ISSN: 07403194. DOI: 10.1002/mrm.22861.

[14] S. P. Sourbron and D. L. Buckley, "Classic models for dynamic contrast-enhanced MRI: CLASSIC MODELS FOR DCE-MRI," en, *NMR Biomed.*, vol. 26, no. 8, pp. 1004-1027, Aug. 2013, ISSN: 09523480. DOI: 10.1002/nbm.2940.

[15] G. Brix, M. L. Bahner, U. Hoffmann, A. Horvath, and W. Schreiber, "Regional Blood Flow, Capillary Permeability, and Compartmental Volumes: Measurement with Dynamic CT—Initial Experience," en, *Radiology*, vol. 210, no. 1, pp. 269-276, Jan. 1999, ISSN: 0033-8419, 1527-1315. DOI: 10.1148/radiology.210.1.r99ja46269.

[16] C. Ancey, P. Bohorquez, and J. Heyman, "Stochastic interpretation of the advection-diffusion equation and its relevance to bed load transport," en, *Journal of Geophysical Research: Earth Surface*, vol. 120, no. 12, pp. 2529-2551, 2015, _eprint: https://onlinelibrary.wiley.com/doi/pdf/10.1002/2014JF003421*, ISSN: 2169-9011. DOI: 10.1002/2014JF003421.

[17] R. Coosemans, W. Dekeyser, and M. Baelmans, "Mean-field transport equations and energy theorem for plasma edge turbulent transport," en, *Journal of Plasma Physics*, vol. 90, no. 2, p. 905 900 202, Apr. 2024, ISSN: 0022-3778, 1469-7807. DOI: 10.1017/S0022377824000163.

[18] J. Du, D. Kim, G. Alhawael, D. N. Ku, and A. L. Fogelson, "Clot Permeability, Agonist Transport, and Platelet Binding Kinetics in Arterial Thrombosis," en, *Biophysical Journal*, vol. 119, no. 10, pp. 2102-2115, Nov. 2020, ISSN: 00063495. DOI: 10.1016/j.bpj.2020.08.041.

[19] J. M. Oakes, S. C. Shadden, C. Grandmont, and I. E. Vignon-Clementel, "Aerosol transport throughout inspiration and expiration in the pulmonary airways," en, *Numer Methods Biomed Eng*, vol. 33, no. 9, e2847, Sep. 2017, ISSN: 2040-7939, 2040-7947. DOI: 10.1002/cnm.2847.
[20] S. Lynch, N. Nama, and C. A. Figueroa, "Effects of non-Newtonian viscosity on arterial and venous flow and transport," English, *Sci Rep*, vol. 12, no. 1, Nov. 2022, Place: Berlin Publisher: Nature Portfolio WOS:000955946400067, ISSN: 2045-2322. DOI: 10.1038/s41598-022-19867-1.

[21] L. Possenti, A. Gallo, P. Vitullo, et al., "A Computational Model of the Tumor Microenvironment Applied to Fractionated Radiotherapy," en, in *Quantitative Approaches to Microcirculation*, A. Linninger, K.-A. Mardal, and P. Zunino, Eds., vol. 36, Series Title: SEMA SIMAI Springer Series, Cham: Springer Nature Switzerland, 2024, pp. 23–47, ISBN: 978-3-031-58518-0 978-3-031-58519-7. DOI: 10.1007/978-3-031-58519-7_2.

[22] I. E. Vignon-Clementel, N. Jagiella, J. Dichamp, et al., "A proof-of-concept pipeline to guide evaluation of tumor tissue perfusion by dynamic contrast-agent imaging: Direct simulation and inverse tracer-kinetic procedures," *Frontiers in Bioinformatics*, vol. 3, 2023, ISSN: 2673-7647.

[23] J. Kowalski, D. Drasdo, P. Kottman, J. Dichamp, and I. Vignon-Clementel, "Microvascular modeling for medical imaging and toxicity assessment," in *Quantitative approaches to microcirculation: mathematical models, computational methods, measurements and data analysis*, ser. SEMA SIMAI Springer Series, Springer, 2023.

[24] S. Sourbron, "A Tracer-Kinetic Field Theory for Medical Imaging," en, *IEEE Trans. Med. Imaging*, vol. 33, no. 4, pp. 935–946, Apr. 2014, ISSN: 0278-0062, 1558-254X. DOI: 10.1109/TMI.2014.2300450.

[25] M. Mescam, M. Kretowski, and J. Bezy-Wendling, "Multiscale Model of Liver DCE-MRI Towards a Better Understanding of Tumor Complexity," en, *IEEE Trans. Med. Imaging*, vol. 29, no. 3, pp. 699–707, Mar. 2010, ISSN: 0278-0062, 1558-254X. DOI: 10.1109/TMI.2009.2031435.

[26] G. K. Batchelor, *An Introduction to Fluid Dynamics*, en, Cambridge University Press, 1967, Google-Books-ID: Rla7OihRvUgC, ISBN: 978-0-521-66396-0.

[27] J. B. Bassingthwaighte, C. A. Goresky, and J. H. Linehan, "Modeling in the Analysis of the Processes of Uptake and Metabolism in the Whole Organ," en, in *Whole Organ Approaches to Cellular Metabolism*, J. B. Bassingthwaighte, J. H. Linehan, and C. A. Goresky, Eds., New York, NY: Springer New York, 1998, pp. 3–27, ISBN: 978-1-4612-7449-0 978-1-4612-2184-5. DOI: 10.1007/978-1-4612-2184-5_1.

[28] E. Lin and A. Alessio, "What are the basic concepts of temporal, contrast, and spatial resolution in cardiac CT?" *J Cardiovasc Comput Tomogr*, vol. 3, no. 6, pp. 403–408, 2009, ISSN: 1934-5925. DOI: 10.1016/jJCCT.2009.07.003.

[29] K. R. Mathura, K. C. Vollebregt, K. Boer, J. C. De Graaff, D. T. Ubbink, and C. Ince, "Comparison of OPS imaging and conventional capillary microscopy to study the human microcirculation," en, *Journal of Applied Physiology*, vol. 91, no. 1, pp. 74–78, Jul. 2001, ISSN: 8750-7587, 1522-1601. DOI: 10.1152/jappl.2001.91.1.74.

[30] C. Somnin, L. Leclercq, J. Chamieh, et al., "Taylor dispersion analysis for measurement of diffusivity and size of gadolinium-based contrast agents," *European Journal of Pharmaceutical Sciences*, vol. 200, p. 106831, Sep. 2024, ISSN: 0928-0987. DOI: 10.1016/j.ejps.2024.106831.

[31] M. L. Huber, R. A. Perkins, A. Laesecke, et al., "New International Formulation for the Viscosity of H2O," en, *Journal of Physical and Chemical Reference Data*, vol. 38, no. 2, pp. 101–125, Jun. 2009, ISSN: 0047-2689, 1529-7845. DOI: 10.1063/1.3088050.

[32] W. Zhu, Y. Kato, and D. Artemov, "Water exchange-minimizing DCE-MRI protocol to detect changes in tumor vascular parameters: Effect of bevacizumab/paclitaxel combination therapy," en, *Magn Reson Mater Phy*, vol. 27, no. 2, pp. 161–170, Apr. 2014, ISSN: 1352-8661. DOI: 10.1007/s10334-013-0389-0.

[33] V. Goh, S. Halligan, and C. I. Bartram, "Quantitative Tumor Perfusion Assessment with Multidetector CT: Are Measurements from Two Commercial Software Packages Interchangeable?" *Radiology*, vol. 242, no. 3, pp. 777–782, Mar. 2007, Publisher: Radiological Society of North America, ISSN: 0033-8419. DOI: 10.1148/radiol.2423060279.

[34] M. S. Reed, M. Ochoa, K. M. Tichauer, A. M. Weichmann, M. M. Doyley, and B. W. Pogue, "Mapping estimates of vascular permeability with a clinical indocyanine green fluorescence imaging system in experimental pancreatic adenocarcinoma tumors," *JBO*, vol. 28, no. 7, p. 076001, Jul. 2023, Publisher: SPIE, ISSN: 1083-3668, 1560-2281. DOI: 10.1117/1.JBO.28.7.076001.

[35] D. R. de Waart, S. Häusler, M. L. H. Vlaming, et al., "Hepatic Transport Mechanisms of Cholyl-l-Lysyl-Fluorescein," *The Journal of Pharmacology and Experimental Therapeutics*, vol. 334, no. 1, pp. 78–86, Jul. 2010, ISSN: 0022-3565. DOI: 10.1124/jpet.110.166991.

[36] N. Hansen, "The CMA Evolution Strategy: A Comparing Review," en, in *Towards a New Evolutionary Computation: Advances in the Estimation of Distribution Algorithms*, ser. Studies in Fuzziness and Soft Computing, J. A. Lozano, P. Larrañaga, I. Inza, and E. Bengoetxea, Eds., Berlin, Heidelberg: Springer, 2006, pp. 75–102, ISBN: 978-3-540-32494-2. DOI: 10.1007/3-540-32494-1_4.

[37] M. Nomura and M. Shibata, *Cmaes: A Simple yet Practical Python Library for CMA-ES*, en, arXiv:2402.01373 [cs], Oct. 2024. DOI: 10.48550/arXiv.2402.01373.

[38] J. M. Hanna, P. Varsos, J. Kowalski, L. Sala, R. Meiburg, and I. E. Vignon-Clementel, *A comparative analysis of metamodels for lumped cardiovascular models, and pipeline for sensitivity analysis, parameter estimation, and uncertainty quantification*, en, arXiv:2410.12654 [math], Oct. 2024. DOI: 10.48550/arXiv.2410.12654.

[39] S. Turco, H. Wijkstra, and M. Mischi, "Mathematical Models of Contrast Transport Kinetics for Cancer Diagnostic Imaging: A Review," en, *IEEE Rev. Biomed. Eng.*, vol. 9, pp. 121–147, 2016, ISSN: 1937-3333, 1941-1189. DOI: 10.1109/REME.2016.2583541.

[40] E. Weidman, C. Foley, O. Kallas, et al., "Evaluating Permeability Surface-Area Product as a Measure of Blood-Brain Barrier Permeability in a Murine Model," *AJNR Am J Neuroradiol*, vol. 37, no. 7, pp. 1267–1274, Jul. 2016, ISSN: 0195-6108. DOI: 10.3174/ajnr.A4712.

[41] M. Aref, X. Ji, J. D. Handbury, K. Bailey, Z. Liang, and E. C. Wiener, "Spatial Resolution Effects on Permeability-Surface Area Estimation," en,

[42] A. Shukla-Dave, N. A. Obuchowski, T. L. Chenevert, et al., "Quantitative imaging biomarkers alliance (QIBA) recommendations for improved precision of DWI and DCE-MRI derived biomarkers in multicenter oncology trials," en, *Journal of Magnetic Resonance Imaging*, vol. 49, no. 7, e101–e121, 2019, _eprint: https://onlinelibrary.wiley.com/doi/pdf/10.1002/jmri.26518, ISSN: 1522-2586. DOI: 10.1002/jmri.26518.

[43] C. Falcó, "From random walks on networks to nonlinear diffusion," en, *Phys. Rev. E*, vol. 106, no. 5, p. 054103, Nov. 2022, arXiv:2206.08859 [physics], ISSN: 2470-0045, 2470-0053. DOI: 10.1103/PhysRevE.106.054103.

[44] D. B. Taulbee and C. P. Yu, "A theory of aerosol deposition in the human respiratory tract," en, *Journal of Applied Physiology*, vol. 38, no. 1, pp. 77–85, Jan. 1975, ISSN: 8750-7587, 1522-1601. DOI: 10.1152/jappl.1975.38.1.77.

[45] C. Darquenne and M. Paiva, "One-dimensional simulation of aerosol transport and deposition in the human lung," en, *J Appl Physiol* (1985), vol. 77, no. 6, pp. 2889–2898, Dec. 1994, ISSN: 8750-7587. DOI: 10.1152/jappl.1994.77.6.2889.

[46] J. Kowalski, L. Sala, D. Drasdo, and I. Vignon-Clementel, "Understanding Organ Vascular Behavior Through their Architecture," in *The 12th International Conference on Complex Networks and their Applications*, Menton, France, Nov. 2023.

# APPENDIX

## A. Mean transit time

We can rewrite (14) thanks to the properties of the Fourier transform

$$ \tau_{t,1} = i\tau \frac{\partial_t \hat{h}(1,0)}{\hat{h}(1,0) = 1} $$
$$
\begin{align*}
\tau_{t,1} &= i\tau \frac{\partial_\nu \hat{h}(1,0) \hat{c}_{\text{in}}(0)}{\hat{h}(1,0) \hat{c}_{\text{in}}(0)} \quad \text{as } \hat{c}_{\text{in}}(0) \neq 0 \\
&= i\tau \frac{(\hat{h}(1,\cdot)\hat{c}_{\text{in}})'(0) - \hat{h}(1,0)\hat{c}_{\text{in}}'(0)}{\hat{h}(1,0)\hat{c}_{\text{in}}(0)} \\
&= i\tau \frac{\hat{c}_{\text{out}}'(0)}{\hat{c}_{\text{out}}(0)} - i\tau \frac{\hat{c}_{\text{in}}'(0)}{\hat{c}_{\text{in}}(0)} \\
&= \frac{\int_t tc_{\text{out}}(t)dt}{\int_t c_{\text{out}}(t)dt} - \frac{\int_t tc_{\text{in}}(t)dt}{\int_t c_{\text{in}}(t)dt}.
\end{align*}
$$

B. Low Péclet number analysis

At low Péclet number, $k^2 = i\nu \frac{\tau_d}{\tau}$, which means that $k = (1 + \text{sign}(\nu)i)\sqrt{|\frac{\nu\tau_d}{2\tau}|}$ and

$$
\begin{align*}
\bar{h}(\nu) &= \frac{\tanh(k)}{k}, \\
\hat{h}(\chi, \nu) &= \frac{\cosh(k(1-\chi))}{\cosh k}.
\end{align*}
$$

In the low frequency range ($|\nu| \ll \frac{\tau}{\tau_d}$), the transfer function for the average concentration, and at the end of the vessel are respectively:

$$
\bar{h}(\nu) = \frac{1}{1 + i \frac{\nu \tau_d}{3\tau}}; \quad \hat{h}(1, \nu) = \frac{1}{1 + i \frac{\nu \tau_d}{2\tau}}. \qquad (30)
$$

By definition, the transfer function $\hat{h}(\nu)$ relates the unknown
concentration $\hat{c}(\nu)$ to the input concentration $\hat{c}_{\text{in}}(\nu)$: $\hat{c}(\nu) =
\hat{h}(\nu)\hat{c}_{\text{in}}(\nu)$. Approximations (30) are valid for frequencies
$|\nu| \ll \frac{\tau}{\tau_d}$. If the frequencies of the input signal are confined
below this threshold, then multiplying $\hat{h}$ by $\hat{c}_{\text{in}}$ leaves only
the effect of low frequencies in $\hat{h}$. This is the case for slowly
varying signals, such as in the case of slow drug uptake. In this
case, going back to the time-space domain through an inverse
Fourier transform gives:

$$
\bar{c}(\theta) + \frac{\tau_d}{3\tau} \frac{d\bar{c}}{d\theta}(\theta) = c_{\text{in}}(\theta) \quad (31)
$$

$$
c_{\text{out}}(\theta) + \frac{\tau_d}{2\tau} \frac{dc_{\text{out}}}{d\theta}(\theta) = c_{\text{in}}(\theta) \quad (32)
$$

Rearranging and redimensionalizing (31) gives

$$
\frac{d\bar{c}}{dt} = \frac{3}{\tau_d} (c_{\text{in}} - \bar{c})
$$

Besides, (31) - 2/3 (32) gives:

$$
\bar{c}(\theta) = \frac{2}{3}c_{\text{out}}(\theta) + \frac{1}{3}c_{\text{in}}(\theta) + \frac{\tau_d}{3\tau} \frac{d(\bar{c} - c_{\text{out}})}{d\theta}(\theta)
$$

Meaning that $\lambda$ equals $\frac{2}{3}$, on average over time.

C. Derivation of the AATH

The goal is to recover (21) and write the TK parameters
in terms of the physical quantities. We integrate (22) by
introducing $\tilde{t}(x, t) = t - \frac{\phi_{p,x}}{F_p L}$, which gives the integral form:

$$
c_p(x,t) = c_{\text{in}}(\tilde{t}(x,t)) + \frac{PS}{F_p} \int_0^x (c_e - c_p) \delta(x^* - L) dx^*.
$$

Introducing $\Delta c(t) = c_e(L,t) - c_p(L,t)$, and $\mathbb{H}$, the Heaviside step function, one gets:

$$
c_p(x, t) = c_{\text{in}}(\tilde{t}(x, t)) + \frac{PS}{F_p}\Delta c(t)\mathbb{H}(x - L). \quad (33)
$$

Equation (23) integrates as follows:

$$
\begin{align}
c_e(x < L, t) &= 0 \nonumber \\
\phi_e \frac{dc_e(L, \cdot)}{dt} &= PS(c_p(L, \cdot) - c_e(L, \cdot)) \tag{34}
\end{align}
$$

In the remainder of this paragraph, the study of the function $c_e(x, t)$ will therefore be restricted to the localization $x = L$. The plasma concentration (33) is evaluated at the vessel end ($x = L$):

$$
c_p(L, t) = c_{\text{in}}(\tilde{t}(L, t)) + \frac{PS}{F_p} (c_e(t) - c_p(L, t)) \mathbb{H}(0),
$$

and thus (regrouping the terms),

$$
c_p(L, t) (F_p + PS) = F_p c_{\text{in}}(\tilde{t}(L, t)) + PS c_e(t). \quad (35)
$$

We replace $c_p(L, t)$ in (34) by its value in (35):

$$
\phi_e \frac{dc_e}{dt} = \frac{PS \cdot F_p}{PS + F_p} (c_{\text{in}}(\tilde{t}(L, t)) - c_e(t)). \quad (36)
$$

Finally, evaluating (33) at $x = L^-$ gives:

$$
c_p(L^{-}, t) = c_{\text{in}}(\tilde{t}(L, t)), \qquad (37)
$$

which results in the AATH model:

$$
\begin{align}
\phi_p \frac{\partial c_p}{\partial t} &= -F_p L_p \frac{\partial c_p}{\partial x} \notag \\
\phi_e \frac{\partial c_e}{\partial t} &= \frac{PS \cdot F_p}{PS + F_p} (c_p(L) - c_e) . \tag{38}
\end{align}
$$

We define $K^{\text{trans}} = \frac{PS \cdot F_p}{PS + F_p}$, to recover (21) for all $x \le L$.

D. The residual error

For demonstration purposes, the residual of the 1C model
(9), for a generic function $f$, is defined as:

$$
\mathcal{R}_{1C}(f) = \frac{df}{d\theta} - \frac{\tau}{\lambda\tau_a}(c_{\text{in}} - f)
$$

$\mathcal{R}_{1C}(c_{1C}) = 0$ if $c_{1C}$ is a solution of (9), but this is not
necessarily the case for $c_f$ (in this case solution of (10)), hence
$\mathcal{R}_{1C}(c_f)$ may differ from zero:

$$
\mathcal{R}_{1C}(c_f) = \frac{d}{d\theta} \left( \int_0^1 c_f d\chi \right) - \frac{\tau}{\lambda\tau_a} \left( c_{\text{in}} - \int_0^1 c_f d\chi \right) \quad (39)
$$

However, as $c_f$ is a solution of (10), the residual of the full-order model $\mathcal{R}_f(c_f) = 0$. Therefore, on can subtract $\mathcal{R}_f(c_f)$ in (39) to get

$$
\mathcal{R}_f(c_f) = \frac{\partial c_f}{\partial \theta} + \frac{\tau}{\tau_a} \frac{\partial c_f}{\partial \chi} - \frac{\tau}{\tau_d} \frac{\partial^2 c_f}{\partial \chi^2} = 0.
$$

$$
\begin{align*}
\mathcal{R}_{1C}(c_f) &= \mathcal{R}_{1C}(c_f) - \int_0^1 \mathcal{R}_f(c_f) d\chi \\
&= \frac{\tau}{\lambda\tau_a} \left( (\lambda - 1)c_{\text{in}} - \lambda c_f(L) + \bar{c}_f \right) + \frac{\tau}{\tau_d} \partial_\chi c_f(0).
\end{align*}
$$