# Deterministic limit of a stochastic individual-based model of biological tissue interacting with diffusing chemicals

Léo Darrigade, Simon Labarthe, Béatrice Laroche

► To cite this version:

Léo Darrigade, Simon Labarthe, Béatrice Laroche. Deterministic limit of a stochastic individual-based model of biological tissue interacting with diffusing chemicals. 2025. hal-04336174v3

HAL Id: hal-04336174

https://hal.science/hal-04336174v3

Preprint submitted on 24 Feb 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
# Deterministic limit of a stochastic individual-based model of biological tissue interacting with diffusing chemicals

Léo Darrigade ⚪<sup>a,c</sup>, Simon Labarthe ⚪<sup>b,d</sup> and Béatrice Laroche ⚪<sup>a,c</sup>

<sup>a</sup> Université Paris-Saclay, INRAE, MaIAGE, 78350, Jouy-en-Josas, France; <sup>b</sup> Univ. Bordeaux, INRAE, BIOGECO, 33610, Cestas, France; <sup>c</sup> Université Paris-Saclay, Inria, MUSCA 91120, Palaiseau, France; <sup>d</sup> Inria, INRAE, PLEIADE, 33400, Talence, France;

## ARTICLE HISTORY

Compiled February 24, 2025

## ABSTRACT

Mathematical models of biological tissues are a promising tool for multi-scale data integration, computational experiments and system biology approaches. While some data and insights originate at the cellular level, macroscopic mechanisms emerge and are observed at the tissue scale, making tissue modeling an inherently multi-scale process. Consequently, tissue models can be broadly categorized into individual-based or continuous population-based.

In this paper, we first introduce a generic individual-based model of biological tissue including the main regulation processes such as cell division, differentiation, migration and death, along with cell-cell mechanical interactions. This cell population interacts with diffusing molecules, which are consumed or produced by cells and, in turn, regulate cellular behavior. This model is a measure-valued, piecewise-deterministic Markov process, coupled with reaction-diffusion PDEs. We establish the well-posedness of the model and rigorously derive its large-population deterministic limit.

## KEYWORDS

individual-based model (IBM); piecewise deterministic Markov process (PDMP); large population asymptotics; biological tissue model; interacting cells

## 1. Introduction

Biological tissues are meso-scale biological structures characterized by interacting cells with similar physiology embedded into an extracellular matrix. Tissue dynamics are governed by key regulatory mechanisms, including cell metabolism, cell-cell communications, cell movement and mechanical interactions [1]. Mathematical models of biological tissues are particularly useful for studying cell populations at the tissue scale, by facilitating multi-scale data integration or for testing hypothesis with high-throughput in silico experiments [2]. Consequently, a wide range of tissue models have been developed across diverse fields such as surgical robotics [3], embryogenesis [4, 5] or mathematical oncology [6].

Tissue models are typically classified into individual-based models (IBM) or continuous population models. IBMs treat cells as distinct entities by defining individual dynamics and interactions with other cells to study the emergence of patterns at the tissue scale [7]. These models align well with biological observations frequently made at the cellular level, making them suitable for integrating data and knowledge. Additionally, IBMs are well-adapted for
high-performance and parallel computing, and the incorporation of inter-individual variability in stochastic events is natural. However, their computational and memory demands scale at least linearly with the number of individual cells taken into account in the model, posing challenges for realistic simulations. In such cases, the synthetic data generated by IBMs can be as complex to analyze as actual experimental data [7].

In contrast, continuous population models consider the fate of the cell distribution at the population level. They focus on the dynamics of averaged quantities across the population, such as the mean of a quantitative trait or higher order momentum, using ordinary or partial differential equations (ODE or PDE). Scaling up the model at the tissue scale is then computationally tractable and ODE or PDE formalism makes them more amenable to mathematical analysis. However, these models lack key features of IBMs. They are most often fully deterministic (although stochastic events can technically be incorporated) whereas stochasticity is an important aspect of cellular behavior [8]. Additionally, they do not facilitate the modeling of individual trajectories hindering the comparison with individual-based data.

A rigorous link can be established between IBMs and continuous population models through large-population approximations [9]. This approach is particularly attractive for several reasons. It establishes a formal link between micro-scale IBM parameters and macro-scale population model parameters, facilitating the connection of measurements across scales. Additionally, the deterministic limit of an IBM can be used to accelerate computations in computationally intensive tasks such as parameter inference or sensitivity analysis. If qualitative mathematical analysis is impractical on IBMs due to technical constraints, it might be feasible on the deterministic approximation and then be extrapolated to the IBM for large enough populations. Such deterministic approximations have been studied in theoretical ecology and evolution [10–12], chemostat [13] or ovarian follicle models [14], but a comparable rigorous analysis remains lacking for tissue models.

The goal of this paper is twofold: we aim to propose a generic stochastic IBM of biological tissues within the PDMP (Piecewise Deterministic Markov Process) framework, that depicts a population of cells interacting with diffusing chemicals, alongside deriving a rigorous deterministic approximation using functional analysis tools. The proposed IBM is designed to be general enough for application across various tissue types, capturing essential aspects of biological tissue dynamics, including: 1) cell fate regulation, 2) spatial movement of cells, 3) stochastic cell division, extrusion and differentiation regulated by chemical concentrations, cell location and local cell density, 4) spatial diffusion of chemicals, 5) consumption or production of chemicals by cell metabolism. Our work extends previous studies by incorporating a broader range of biological phenomena and providing a strong theoretical foundation for models commonly used in mathematical oncology, morphogenesis and gut microbiota research [15].

The remainder of this paper is structured as follows. Section 2 introduces the custom notations used throughout. Section 3.1 presents the IBM model, with its stochastic part presented in section 3.1.1 and its deterministic part in section 3.1.2. Section 3.2 introduces a rescaled version of the model, followed by the derivation of its deterministic limit in the large population regime in section 3.3. This derivation builds on previous works [10, 12, 13] with adaptations to take into account particular features of biological tissues such as mechanical interactions between cells and the strong coupling of the cell population with reaction-diffusion PDEs.
## 2. Notations

For an integer $n \in \mathbb{N}^*$, let $[n] = \{1, \dots, n\}$. On the domain $\mathcal{X} = \mathcal{Z} \times [N_c]$, where $\mathcal{Z}$ represents space and $[N_c]$ the discrete cell types, we use the metric $d_{\mathcal{Z}}$:

$$d_{\mathcal{Z}}((z_1, l_1), (z_2, l_2)) = \begin{cases} |z_1 - z_2| & \text{if } l_1 = l_2, \\ \beta & \text{otherwise,} \end{cases} \quad (1)$$

with $\beta > \text{diam}(\mathcal{Z})$ and $|\cdot|$ the euclidean norm on $\mathcal{Z}$. Let $\mathcal{B}(\mathcal{X})$ be the Borel $\sigma$-algebra on $(\mathcal{X}, d_{\mathcal{Z}})$. We sometimes use the notation $\hat{\mathcal{Z}} := \{z-y : (z,y) \in \mathcal{Z}^2\}$.

Let $f$ be a function defined on $\mathcal{X}$. For all $l$ in $[N_c]$, we denote $f_l : z \mapsto f(z, l)$ the restriction of $f$ to $\mathcal{Z} \times \{l\}$. We write $f \in \mathcal{C}^{k, \bullet}(\mathcal{X}, \mathbb{R})$ iff, for all $l \in [N_c]$, $f_l$ is $k$ times differentiable on $\mathcal{Z}$, with $k$-th derivative continuous on $\mathcal{Z}$. When $k \ge 1$, we denote $\nabla_z f(x) \in \mathbb{R}^d$ the gradient operator with respect to the space variable $z$. We write $f \in \mathcal{C}^{k_1, k_2, \bullet}(\mathbb{R}^+ \times \mathcal{X}, \mathbb{R})$ iff, for all $l \in [N_c]$, $f_l(t, z)$ is continuous of class $k_1$ with respect to $t \in \mathbb{R}^+$ and continuous of class $k_2$ with respect to $z \in \mathcal{Z}$.

For $f: \mathcal{X} \to \mathbb{R}$, let $|f|_{\infty} := \sup_{x \in \mathcal{X}} |f(x)|$ and $|f|_{\text{Lip}} = \sup_{(x,y) \in \mathcal{X}, x \neq y} \frac{|f(x)-f(y)|}{d_{\mathcal{X}}(x,y)}$. Note that $f: \mathcal{X} \to \mathbb{R}$ is Lipschitz on $\mathcal{X}$ iff, $\forall l \in [N_c]$, $f_l$ is Lipschitz on $\mathcal{X}$. Moreover, $|f|_{\text{Lip}} \le \max(\max_{l \in [N_c]} |f_l|_{\text{Lip}(\mathcal{X})}, \frac{2|f|_{\infty}}{\beta})$ where $|f_l|_{\text{Lip}(\mathcal{X})}$ is the Lipschitz constant of $f_l$ on $\mathcal{X}$ and $\beta$ is defined in (1). For $f$ Lipschitz on $\mathcal{X}$, let $|f|_{\vee} = \max(|f|_{\infty}, |f|_{\text{Lip}})$ and let $\mathcal{L}ip_{1,1}(\mathcal{X}, \mathbb{R}) = \{f: \mathcal{X} \to \mathbb{R}, |f|_{\vee} \le 1\}$ be the space of functions bounded by 1 with Lipschitz constant lower than 1.

Let $\mathcal{M}_F(\mathcal{X})$ ($\mathcal{M}_{F^+}(\mathcal{X})$) be the space of finite (resp. finite positive) measures on $(\mathcal{X}, \mathcal{B}(\mathcal{X}))$, and $\mathcal{M}_P(\mathcal{X})$ be the subset of finite positive point measures $\mathcal{M}_P(\mathcal{X}) = \{\sum_{i=1}^n \delta_{x^i}; n \in \mathbb{N}, x^i = (z^i, l^i) \in \mathcal{Z} \times [N_c]\}$, with $\sum_{i=1}^0 \delta_{x^i}(x)$ the null measure by convention. For $v \in \mathcal{M}_F(\mathcal{X})$, $|\mu|_{\text{tv}}$ is the total variation norm; for $\mu, v \in \mathcal{M}_{F^+}(\mathcal{X})$, $d_{\text{PR}}(\mu, v)$ is the Prokhorov distance (see appendix D for definitions).

Let $v$ be a measure of $\mathcal{M}_F(\mathcal{X})$, then, for all $l \in [N_c]$, we define $v^l \in \mathcal{M}_F(\mathcal{X})$ as $v^l(Z) = v(Z \times \{l\})$, with $Z$ any Borel set of $\mathcal{Z}$. Then $v(dz, dl) = \sum_{k \in [N_c]} \delta_k(dx) v^k(dz)$ as for any set $X$ of $\mathcal{B}(\mathcal{X})$, $X = \cup_{l \in [N_c]} Z_l \times \{l\}$ with $Z_l$ a Borel set of $\mathcal{Z}$.

For a measurable function $f: \mathcal{X} \to \mathbb{R}$ and a measure $v \in \mathcal{M}_F(\mathcal{X})$, we define

$$\langle v, f \rangle := \int_{\mathcal{Z}} f(x) v(dx) = \sum_{l=1}^{N_c} \int_{\mathcal{Z}} f(z, l) v(dz, dl). \quad (2)$$

If $v \in \mathcal{M}_P(\mathcal{X})$, then $\langle v, f \rangle = \sum_{i=1}^n f(l^i, z^i)$.

For a measurable function $D: \mathbb{R}^d \times [N_c] \to \mathbb{R}$ and $v \in \mathcal{M}_F(\mathcal{X})$, we set for $x = (z, l) \in \mathcal{X}$, $v \ast D(x) := \sum_{l=1}^{N_c} \int_{\mathcal{Z}} D(z-y, l) v(dy, dl)$. If $v \in \mathcal{M}_P(\mathcal{X})$, $v \ast D(x) = \sum_{i=1}^{\langle v, 1 \rangle} D(z-z^i, l^i)$.
## 3. Main Results

### 3.1. A generic stochastic model of biological tissues

We now describe our generic individual based model of a biological tissue, expressed as a piecewise-deterministic (or hybrid) Markov process.

The spatial domain $\mathcal{X}$ is a compact subset of $\mathbb{R}^d$, with $d \le 3$. Cells are characterized by their center position in $\mathcal{X}$ and their type in $[N_c]$, the finite set of cell types. We model the population of cells $v$ as a point measure on $\mathcal{X} = \mathcal{X} \times [N_c]$:

$$v = \sum_{i=1}^{n} \delta_{x^i}(dx) \qquad (3)$$

where $x_i = (z_i, l_i) \in \mathcal{X}$ is the position and type the $i^{th}$ cell, $\delta_x(dx)$ is the Dirac measure on $x$ and $n = \langle v, 1 \rangle$ is the total number of cells in the population. Without explicitly specifying the shape of cells, we introduce $a$ as a typical length associated to cells, which may represent the cell radius or more generally the cell neighborhood.

The concentration of chemicals is modeled through a vector $c$ of $N_m \in \mathbb{N}^*$ non negative functions on the spatial domain $\mathcal{L}_a$ such that

(H1) $\mathcal{L}_a$ is a bounded open set of $\mathbb{R}^d$ of regularity $\mathcal{C}^\infty$, with boundary $\Gamma$, that satisfies

$$\mathcal{L}_a \supset \bigcup_{z \in \mathcal{X}} B_{\mathbb{R}^d}(z, a),$$

where $B_{\mathbb{R}^d}(z, r)$ is the open ball of $\mathbb{R}^d$ centered on $z$ of radius $r$.

$\mathcal{L}_a$ is an enlarged version of $\mathcal{X}$ which contains not only the position of the center of cells, as $\mathcal{X}$, but also the whole cells neighborhood and the bulk media where cells live.

#### 3.1.1. Stochastic model of cell fate

Cells can divide symmetrically, differentiate or die. These phenomena are modeled as independent random events, indexed by $[N_e]$ the finite set of event types. When at $t$ the cell $i$ undergoes a jump of type $k$, the population is modified according to

$$v_t = v_{t-} + \mu_k(x^i(t^-))$$

where $\mu_k: \mathcal{X} \to \mathcal{M}_{F^+}(\mathcal{X})$ and if $k$ is a

* death: $\mu_k(x) = -\delta_x$,

* differentiation of a cell of type $l$ into a type $l'$: $\mu_k(z, l) = -\delta_{(z,l)} + \delta_{(z,l')}$,

* division: a new daughter cell appears, located at $\lambda(z)$ from the mother cell: $\mu_k(z, l) = \delta_{(z+\lambda(z),l)}$.

For simplicity, we consider a deterministic division shift $\lambda$, however our results could be easily extended to a stochastic division shift, as done in [10, 12, 13].

For a cell at $x \in \mathcal{X}$, in a population $v \in \mathcal{M}_{F^+}(\mathcal{X})$, the individual jump rate associated to jump $k \in [N_e]$ takes the form

$$q_k^\infty \times q_k(x, v * D(z), c * \psi(z)) \qquad (4)$$
where $q_k : \mathcal{X} \times \mathbb{R} \times \mathbb{R}^{N_m} \to \mathbb{R}^+$. Here $D : \mathbb{R}^d \to \mathbb{R}^+$ is a kernel so that $v * D(z)$ accounts for the local cell density, $\psi$ is a kernel representing the cell neighborhood and $c * \psi(z) = (\int_{\mathcal{Z}_a} c_1(y) \psi(z-y) dy, \dots, \int_{\mathcal{Z}_a} c_{N_m}(y) \psi(z-y) dy)$ accounts for the concentrations of chemicals near the cell. $q_k^\infty < +\infty$ are constants. Therefore, we model a cell population where jump rates depend of the cell state $x$, the local density of the population and concentrations of chemicals near the cell. For instance, a high local density could promote death events as observed experimentally [16].

We assume the following regularity conditions.

(H2) $\lambda : \mathcal{L} \to \mathbb{R}^d$ is continuous and such that, $\forall z \in \mathcal{L}$, $z + \lambda(z) \in \mathcal{L}$.

(H3) $\lambda$ is locally Lipschitz.

(H4) For all $k \in [N_e]$, $q_k \le 1$.

(H5) For all $k \in [N_e]$, for all $l \in [N_c]$, $(z, y, c) \mapsto q_k(z, l, y, c)$ is Lipschitz on $\mathcal{L} \times \mathbb{R} \times \mathbb{R}$ with constant $K_{k,l}$. We define $K(q_k) := \max_{l \in [N_c]} K_{k,l}$.

(H6) $D \in \mathcal{C}^1(\mathbb{R}^d, \mathbb{R}^+)$.

(H7) $\psi$ is a positive function, Lipschitz on $\mathbb{R}^d$ with support contained in $B_{\mathbb{R}^d}(0, a)$ and its integral is 1.

### 3.1.2. Deterministic evolution of cells and chemical molecules

Between stochastic jumps, the state variables of the PDMP (cells and concentrations of chemicals) undergo a deterministic evolution that we detail now.

**Movement of cells.** For each cell $i$, the position of its center $z^i(t)$ evolves continuously in the spatial domain $\mathcal{L}$, in a non-inertial regime

$$ \frac{dz^i}{dt}(t) = \mathcal{F}(z^i(t), v_t * F(z^i(t)) - F(0)) \quad (5) $$

with

$$ F : \mathbb{R}^d \to \mathbb{R}^d \text{ and } \mathcal{F} : \mathcal{L} \times \mathbb{R}^d \to \mathbb{R}^d. $$

The first variable $z^i(t)$ in the resulting force $\mathcal{F}$ is introduced to take into account local modulations (if any). The second variable is the total force exerted by other cells. The cell $j$ at position $z^j$ exerts on cell $i$ at position $z^i$ a force $F(z^i - z^j)$. This force only depends on spatial distance and direction between cells, and not on cell types. The total force exerted on cell $i$ is

$$ \sum_{j=1, j \neq i}^{n} F(z^i - z^j) = \sum_{j=1}^{n} F(z^i - z^j) - F(0) = v * F(z^i) - F(0). \quad (6) $$

We impose the following regularity assumptions.

(H8) $F \in \mathcal{C}^1(\mathbb{R}^d, \mathbb{R}^d)$, $\mathcal{F} \in \mathcal{C}^1(\mathcal{L} \times \mathbb{R}^d, \mathbb{R}^d)$ and there exists continuous $K_{\mathcal{F}} : \mathcal{L} \to \mathbb{R}_+$ such that $|\mathcal{F}(z, y_1) - \mathcal{F}(z, y_2)| \le K_{\mathcal{F}}(z)|y_1 - y_2|$.

(H9) $\mathcal{F}(z, \cdot) = 0 \quad \forall z \in \partial \mathcal{L}$.

Thanks to hypothesis (H8) and (H9) the system of coupled equations (5) has a unique solution on $\mathbb{R}^+$, taking values in $\mathcal{L}^n$, for any initial condition in $\mathcal{L}^n$. Note that (H9) could be weakened, by imposing tangent forces on the domain boundary for instance, here we chose null forces on the boundary for simplicity.
**Reaction-diffusion of molecules.** Let $c(z)$ be the vector of concentrations of the $N_m$ molecules located at $z$. Each cell type $l \in [N_c]$ can degrade or produce chemical $m \in [N_m]$, at speed $\gamma_{ml}(c(z))$. A cell $i$ located at $z^i$ can catalyze a reaction at $z$ only if it is physically close to $z$. Therefore, we weight the reaction speed $\gamma_{ml}(c(z))$ with a term $\psi(z-z^i)$ where $\psi$ was introduced in hypothesis (H7). On the boundary of $\mathcal{X}_a$, concentrations are influenced by the external media through quantities $f_m$ which depend of time and space.

Then, concentrations $c$ on $\mathcal{X}_a$ evolve according to the following system of coupled reaction-diffusion equations

$$ 
\begin{array}{l}
\forall 1 \le m \le N_m, \quad \partial_t c_m - \sigma \Delta c_m = \sum_{l=1}^{N_c} \gamma_{ml}(c) v_l^l * \psi(z), \qquad (7) \\
\qquad \partial_n c_{ml}[0,T] \times \Gamma + b_m c_{ml}[0,T] \times \Gamma = f_m, \\
\qquad c_m(\cdot, 0) = c_{0m},
\end{array}
$$

where $\partial_n$ is the normal derivative on $\Gamma$.

For an open set $U \subset \mathbb{R}^d$ and $\alpha \ge 0$, let

$$ W_{L^2(U)}^{2\alpha, [\alpha]} := \left\{ f \in L^2([0,T], H^{2\alpha}(U)) ; \forall j \in \{1, \dots, [\alpha]\} \quad \frac{\partial^j f}{\partial t^j} \in L^2([0,T], H^{2\alpha-2j}(U)) \right\} $$

with $[\cdot]$ the ceiling function.

We assume that for all $m \in [N_m]$,

(H10) $b_m$ is in $C^\infty(\mathcal{X}_a)$,

(H11) $\forall l \in [N_c]$, $\gamma_{ml}$ is bounded and uniformly Lipschitz from $\mathbb{R}_+^{N_m}$ to $\mathbb{R}$,

(H12) the initial condition $c_{0m}$ is in $H^1(\mathcal{X}_a)$,

(H13) $f_m$ is in $W_{L^2(\Gamma)}^{\frac{1}{2},1}$.

Moreover, to ensure nonnegative chemical concentrations, we assume that

(H14) $c_{0m} \ge 0$ a.e. and $\gamma_{ml}(c) \ge 0$ for all $c$ such that $c_m = 0$.

**Proposition 3.1.** (1) Let $T \ge 0$, for all $v$ measurable on $[0,T]$ with values in $M_{F^+}(\mathcal{X})$ endowed with the weak convergence topology and the associated $\sigma$-algebra, and such that $\sup_{t \in [0,T]} |v_t|_{tv} < +\infty$, the non-linear evolution problem (7) has a unique non-negative solution $c$ in $(W_{L^2(\mathcal{X}_a)}^{2,1})^{N_m} \cap \mathcal{C}([0,T], H^1(\mathcal{X}_a)^{N_m})$.

(2) Let $T \ge 0$, $v^1$, $v^2$ be two measurable functions from $[0,T]$ to $M_{F^+}(\mathcal{X})$ such that $\sup_{t \in [0,T]} |v_t^1|_{tv}$ and $\sup_{t \in [0,T]} |v_t^2|_{tv}$ are bounded by a positive constant $C$, let $c^1$ et $c^2$ be the associated solutions of (7), then

$$ \sup_{t \le T} \|c^1(t) - c^2(t)\|_{H^1(\mathcal{X}_a)^{N_m}} \le K_1 e^{K_2 CT} \sum_{l=1}^{N_c} \int_0^T \|(v_l^{1l} - v_l^{2l}) * \psi\|_{L^2(\mathcal{X}_a)} dt \quad (8) $$

with $K_1, K_2$ two positive constant independent of $v^1, v^2$.

**Proof.** The proof of proposition 3.1 is done in Appendix A. $\square$
### 3.2. Rescaled process

We introduce a sequence of processes $(v_t^N, c_t^N)_{N \in \mathbb{N}^*}$ similar to the model described in section 3.1 with parameters $F^N, D^N, \gamma^N$ such that

$$F^N = \frac{F}{N}, \quad D^N = \frac{D}{N}, \quad \gamma^N = \frac{\gamma}{N}, \qquad (9)$$

leaving the other model parameters unchanged. Note that $N = 1$ corresponds to the process described in section 3.1. The sequence of rescaled processes is defined as $(\bar{v}_t^N, c_t^N)_{N \in \mathbb{N}^*}$, with

$$\bar{v}_t^N := \frac{1}{N} v_t^N.$$

The global existence and uniqueness properties of cell spatial trajectories and concentrations of chemicals as stated in section 3.1.2 readily extend to $(v_t^N, c_t^N)$ and $(\bar{v}_t^N, c_t^N)$, for all $N \in \mathbb{N}^*$, so that without stochastic jumps (or between jumps), the deterministic flow of these processes is well defined. Moreover, from (9), we have $v_t^N * F^N(z) - F^N(0) = \bar{v}_t^N * F(z) - F(0)/N$, $v_t^N * D^N(z) = \bar{v}_t^N * D(z)$ for all $z \in \mathcal{Z}$, and $\gamma_{ml}^N(c_t^N) v_t^{IN} * \psi(z) = \gamma_{ml}(c_t^N) \bar{v}_t^{IN} * \psi(z)$, for all $(m, l) \in [N_m] \times [N_c]$.

We now introduce independent Poisson random measure on $[0, 1] \times \mathbb{N}^* \times \mathbb{R}^+$, one for each type of jump $k \in [N_e]$, denoted $N_k(d\theta, di, dt)$, with intensity $n_k(d\theta, di, dt) = q_k^\infty d\theta dt \sum_{k \in \mathbb{N}^*} \delta_k(di)$. Let us consider $T > 0$, provided trajectories exist on $[0, T]$, the process $(\bar{v}_t^N, c_t^N) \in \mathcal{D}([0, T], \mathcal{M}_{F^+}(\mathcal{X})) \times \mathcal{C}([0, T], H^1(\mathcal{L}_a)^{N_m})$ is such that for all $f \in \mathcal{C}^{1,1,\bullet}(\mathbb{R}^+ \times \mathcal{X}, \mathbb{R})$,

$$\langle \bar{v}_t^N, f_t \rangle = \langle \bar{v}_0^N, f_0 \rangle \\ + \int_0^t \langle \bar{v}_s^N, \mathcal{F}(\cdot, \bar{v}_s^N * F(\cdot) - F(0)/N) \nabla_z f_s(\cdot) + \partial_s f_s(\cdot) \rangle ds \\ + \sum_{k=1}^{N_e} \int_0^t \int_{\mathbb{N}^*} \int_0^1 \left[ \langle \bar{v}_{s^{-}}^N + \mu_k(x_{s^{-}}^{Ni})/N, f_s \rangle - \langle \bar{v}_{s^{-}}^N, f_s \rangle \right] \\ \times \mathbb{1}\{i \le \langle v_{s^{-}}, 1 \rangle\} (i) \mathbb{1}\{\theta \le q_k(x_{s^{-}}^{Ni}, \bar{v}_{s^{-}}^N * D(z_{s^{-}}^{Ni}), c_s^N * \psi(z_{s^{-}}^{Ni}))\} (\theta) N_k(d\theta, di, ds), \qquad (10)$$

where $\mathcal{F}(\cdot) \nabla_z f_s(\cdot)$ is a scalar product between vectors of $\mathbb{R}^d$, and, for chemical $m \in [N_m]$,

$$\partial_t c_m^N - \sigma \Delta c_m^N = \sum_{l=1}^{N_c} \gamma_{ml}(c_t^N) \bar{v}_t^{IN} * \psi(z), \qquad (11)$$

$$\partial_n c_{m|\lfloor 0,T \rfloor \times \Gamma}^N + b_m c_{m|\lfloor 0,T \rfloor \times \Gamma}^N = f_m,$$

$$c_m^N(\cdot, 0) = c_{0m}.$$

For $(\bar{v}_t^N, c_t^N)$ to have well defined trajectories on $\mathbb{R}^+$ when stochastic events occur, the process needs to be stable, meaning that an infinite number of jumps cannot accumulate in finite time. This is proved in the next proposition.

**Proposition 3.2.** *The following assertions are true:*
(1) Let $t \ge 0$ and $p \in \mathbb{N}^*$. If $\sup_{N \in \mathbb{N}^*} \mathbb{E}[\langle \bar{\mathbf{v}}_0^N, 1 \rangle^p] < +\infty$, then

$$ \sup_{N \in \mathbb{N}^*} \mathbb{E} \left[ \sup_{s \le t} \langle \bar{\mathbf{v}}_s^N, 1 \rangle^p \right] \le C(p,t) < +\infty \quad (12) $$

with $C(p,t)$ a constant that depends only on $p$ and $t$.

(2) For all $N \in \mathbb{N}^*$, the process $(c_t^N, \bar{\mathbf{v}}_t^N)$ solution of (10-11) is stable and trajectories exist on $\mathbb{R}^+$.

**Proof.** For the first assertion, let $p, N \in \mathbb{N}^*$. For all $n \in \mathbb{N}$, we define the stopping-time $\tau_n^N = \inf\{t \ge 0, \langle \mathbf{v}_t^N, 1 \rangle = n\}$. The process $(\bar{\mathbf{v}}_{t \wedge \tau_n^N}^N, c_{t \wedge \tau_n^N}^N)_t$ is well-defined as the cell population is bounded by $n$ and, by proposition 3.1, $c^N$ exists.

Note that the deterministic evolution does not change population size, and jumps modify it at most by one. Let $t \ge 0$, applying Itô formula (theorem 5.1 chapter II in Ref. [17]) to $\langle \bar{\mathbf{v}}_{t \wedge \tau_n^N}^N, 1 \rangle^p$, we get

$$ \langle \bar{\mathbf{v}}_{t \wedge \tau_n^N}^N, 1 \rangle^p \le \langle \bar{\mathbf{v}}_0^N, 1 \rangle^p + \sum_{k=1}^{N_e} \int_0^{t \wedge \tau_n^N} \int_{\mathbb{N}^*} \int_0^1 \left[ \left( \langle \bar{\mathbf{v}}_{s-}^N, 1 \rangle + \frac{1}{N} \right)^p - \langle \bar{\mathbf{v}}_{s-}^N, 1 \rangle^p \right] \mathbb{1}_{\{i \le \langle \mathbf{v}_{s-}^N, 1 \rangle\}} (i) \times \mathbb{1}_{\{\theta \le q_k(\cdot)\}} (\theta) N_k(d\theta, di, ds). $$

For positive $x$, $(x + \frac{1}{N})^p - x^p \le \frac{C(p)}{N^p}((Nx)^{p-1} + 1) \le \frac{C(p)}{N}(x^{p-1} + 1)$ with $C(p) = \sum_{k=1}^{p-1} \binom{p}{k}$. Therefore,

$$ \sup_{s \le t \wedge \tau_n} \langle \bar{\mathbf{v}}_s^N, 1 \rangle^p \le \langle \bar{\mathbf{v}}_0^N, 1 \rangle^p + \sum_{k=1}^{N_e} \int_0^{t \wedge \tau_n^N} \int_{\mathbb{N}^*} \int_0^1 \frac{C(p)}{N} \left( \langle \bar{\mathbf{v}}_{s-}^N, 1 \rangle^{p-1} + 1 \right) \mathbb{1}_{\{i \le \langle \mathbf{v}_{s-}^N, 1 \rangle\}} (i) \times \mathbb{1}_{\{\theta \le q_k(\cdot)\}} (\theta) N_k(d\theta, di, ds) \\ \Rightarrow \mathbb{E} \left[ \sup_{s \le t \wedge \tau_n} \langle \bar{\mathbf{v}}_s^N, 1 \rangle^p \right] \le \mathbb{E} \left[ \langle \bar{\mathbf{v}}_0^N, 1 \rangle^p \right] + \sum_{k=1}^{N_e} \mathbb{E} \left[ \int_0^t q_k^\infty C(p) \left( \langle \bar{\mathbf{v}}_{s \wedge \tau_n^N}^N, 1 \rangle^p + \langle \bar{\mathbf{v}}_{s \wedge \tau_n^N}^N, 1 \rangle \right) ds \right] \\ \le \mathbb{E} \left[ \langle \bar{\mathbf{v}}_0^N, 1 \rangle^p \right] + \sum_{k=1}^{N_e} \mathbb{E} \left[ \int_0^t q_k^\infty C(p) \left( \langle \bar{\mathbf{v}}_{s \wedge \tau_n^N}^N, 1 \rangle^p + \langle \bar{\mathbf{v}}_{s \wedge \tau_n^N}^N, 1 \rangle \right) ds \right] \\ \le \mathbb{E} \left[ \langle \bar{\mathbf{v}}_0^N, 1 \rangle^p \right] + 2C(p) \left( \sum_{k=1}^{N_e} q_k^\infty \right) \int_0^t \mathbb{E} \left[ \sup_{u \le s \wedge \tau_n^N} \langle \bar{\mathbf{v}}_u^N, 1 \rangle^p \right] ds. $$

as for $p \ge 1$, $n + n^p \le 2n^p$. Using Gronwall lemma and our hypotheses, we can find a bound which does not depend of $n$ nor $N$:

$$ \mathbb{E} \left[ \sup_{s \le t \wedge \tau_n^N} \langle \bar{\mathbf{v}}_s^N, 1 \rangle^p \right] \le \mathbb{E} \left[ \langle \bar{\mathbf{v}}_0^N, 1 \rangle^p \right] \exp \left( 2C(p) \left( \sum_{k=1}^{N_e} q_k^\infty \right) t \right) = C(p,t). \quad (13) $$
Therefore $\tau_n \to +\infty$ a.s. Then, by Fatou lemma

$$ \mathbb{E} \left[ \sup_{s \le t} \langle \bar{\mathbf{v}}_s^N, 1 \rangle^p \right] = \mathbb{E} \left[ \liminf_{n \to \infty} \sup_{s \le t \wedge \tau_n^N} \langle \bar{\mathbf{v}}_s^N, 1 \rangle^p \right] \le \liminf_{n \to \infty} \mathbb{E} \left[ \sup_{s \le t \wedge \tau_n^N} \langle \bar{\mathbf{v}}_s^N, 1 \rangle^p \right] \le C(p,t), $$

which proves the desired result.

A direct consequence of assertion 1 is that the population does not explode in finite time a.s. Individual jump rates being bounded, as long as the population is finite, the total jump rate is finite and only a finite number of jumps occurs on finite time interval, which proves the second assertion. $\square$

We now establish a martingale property which will be used in our convergence study. To do so, we introduce for each event type $k \in [N_e]$ the compensated random Poisson measure $\tilde{N}_k(d\theta, di, ds) = N_k(d\theta, di, ds) - n_k(d\theta, di, ds)$.

**Proposition 3.3.** Assume that $\mathbb{E}[\langle \bar{\mathbf{v}}_0^N, 1 \rangle^2] < +\infty$. For all $f \in \mathcal{C}^{1,1,\bullet}(\mathbb{R}^+ \times \mathcal{X}, \mathbb{R})$ and $t \ge 0$,

$$ \begin{aligned} \langle \bar{\mathbf{v}}_t^N, f_t \rangle = \langle \bar{\mathbf{v}}_0^N, f_0 \rangle &+ \int_0^t \langle \bar{\mathbf{v}}_s^N, \mathcal{F}(\cdot, \bar{\mathbf{v}}_s^N * F(\cdot) - F(0)/N) \nabla_z f_s + \partial_s f_s \rangle ds \quad (14) \\ &+ \sum_{k=1}^{N_e} q_k^\infty \int_0^t \int_{\mathcal{X}} \langle \mu_k(x), f_s \rangle q_k(x, \bar{\mathbf{v}}_s^N * D(z), c_s^N * \psi(z)) \bar{\mathbf{v}}_s^N(dx) ds \\ &+ Z_t^{f,N} \end{aligned} $$

where $Z_t^{f,N}$ is the square integrable martingale

$$ Z_t^{f,N} = \sum_{k=1}^{N_e} \int_0^t \int_{\mathbb{N}^*} \int_0^1 \langle \mu_k(x_s^{Ni}), f_s/N \rangle \mathbb{1}_{\{i \le \langle v_s^N, 1 \rangle\}} (i) \mathbb{1}_{\{\theta \le q_k(\cdot)\}} (\theta) \tilde{N}_k(d\theta, di, ds), \quad (15) $$

with quadratic variation

$$ \langle Z_t^{f,N} \rangle_t = \frac{1}{N} \sum_{k=1}^{N_e} q_k^\infty \int_0^t \int_{\mathcal{X}} \langle \mu_k(x), f_s \rangle^2 q_k(x, \bar{\mathbf{v}}_s^N * D(z), c_s^N * \psi(z)) \bar{\mathbf{v}}_s^N(dx) ds. \quad (16) $$

**Proof.** The integrands in (15) are left-continuous and therefore predictable. From proposition 3.2, for all $t \ge 0$,

$$ \begin{aligned} \mathbb{E} \left[ \sum_{k=1}^{N_e} \int_0^t \int_{\mathbb{N}^*} \int_0^1 \left| \langle \mu_k(x_s^{Ni}), f_s/N \rangle \right| \mathbb{1}_{\{i \le \langle v_s^N, 1 \rangle\}} (i) \mathbb{1}_{\{\theta \le q_k(\cdot)\}} (\theta) n_k(d\theta, di, ds) \right] \\ = \sum_{k=1}^{N_e} \int_0^t \int_{\mathcal{X}} \mathbb{E} \left[ q_k^\infty \left| \langle \mu_k(x), f_s/N \rangle \right| q_k(x, \bar{\mathbf{v}}_s^N * D(z), c_s^N * \psi(z)) \bar{\mathbf{v}}_s^N(dx) ds \right] < +\infty \end{aligned} $$

and

$$ \mathbb{E} \left[ \sum_{k=1}^{N_e} \int_0^t \int_{\mathbb{N}^*} \int_0^1 \left| \langle \mu_k(x_s^{Ni}), f_s/N \rangle \right|^2 \mathbb{1}_{\{i \le \langle v_s^N, 1 \rangle\}} (i) \mathbb{1}_{\{\theta \le q_k(\cdot)\}} (\theta) n_k(d\theta, di, ds) \right] < +\infty. $$
We can then apply the classical result given p.62 of Ref. [17].

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

### 3.3. Deterministic limit of the rescaled model with a large population approximation.

Let $T$ be a positive time. Let $\mathcal{D}([0,T], \mathcal{M}_{F^+}(\mathcal{X}))$ be the space of càdlàg functions of $[0,T]$ to $\mathcal{M}_{F^+}(\mathcal{X})$. We now prove the weak convergence of the rescaled process sequence $(\bar{v}_t^N, c_t^N),_{N \in \mathbb{N}^*}$ defined in section 3.2 towards a deterministic process in the space

$$ \mathcal{D}([0,T], \mathcal{M}_{F^+}(\mathcal{X})) \times \mathcal{C}([0,T], H^1(\mathcal{Z}_a)^{N_m}). $$

We endow $\mathcal{D}([0,T], \mathcal{M}_{F^+}(\mathcal{X}))$ with the Skorokhod metric, using the Prokhorov distance $d_{PR}$ on $\mathcal{M}_{F^+}(\mathcal{X})$. The topology associated to the Prokhorov metric is the topology of the weak convergence on $\mathcal{M}_{F^+}(\mathcal{X})$. For more details, see appendix D.

We suppose that regularity hypotheses (H1)–(H14) are verified. Additionally, we assume that:

(H15) the sequence $\bar{v}_0^N = \frac{1}{N}v_0^N$, seen as a sequence of random variables on $\mathcal{M}_{F^+}(\mathcal{X})$, converges weakly to $\xi_0 \in \mathcal{M}_{F^+}(\mathcal{X})$.

(H16) There exists an integer $p \ge 2$ such that $\sup_{N \in \mathbb{N}^*} \mathbb{E}[\langle \bar{v}_0^N, 1 \rangle^p] < +\infty$.

We can now state the main result of the paper.

**Theorem 3.4.** Under all the hypotheses mentioned above, the process $(\bar{v}_t^N, c_t^N)_{t \ge 0}$ converges weakly in $\mathcal{D}([0,T], \mathcal{M}_{F^+}(\mathcal{X})) \times \mathcal{C}([0,T], H^1(\mathcal{Z}_a)^{N_m})$ to a limit $(\xi_t, c_t)_{t \ge 0}$, where $(\xi_t, c_t)_{t \ge 0}$ is a. s. the unique solution in $\mathcal{C}([0,T], \mathcal{M}_{F^+}(\mathcal{X})) \times \mathcal{C}([0,T], H^1(\mathcal{Z}_a)^{N_m})$ of

$$ \langle \xi_t, f_t \rangle = \langle \xi_0, f_0 \rangle + \int_0^t \langle \xi_s, \mathcal{F}(\cdot, \xi_s * F(\cdot)) \nabla_z f_s + \partial_t f_s \rangle ds \\ + \sum_{k=1}^{N_c} \int_0^t \sum_{l=1}^{N_c} \int_{\mathcal{Z}} \langle \mu_k(z,l), f_s \rangle q_k^\infty q_k((z,l), \xi_s * D(z), c_s * \psi(z)) \xi_s(dz, dl) ds \quad (17) $$

for all $f \in \mathcal{C}^{1,1,\bullet}(\mathbb{R}^+ \times \mathcal{X}, \mathbb{R})$, and $\forall 1 \le m \le N_m$,

$$ \partial_t c_m - \sigma \Delta c_m = \sum_{l=1}^{N_c} \gamma_{ml}(c) \xi_l^l * \psi(z), \quad (18) $$

$$ \partial_n c_{m|_{[0,T] \times \Gamma}} + b_m c_{m|_{[0,T] \times \Gamma}} = f_m, $$

$$ c_m(\cdot, 0) = c_{0m}. $$

This result deserves some comments. We can first note that the result is still valid if the initial condition $\xi_0$ is random and $(\bar{v}_0^N)_N$ converges weakly to $\xi_0$. We can also observe that, setting $Q_{k,l}(z,l) = q_k^\infty q_k((z,l), \xi_l * D(z), c_l * \psi(z))$, the strong formulation of (17) takes the form of a system of nonlocal PDEs, where for all $l \in [N_c]$,

$$ \begin{aligned} \partial_t \xi_l^l(z) - \mathrm{Div} \left( \mathcal{F}(z, \xi_l * F(z)) \xi_l^l(z) \right) = & - \left( Q_{(\mathrm{death} \ l),l}(\cdot, \cdot) \times \xi_l^l \right) (z) \\ & + \left( Q_{(\mathrm{div.} \ l),l}(\cdot, \cdot) \times \xi_l^l \right) (z - \lambda(z)) \\ & - \sum_{l'=1}^{N_c} \left( Q_{(l \ \mathrm{to} \ l'),l}(\cdot, \cdot) \times \xi_l^{l'} \right) (z) \\ & + \sum_{l'=1}^{N_c} \left( Q_{(l' \ \mathrm{to} \ l),l}(\cdot, \cdot) \times \xi_l^{l'} \right) (z) \end{aligned} $$
where (death $l$), (div. $l$) and ($l$ to $l'$) denotes respectively the events (death of a cell of type $l$), (division of a cell of type $l$) and (differentiation of a cell of type $l$ into type $l'$). This system is reminiscent of aggregation equations through its left handside [18], frequently used to model collective movements such as particle aggregation, repulsion or swarming. An original feature comes from the presence of reaction terms on the right handside and the strong coupling with reaction-diffusion equations, which is not commonly encountered in the literature to our knowledge.

The proof of theorem 3.4 is established through three successive steps. First, let us prove that the sequence of laws of $(\bar{v}^N, c^N)_N$ has an accumulation point (lemma 3.5 and 3.6).

**Lemma 3.5.** The sequence of laws of $(\bar{v}^N)_{N \in \mathbb{N}^*}$, considered as probability measures on $\mathcal{D}([0, T], \mathcal{M}_{F^+}(\mathcal{X}))$, is tight.

*Proof.* Based on theorem 2.1 in [19], the tightness of $(\bar{v}^N)_{N \in \mathbb{N}^*}$ is equivalent to the tightness of laws of $(\langle \bar{v}^N, f \rangle)_{N \in \mathbb{N}^*}$ on $\mathcal{D}([0, T], \mathbb{R})$, for all $f$ in a dense set of $\mathcal{C}(\mathcal{X}, \mathbb{R})$ endowed with the topology of uniform convergence. Here, we choose $f \in \mathcal{C}^{1, \bullet}(\mathcal{X}, \mathbb{R})$. A sufficient condition for this property is given by the Aldous-Rebolledo criteria (corollary 2.3.3 in [20]) according to which we have to establish that

(1) for $t \ge 0$, the sequence $(\langle \bar{v}_t^N, f \rangle)_{N \in \mathbb{N}^*}$ is tight on $\mathbb{R}$,

(2) considering the semi-martingale decomposition defined at prop. 3.3 $\langle \bar{v}_t^N, f \rangle = \langle \bar{v}_0^N, f \rangle + A_t^{f,N} + Z_t^{f,N}$ with $A_t^{f,N}$ a finite variation process and $Z_t^{f,N}$ a martingale, then for all $t \ge 0$, $\varepsilon > 0$, $\eta > 0$, there exists $\delta > 0$ and $N_0 \in \mathbb{N}^*$ such that for all sequence of stopping times $(\tau_N)_{N \in \mathbb{N}^*}$ ($\tau_N$ is a $\mathcal{F}^N$-stopping time) with $\tau_N \le t$ we have

$$ \sup_{N \ge N_0} \sup_{0 \le \theta \le \delta} \mathbb{P}(|A_{\tau_N+\theta}^{f,N} - A_{\tau_N}^{f,N}| \ge \eta) \le \varepsilon, \quad (19) $$

$$ \sup_{N \ge N_0} \sup_{0 \le \theta \le \delta} \mathbb{P}(|\langle Z^{f,N} \rangle_{\tau_N+\theta} - \langle Z^{f,N} \rangle_{\tau_N}| \ge \eta) \le \varepsilon. \quad (20) $$

To prove assertion 1, we consider some $f \in \mathcal{C}^{1, \bullet}(\mathcal{X}, \mathbb{R})$ and $t \ge 0$. Under hypotheses (H15)–(H16) and from proposition 3.2.1, $\sup_{N \in \mathbb{N}^*} \mathbb{E}[\langle \bar{v}_t^N, 1 \rangle] \le C(p,t) + 1$. From the Markov inequality, for any $N$,

$$ \mathbb{P}(\langle \bar{v}_t^N, f \rangle \ge K) \le \frac{1}{K} |f|_{\infty} \mathbb{E}[\langle \bar{v}_t^N, 1 \rangle] \le \frac{1}{K} |f|_{\infty} (C(p,t) + 1). $$

Therefore, the sequence $(\langle \bar{v}_t^N, f \rangle)_{N \in \mathbb{N}^*}$ is tight.

To prove assertion 2, let us consider $t, \eta, \varepsilon, \delta \ge 0$, and $\theta \le \delta$, $\tau_N < t$. From proposition 3.3,

$$ \mathbb{E} \left[ |A_{\tau_N+\theta}^{f,N} - A_{\tau_N}^{f,N}| \right] \le \mathbb{E} \left[ \left| \int_{\tau_N}^{\tau_N+\theta} \langle \bar{v}_s^N, \mathcal{F}(\cdot, \bar{v}_s^N * F(\cdot) - F(0)/N) \nabla_z f \rangle ds \right| \right. \\ \left. + \sum_{k=1}^{N_e} \mathbb{E} \left[ \left| \int_{\tau_N}^{\tau_N+\theta} \int_{\mathcal{X}} q_k(x, \bar{v}_s^N * D(z), c_s^N * \psi(z)) \langle \mu_k(x), f \rangle \bar{v}_s^N(dx) ds \right| \right] \right] $$
From hypothesis (H8),

$$ \begin{aligned} \left| \langle \bar{\nu}_s^N, \mathcal{F}(\cdot, \bar{\nu}_s^N * F(\cdot) - F(0)/N) \nabla_z f \rangle \right| &\leq \left| \nabla_z f \right|_{\infty} \left| F_{|\mathcal{F}|} \right|_{\infty} \left| K_{\mathcal{F}} \right|_{\infty} \langle \bar{\nu}_s^N, 1 \rangle^2 \\ &\quad + \left| \nabla_z f \right|_{\infty} \left| \mathcal{F}(\cdot, -F(0)/N) \right|_{\mathcal{Z}} \left| \bar{\nu}_s^N, 1 \right\rangle. \end{aligned} $$

Setting $C_1 := |\nabla_z f|_{\infty} |F_{|\mathcal{F}|}|_{\infty} |K_{\mathcal{F}}|_{\infty}$ and $C_2 := |\nabla_z f|_{\infty} |\mathcal{F}(\cdot, -F(0)/N)|_{\mathcal{Z}}|_{\infty} + \sum_{k=1}^{N_c} 2q_k^{\infty} |f|_{\infty}$, we obtain

$$ \begin{aligned} \mathbb{E} \left[ |A_{\tau_N+\theta}^{f,N} - A_{\tau_N}^{f,N}| \right] &\leq C_1 \mathbb{E} \left[ \int_{\tau_N}^{\tau_N+\theta} \langle \bar{\nu}_s^N, 1 \rangle^2 ds \right] \\ &\quad + C_2 \mathbb{E} \left[ \int_{\tau_N}^{\tau_N+\theta} \langle \bar{\nu}_s^N, 1 \rangle ds \right]. \end{aligned} $$

We have that $\int_{\tau_N}^{\tau_N+\theta} \langle \bar{\nu}_s^N, 1 \rangle^2 ds \leq \delta \sup_{s \leq t} \langle \bar{\nu}_s^N, 1 \rangle^2$. So, applying proposition 3.2.1, we get that

$$ \mathbb{E} \left[ |A_{\tau_N+\theta}^{f,N} - A_{\tau_N}^{f,N}| \right] \leq \delta(C_1 + C_2)(C(p,t+\delta) + 1). $$

As $\delta$ can be made arbitrary small and the upper bound is uniform in $N$ and $\theta$, we get (19) by the Markov inequality. With the same reasoning, we get that

$$ \mathbb{E} \left[ |\langle Z^N \rangle_{\tau_N+\theta} - \langle Z^N \rangle_{\tau_N}| \right] \leq \frac{\delta}{N} \times (C(p,t+\delta)+1) \sum_{k \in [N_c]} 2q_k^{\infty} |f|^2_{\infty} $$

and we deduce (20) which ends the proof of the lemma. $\square$

From lemma 3.5 and the Prokhorov theorem, we deduce that we can find a subsequence $(\bar{\nu}^{N(n)})_n$ of $(\bar{\nu}^N)_{N \in \mathbb{N}^*}$ which converges weakly in $\mathcal{D}([0,T], \mathcal{M}_{F^+}(\mathcal{X}))$ to a limit denoted $\tilde{\nu}$. Indeed, the following stronger result hold:

**Lemma 3.6.** There exist a subsequence $(\bar{\nu}^{N(n)}, c^{N(n)})$ such that $c^{N(n)}$ is the solution of (18) associated to $\bar{\nu}^{N(n)}$ and $(\bar{\nu}^{N(n)}, c^{N(n)})$ converges in $\mathcal{D}([0,T], \mathcal{M}_{F^+}(\mathcal{X})) \times \mathcal{C}([0,T], H^1(\mathcal{L}_a)^{N_m})$ to an accumulation point $(\tilde{\nu}, \tilde{c})$. Moreover, $\tilde{c}$ is a solution of (18) associated to $\tilde{\nu}$.

**Proof.** We define the mapping $\Phi : \mathcal{D}([0,T], \mathcal{M}_{F^+}(\mathcal{X})) \to (W_{L^2(\mathcal{Z}_a)}^{2,1})^{N_m} \cap \mathcal{C}([0,T], H^1(\mathcal{L}_a)^{N_m})$ such that $\Phi(\nu) = c$, where $c$ is the solution of the reaction-diffusion equation (7) associated to $\nu$, with fixed $c_0$ and boundary conditions satisfying hypothesis (H10)-(H13). Let us show that $\Phi$ is continuous.

Let $(\nu^N)_N$ be a sequence which converges to a given $\xi$ in $\mathcal{D}([0,T], \mathcal{M}_{F^+}(\mathcal{X}))$, let $c^N = \Phi(\nu^N)$ and $\tilde{c} = \Phi(\xi) \in \mathcal{C}([0,T], H^1(\mathcal{L}_a)^{N_m})$. As $\xi_t$ is cadlag $\sup_{[0,T]} |\xi_t|_{tv} < +\infty$, so, from proposition 3.1.2, there exists $C(\xi, T) < +\infty$ such that

$$ \sup_{t \in [0,T]} \|c^N(t) - \tilde{c}(t)\|_{H^1(\mathcal{L}_a)^{N_m}} \leq K_2 e^{K_1 C(\xi, T)T} \sum_{l=1}^{N_c} \int_0^T \|(\nu_t^{lN} - \xi_t^l) * \psi\|_{L^2(\mathcal{Z}_a)} dt. \quad (21) $$
From proposition D.4, for almost every $t \in [0, T]$, $d_{PR}(v_t^N, \xi_t) \xrightarrow{n} 0$. So, for almost every $t$, $d_{PR}(v_t^{IN}, \xi_t^l) \xrightarrow{n} 0$. Then, from proposition D.7, for almost every $t$, $\sup_{z \in \mathcal{Z}_a} |(v_t^{IN} - \xi_t^l) * \psi(z)| \xrightarrow{n} 0$ and so

$$ \|(v_t^{IN} - \xi_t^l) * \psi\|_{L^2(\mathcal{Z}_a)} \xrightarrow{n} 0. $$

From the dominated convergence theorem, we conclude that the right hand side of (21) goes to 0 when $N \to +\infty$, which shows that $\Phi$ is continuous from $\mathcal{D}([0, T], \mathcal{M}_{F^+}(\mathcal{X}))$ to $\mathcal{C}([0, T], H^1(\mathcal{Z}_a)^{N_m})$, hence we obtain the weak convergence of $\Phi(\bar{v}^{N(n)})$ to $\Phi(\bar{v})$. $\square$

We proceed with the following lemma.

**Lemma 3.7.** Any accumulation point $(\bar{v}, \tilde{c})$ of $(\bar{v}^N, c^N)$ in $\mathcal{D}([0, T], \mathcal{M}_{F^+}(\mathcal{X})) \times \mathcal{C}([0, T], H^1(\mathcal{Z}_a)^{N_m})$ is almost surely a solution of (17-18) in $\mathcal{C}([0, T], \mathcal{M}_{F^+}(\mathcal{X})) \times \mathcal{C}([0, T], H^1(\mathcal{Z}_a)^{N_m})$.

**Proof.** The proof is provided in appendix B $\square$

Finally, we tackle uniqueness of the solution of (17-18).

**Lemma 3.8.** If (17-18) has a solution in $\mathcal{C}([0, T], \mathcal{M}_{F^+}(\mathcal{X})) \times \mathcal{C}([0, T], H^1(\mathcal{Z}_a)^{N_m})$ then this solution is unique.

**Proof.** The proof is provided in appendix C $\square$

We conclude from lemma (3.7-3.8) that all accumulation points of $(\bar{v}^N, c^N)$ are equal a.s., therefore the sequence is weakly convergent and its limit is a.s. equal to the unique solution of (17-18) in $\mathcal{C}([0, T], \mathcal{M}_{F^+}(\mathcal{X})) \times \mathcal{C}([0, T], H^1(\mathcal{Z}_a)^{N_m})$. This ends the proof of theorem 3.4.

# 4. Conclusive discussion

In this study, we introduced a generic IBM of biological tissues interacting with diffusive compounds, formulated within a rigorous PDMP framework. The main achievement of this work is the proof of convergence of this stochastic IBM model under a large population assumption towards a non local interacting particle like PDE model. Although the general functional analysis framework is well established, in particular in the theoretical ecology and evolution community [10–12], our study addresses technical challenges specific to biological tissues. Notably, cell-to-cell interactions introduce nonlinearities in the deterministic limit, and complex regulatory mechanisms governing cell fate add further technical difficulties to the proofs.

Interestingly, our work establishes a link between IBMs, widely used for tissue modeling, and non local PDE models for cell movements (see [21] and references therein), demonstrating that the latter can emerge as large-population limits of the former. Although we do not provide explicit convergence rate estimates, making it difficult to provide ranges on the number of cells required for convergence in real tissues, our results emphasize the importance of considering population size when selecting a modeling approach.

Despite this limitation, our work is of special interest for tissue modelers by bridging stochastic IBM and deterministic continuous population models. This connection enables a complementary modeling strategy: micro-scale data, such as microscopy imaging, can be
used to calibrate the IBM, while the deterministic limit allows upscaling this calibration to the
continuous population model. Conversely, parameters inferred from macro-scale data, such
as metabolomics, can be downscaled to the IBM. Similarly, a stochastic IBM can be used to
model accurately intra-specific variability within cell populations, but the deterministic limit
would be more appropriate for intensive numerical exploration such as sensitivity analysis.
For example, [15] employed this approach to analyze cell turnover dynamics in the intestinal
crypt in presence of microbial metabolites.

References

[1] M. Meier-Schellersheim, I.D. Fraser, and F. Klauschen, Multiscale modeling for biologists, Wiley Interdisciplinary Reviews: Systems Biology and Medicine 1 (2009), pp. 4–14.

[2] A.G. Fletcher and J.M. Osborne, *Seven challenges in the multiscale modeling of multicellular tissues*, WIREs mechanisms of disease 14 (2022), p. e1527.

[3] N. Famaey and J.V. Sloten, *Soft tissue modelling for applications in virtual surgery and surgical robotics*, Computer methods in biomechanics and biomedical engineering 11 (2008), pp. 351–366.

[4] M. Ben Amar, P. Nassoy, and L. LeGoff, *Physics of growing biological tissues: the complex cross-talk between cell activity, growth and resistance*, Philosophical Transactions of the Royal Society A 377 (2019), p. 20180070.

[5] J. Rosenbauer, C. Zhang, B. Mattes, I. Reinartz, K. Wedgwood, S. Schindler, C. Sinner, S. Scholpp, and A. Schug, *Modeling of wnt-mediated tissue patterning in vertebrate embryogenesis*, PLOS Computational Biology 16 (2020), p. e1007417.

[6] R.C. Rockne, A. Hawkins-Daarud, K.R. Swanson, J.P. Sluka, J.A. Glazier, P. Macklin, D.A. Hornuth, A.M. Jarrett, E.A. Lima, J.T. Oden, et al., *The 2019 mathematical oncology roadmap*, Physical biology 16 (2019), p. 041005.

[7] J.M. Osborne, A.G. Fletcher, J.M. Pitt-Francis, P.K. Maini, and D.J. Gavaghan, *Comparing individual-based approaches to modelling the self-organization of multicellular tissues*, PLoS computational biology 13 (2017), p. e1005387.

[8] R. Losick and C. Desplan, *Stochasticity and Cell Fate*, Science 320 (2008), pp. 65–68, Available at https://www.science.org/doi/10.1126/science.1147888.

[9] S.N. Ethier and T.G. Kurtz, *Markov Processes: Characterization and Convergence*, Wiley Series in Probability and Its Applications, Wiley, 1986.

[10] N. Fournier and S. Méléard, *A microscopic probabilistic description of a locally regulated population and macroscopic approximations*, Annals of Applied Probability 14 (2004), pp. 1880–1919, Available at https://doi.org/10.1214/105051604000000882.

[11] N. Champagnat, R. Ferrière, and S. Méléard, *Individual-based probabilistic models of adaptive evolution and various scaling approximations*, in *Seminar on stochastic analysis, random fields and applications V*. Springer, 2007, pp. 75–113.

[12] V.C. Tran, *Stochastic particle models for problems of adaptive evolution and for the approximations of statistical solutions*, Theses, Université de Nanterre - Paris X, 2006, Available at https://tel.archives-ouvertes.fr/tel-00125100.

[13] F. Campillo and C. Fritsch, *Weak convergence of a mass-structured individual-based model*, Applied Mathematics & Optimization 72 (2014), p. 37–73.

[14] P. Michel, *Large population and size scale limit of a stochastic particle model*, Mathematical Models and Methods in Applied Sciences 27 (2017), pp. 581–615, Available at https://doi.org/10.1142/S0218202517500105.

[15] L. Darrigade, M. Haghebaert, C. Cherbuy, S. Labarthe, and B. Laroche, *A PDMP model of the epithelial cell turn-over in the intestinal crypt including microbiota-derived regulations*, Journal of Mathematical Biology 84 (2022), pp. 1–67.

[16] G.T. Eisenhoffer, P.D. Loftus, M. Yoshigi, H. Otsuna, C.B. Chien, P.A. Morcos, and J. Rosenblatt, *Crowding induces live cell extrusion to maintain homeostatic cell numbers in epithelia*, Nature
484 (2012), pp. 546–549.

[17] N. Ikeda and S. Watanabe, Stochastic differential equations and diffusion processes, North-Holland Publishing Company, 1989.

[18] G. Crippa and M. Lécureux-Mercier, Existence and uniqueness of measure solutions for a system of continuity equations with non-local flow, Nonlinear Differential Equations and Applications NoDEA 20 (2013), pp. 523–537.

[19] S.R. Coppoletta, A criterion of convergence of measure-valued processes: application to measure branching processes, Stochastics 17 (1986), pp. 43–65, Available at https://doi.org/10.1080/17442508608833382.

[20] A. Joffe and M. Metivier, Weak convergence of sequences of semimartingales with applications to multitype branching processes, Advances in Applied Probability 18 (1986), pp. 20–65, Available at http://www.jstor.org/stable/1427238.

[21] L. Chen, K. Painter, C. Surulescu, and A. Zhigun, Mathematical models for cell migration: a non-local perspective, Philosophical Transactions of the Royal Society B 375 (2020), p. 20190379.

[22] R.F. Curtain and H. Zwart, An introduction to infinite-dimensional linear systems theory, Vol. 21, Springer Science & Business Media, 2012.

[23] H. Amann, Nonhomogeneous linear and quasilinear elliptic and parabolic boundary value problems, in Function spaces, differential operators and nonlinear analysis, Springer, 1993, pp. 9–126.

[24] P. Billingsley, Convergence of probability measures, Wiley, New York, 1999.

[25] D.J. Daley and D. Vere-Jones, An Introduction to the Theory of Point Processes, Probability and Its Applications, Springer, 2003, second edition, Available at https://books.google.fr/books?id=GqRXYFxe0l0C.

[26] A. Pazy, Semigroups of linear operators and applications to partial differential equations, Vol. 44, Springer Science & Business Media, 2012.

[27] J.L. Lions and E. Magenes, Non-homogeneous boundary value problems and applications, Die Grundlehren der mathematischen Wissenschaften Vol. 1, Berlin: Springer-Verlag, 1972.

[28] I. Lasiecka, Unified theory for abstract parabolic boundary problems–a semigroup approach, Applied Mathematics and Optimization 6 (1980), pp. 287–333.

[29] B. Perthame, Transport Equations in Biology, Birkhäuser Basel, Basel, 2007.

[30] L. Darrigade, Modélisation du dialogue hôte-microbiote au voisinage de l'épithélium de l'intestin distal, Ph.D. diss., Ecole Doctorale de Mathématiques Hadamard, Université Paris-Saclay, 2020, Available at http://www.theses.fr/2020UPASM0008.

## Appendix A. Proof of proposition 3.1

### Assertion 3.1.1. Let us define the function $G$ on $[0, T] \times L^2(\mathcal{Z}_a)^{N_m}$ as

$$G(t, c)(z) := \left( \sum_{l=1}^{N_c} \gamma_{1l}(c(z)) v_l^l * \psi(z), \dots, \sum_{l=1}^{N_c} \gamma_{N_m l}(c(z)) v_l^l * \psi(z) \right) \in L^\infty([0, T], L^2(\mathcal{Z}_a)^{N_m})$$

and show that $G$ is Lipschitz with respect to $c$ on $L^2(\mathcal{Z}_a)^{N_m}$, uniformly for $t \in [0, T]$. $\gamma_{ml}$ are Lipschitz under hypothesis (H11), let $K(\gamma) := \sup_{m,l} |\gamma_{ml}|_{\text{Lip}}$. Let $c_1, c_2$ in $L^2(\mathcal{Z}_a)^{N_m}$, then

$$ \begin{aligned} \sup_{t \in [0,T]} \|G(t, c_1) - G(t, c_2)\|_{L^2(\mathcal{Z}_a)^{N_m}} &\le \sup_{t \in [0,T]} |v_t|_{\text{tv}} \times |\psi|_{\infty} \sum_{m=1}^{N_m} \sum_{l=1}^{N_c} \|\gamma_{ml}(c_1) - \gamma_{ml}(c_2)\|_{L^2(\mathcal{Z}_a)} \\ &\le \sup_{t \in [0,T]} |v_t|_{\text{tv}} \times |\psi|_{\infty} N_m N_c K(\gamma) \|c_1 - c_2\|_{L^2(\mathcal{Z}_a)^{N_m}}. \end{aligned} $$
With the notations introduced in appendix E, we denote $S_m$ the analytic semi-group $S$ on $L^2(\mathcal{Z}_a)$ generated by $\sigma\Delta$ with boundary conditions of chemical $m$, $D_m$ the related lifting operator ($D$ in appendix E) and $\lambda_m$ the associated scalar ($\lambda$ in appendix E). With a standard fixed point argument (see [22, 23]) on the system

$$c_m(z,t) = S_m(t)c_{0m}(z) + (\lambda_m I - \sigma\Delta)(\int_0^t S_m(t-s)D_m(f_m)(z,s)ds) \quad (A1) \\ + \int_0^t S_m(t-s)\sum_{l=1}^{N_c} \gamma_{ml}(c(s)) v_s^l * \psi(z)ds$$

and the semi-group of $L^2(\mathcal{Z}_a)^{N_m}$ defined by $S = \text{diag}(S_1, \dots, S_{N_m})$, we prove the existence of a unique weak solution $c$ in $\mathcal{C}([0, T], L^2(\mathcal{Z}_a)^{N_m})$, such that for $1 \le m \le N_m$, $c_m$ satisfies (7). Then

$$u(z,t) = \sum_{l=1}^{N_c} \gamma_{ml}(c(z,t)) v_t^l * \psi(z) \quad \in L^\infty([0,T], L^2(\mathcal{Z}_a)) \subset W_{L^2(\mathcal{Z}_a)}^{0,0}.$$

Indeed, $\gamma_{ml}$ are bounded lipschitz on $\mathbb{R}_+^{N_m}$ under hypothesis (H11), so $\gamma_{ml}(c(z,t)) \in \mathcal{C}([0,T], L^2(\mathcal{Z}_a))$ as $c \in \mathcal{C}([0,T], L^2(\mathcal{Z}_a)^{N_m})$ according to the previous step. Moreover, $v_t^l * \psi(z) \in L^\infty([0,T], L^2(\mathcal{Z}_a))$ because $\psi$ is bounded and $\sup_{t \in [0,T]} |v_t|_{tv} < \infty$.

We obtain the adequate regularity result by observing that $c_m$ is a solution of a problem of type (E4), with $c_{0m} \in H^1(\mathcal{Z}_a)$, $f_m$ verifying the regularity conditions stated in proposition E.1 and $u \in W_{L^2(\mathcal{Z}_a)}^{0,0}$. Therefore $c_m \in W_{L^2(\mathcal{Z}_a)}^{2,1} \cap \mathcal{C}([0,T], H^1(\mathcal{Z}_a))$.

**Assertion 3.1.2.** Let $v^1, v^2$ be two measurable functions from $[0,T]$ to $\mathcal{M}_{F^+}(\mathcal{X})$ and $c^1, c^2$ be the associated solutions of (7). By the continuous injection of $(W_{L^2(\mathcal{Z}_a)}^{2,1})^{N_m}$ in $\mathcal{C}([0,T], H^1(\mathcal{Z}_a)^{N_m})$, we have

$$\sup_{t \in [0,T]} \|c^1(t) - c^2(t)\|_{H^1(\mathcal{Z}_a)^{N_m}} \le \beta_1 \|c^1 - c^2\| \left(\mathcal{W}_{L^2(\mathcal{Z}_a)}^{2,1}\right)^{N_m} \quad (A2)$$

and applying proposition E.1 with $u = G_1 - G_2$, we obtain

$$\|c^1 - c^2\| \left(\mathcal{W}_{L^2(\mathcal{Z}_a)}^{2,1}\right)^{N_m} \le \beta_2 \|G_1 - G_2\| \left(\mathcal{W}_{L^2(\mathcal{Z}_a)}^{0,0}\right)^{N_m}. \quad (A3)$$

We can now bound $\|G_1 - G_2\| \left(\mathcal{W}_{L^2(\mathcal{Z}_a)}^{0,0}\right)^{N_m}$ as

$$\begin{aligned} \|G_1 - G_2\| \left(\mathcal{W}_{L^2(\mathcal{Z}_a)}^{0,0}\right)^{N_m} &\le \sum_{m=1}^{N_m} \left( \int_0^T \left\| \sum_{l=1}^{N_c} \left( \gamma_{ml}(c^1) v_t^{1l} * \psi - \gamma_{ml}(c^2) v_t^{2l} * \psi \right) \right\|_{L^2(\mathcal{Z}_a)}^2 dt \right)^{1/2} \\ &\le N_m \left[ |\gamma|_{\infty} \sum_{l=1}^{N_c} \int_0^T \|(v_t^{1l} - v_t^{2l}) * \psi\|_{L^2(\mathcal{Z}_a)} dt \right. \\ &\qquad \left. + C |\psi|_{\infty} K(\gamma) \int_0^T \sup_{s \in [0,T]} \|c^1(s) - c^2(s)\|_{H^1(\mathcal{Z}_a)^{N_m}} dt \right]. \quad (A4) \end{aligned}$$
Combining (A2), (A3), (A4) and applying Gronwall lemma, we obtain the desired estimate with $K_1 := N_m\beta_1\beta_2 |\psi|_\infty K(\gamma)$ and $K_2 := N_m\beta_1\beta_2 |\gamma|_\infty$.

## Appendix B. Proof of lemma 3.7

For a fixed $t \in [0, T]$ and $f \in \mathcal{C}^{1,1,\bullet}([0, T] \times \mathcal{X}, \mathbb{R})$, let us define the following mapping on $\mathcal{D}([0, T], \mathcal{M}_{F^+}(\mathcal{X}))$ as

$$ \Psi_t(v) = \langle v_t, f_t \rangle - \langle v_0, f_0 \rangle - \int_0^t \langle v_s, \mathcal{F}(\cdot, v_s * F(\cdot)) \nabla_z f_s + \partial_s f_s \rangle ds \\ - \sum_{k=1}^{N_c} \int_0^t \int_{\mathcal{X}} \langle \mu_k(x), f_s \rangle q_k(x, v_s * D(z), \Phi(v)_s * \psi(z)) v_s(dx) ds, $$

and establish two technical lemmas.

**Lemma B.1.** Consider a sequence $(v^n)_{n \in \mathbb{N}}$ of $\mathcal{D}([0, T], \mathcal{M}_{F^+}(\mathcal{X}))$ that converges towards a limit $v$ with $v \in \mathcal{C}([0, T], \mathcal{M}_{F^+}(\mathcal{X}))$, then $\lim_{n \to \infty} \Psi_t(v^n) = \Psi_t(v)$.

*Proof.* Let $C(T, v) = \sup_{[0,T]} \langle v_t, 1 \rangle < +\infty$ as $v_t$ is cadlag. The convergence of $v^n$ to $v$ is uniform in time (prop. D.3). From the definition of the Prokhorov distance (D1), we deduce that for all $\varepsilon$ and $n$ big enough, $\sup_{t \in [0,T]} \langle v_t^n, 1 \rangle \le \sup_{t \in [0,T]} \langle v_t, 1 \rangle + \varepsilon \le C(T, v) + \varepsilon$. Therefore, we restrict ourselves without loss of generality to sequences such that, for all $n \in \mathbb{N}$, $\sup_{t \in [0,T]} \langle v_t^n, 1 \rangle \le 2C(T, v)$.

The uniform convergence of $(v^n)_{n \in \mathbb{N}}$ to $v$ implies that for all $t \le T$, the projection at $t$ $\pi_t: v \mapsto v_t$ is continuous at every element $v$ continuous in time. Then, by definition of the weak topology, we see that the mapping $v \mapsto \langle v_t, f \rangle$ is continuous for any $f \in \mathcal{C}_b(\mathcal{X}, \mathbb{R})$ at each element $v$ of $\mathcal{C}([0, T], \mathcal{M}_{F^+}(\mathcal{X}))$.

We can now show the continuity of $\Psi_t$ at $v$. We denote $c^n = \Phi(v^n)$ and $c = \Phi(v)$, then

$$ \Psi_t(v^n) - \Psi_t(v) = \langle v_t^n - v_t, f_t \rangle - \langle v_0^n - v_0, f_0 \rangle - \int_0^t \langle v_s^n - v_s, \partial_s f_s \rangle ds \\ \underbrace{\qquad - \left( \int_0^t \langle v_s^n, \mathcal{F}(\cdot, v_s^n * F(\cdot)) \nabla_z f_s \rangle ds - \int_0^t \langle v_s, \mathcal{F}(\cdot, v_s * F(\cdot)) \nabla_z f_s \rangle ds \right)}_{(B_n)} \\ \underbrace{\qquad - \sum_{k=1}^{N_c} \int_0^t \left( \int_{\mathcal{X}} \langle \mu_k(x), f_s \rangle q_k(x, v_s^n * D(z), c_s^n * \psi(z)) v_s^n(dx) ds - \int_{\mathcal{X}} \langle \mu_k(x), f_s \rangle q_k(x, v_s * D(z), c_s * \psi(z)) v_s(dx) ds \right) ds}_{(C_n)} $$

According to the foregoing, regularity of $f$ and compactness of $\mathcal{X}$, we have $\langle v_t^n - v_t, f_t \rangle \underset{n \to +\infty}{\to} 0$, $\langle v_0^n - v_0, f_0 \rangle \underset{n \to +\infty}{\to} 0$ and $s \mapsto \langle v_s^n - v_s, \partial_s f_s \rangle$ converges pointwise to 0.

By dominated convergence theorem, $\int_0^t \langle v_s^n - v_s, \partial_s f_s \rangle ds \underset{n \to +\infty}{\to} 0$. Therefore, $\lim_{n \to +\infty} A_n = 0$.
We have $B_n = B_n^1 + B_n^2$ with

$$B_n^1 = \int_0^t \langle v_s^n - v_s, \mathcal{F}(\cdot, v_s * F(\cdot)) \nabla_z f_s \rangle ds,$$

$$B_n^2 = \int_0^t \langle v_s^n, (\mathcal{F}(\cdot, v_s^n * F(\cdot)) - \mathcal{F}(\cdot, v_s * F(\cdot))) \nabla_z f_s \rangle ds.$$

$\mathcal{F}(\cdot, v_s * F(\cdot))\nabla_z f_s$ is continuous bounded as $\mathcal{X}$ is compact hence, by dominated convergence, $B_n^1 \underset{n \to +\infty}{\longrightarrow} 0$. Moreover, $\mathcal{F}$ is Lipschitz on $\mathcal{X}$ compact, therefore

$$B_n^2 \le \int_0^t \langle v_s^n, |\mathcal{F}|_{\text{Lip}} \times |v_s^n * F(\cdot) - v_s * F(\cdot)| \times |\nabla_z f_s| \rangle ds.$$

As $F$ is Lipschitz and bounded on $\hat{\mathcal{X}}$ (H8), we can apply proposition D.7 : $\sup_{z \in \mathcal{X}} |(v_t^n - v_t) * F(z)| \underset{n \to +\infty}{\longrightarrow} 0$, which proves that $B_n^2 \underset{n \to +\infty}{\longrightarrow} 0$.

Finally, let $C_n = C_n^1 + C_n^2$ with

$$C_n^1 = \sum_{k=1}^{N_e} \int_0^t \int_{\mathcal{X}} \langle \mu_k(x), f_s \rangle q_k(x, v_s^n * D(z), c_s^n * \psi(z)) (v_s^n(dx) - v_s(dx)) ds,$$

$$C_n^2 = \sum_{k=1}^{N_e} \int_0^t \int_{\mathcal{X}} \langle \mu_k(x), f_s \rangle (q_k(x, v_s * D(z), c_s * \psi(z)) - q_k(x, v_s^n * D(z), c_s^n * \psi(z))) v_s(dx) ds.$$

$q_k$ are continuous bounded (H4) so $C_n^1 \underset{n \to +\infty}{\longrightarrow} 0$. Regarding $C_n^2$, for all $t \in [0, T]$,

$$\sup_{z \in \mathcal{X}} |(c_t^n - c_t) * \psi(z)| \le \|\psi\|_{L^2} \|c_t^n - c_t\|_{L^2(\mathcal{X}_a)^{N_m}} \\ \le \|\psi\|_{L^2} \sup_{t \in [0, T]} \|\Phi(v_t^n) - \Phi(v)_t\|_{H^1(\mathcal{X}_a)^{N_m}} \underset{n \to +\infty}{\longrightarrow} 0$$

where we used continuity of $\Phi$ obtained in the proof of lemma 3.6. From regularity of $\psi$ (H7) and proposition D.7, $\sup_{z \in \mathcal{X}} |(v_t^n - v_t) * D(z)| \underset{n \to +\infty}{\longrightarrow} 0$. Finally, using that $q_k$ is Lipschitz with respect to its second and third variable (H5) and dominated convergence, we get that $C_n^2 \underset{n \to +\infty}{\longrightarrow} 0$. This ends the proof of the continuity of $\Psi_t$ at $v$. $\square$

Let $(\tilde{v}, \tilde{c})$ be an accumulation point of $(\tilde{v}^N, c^N)_{N \in \mathbb{N}^*}$ and $(\tilde{v}^{N(n)}, c^{N(n)})_{n \in \mathbb{N}}$ be a subsequence which converges to $(\tilde{v}, \tilde{c})$. We recall that from lemma 3.6, $\tilde{c}$ is the solution of (7) associated with $\tilde{v}$. Then the following holds true.

**Lemma B.2.** $(\tilde{v}, \tilde{c})$ is almost surely a solution of (17-18).

**Proof.** We observe that $\tilde{v}$ is solution of (17) if and only if $\Psi_t(\tilde{v}) = 0$ a.s. for all $t$ and $f$. We therefore prove the latter. We first show that

$$\tilde{v} \in \mathcal{C}([0, T], \mathcal{M}_{F^+}(\mathcal{X})) \text{ a. s.} \qquad (B1)$$

$\tilde{v}$ is right continuous by definition. Let us prove the left continuity. Considering the total variation norm on $\mathcal{M}_{F^+}(\mathcal{X})$ and its characterization given in proposition D.1, let $f \in \mathcal{C}(\mathcal{X})$ with $|f|_{\infty} \le 1$. For all $n$, random jumps of $\tilde{v}^{N(n)}$ happen one at a time a. s. so $|\langle \tilde{v}_t^{N(n)}, f \rangle|$ -
$$ \langle \bar{\mathbf{v}}_{t-}^{N(n)}, f \rangle \le 2 \frac{|f|_{\infty}}{N(n)} \le \frac{2}{N(n)}, \text{ and } \sup_{0 \le t \le T} |\bar{\mathbf{v}}_{t}^{N(n)} - \bar{\mathbf{v}}_{t-}^{N(n)}|_{\text{tv}} \le \frac{2}{N(n)} \text{ a. s.} $$

From proposition D.5, we conclude that $\lim_{n \to \infty} \sup_{0 \le t \le T} d_{PR}(\bar{\mathbf{v}}_t^{N(n)}, \bar{\mathbf{v}}_{t-}^{N(n)}) = 0$ a. s. and deduce (B1) (see [9], theorem 10.2).

We now observe that $\Psi_t(\bar{\mathbf{v}}^N) = Z_t^{f,N} - R_t^{f,N}$ with $Z_t^{f,N}$ the martingale defined in proposition 3.3 and $R_t^{f,N} := \int_0^t \langle \bar{\mathbf{v}}_s^N, (\mathcal{F}(\cdot, \bar{\mathbf{v}}_s^N * F(\cdot) - F(0)/N) - \mathcal{F}(\cdot, \bar{\mathbf{v}}_s^N * F(\cdot))) \nabla_z f_s \rangle ds$. Then, under (H16):

$$ \begin{aligned} \mathbb{E} \left[ \left( \Psi_t(\bar{\mathbf{v}}^{N(n)}) \right)^2 \right] &\le 2 \mathbb{E} \left[ \left\langle Z^{f,N(n)} \right\rangle_t \right] + 2 \mathbb{E} \left[ \left( R_t^{f,N(n)} \right)^2 \right] \\ &\le \frac{2}{N(n)} (C(p,t) + 1) t |f|_{\infty}^2 \sum_{k=1}^{N_e} q_k^{\infty} + \frac{2}{N(n)^2} (C(p,t) + 1) t \left( \sup_{s \le t} |K_{\mathcal{F}}|_{\infty} |F(0)| |\nabla_z f_s|_{\infty} \right)^2 \\ &\xrightarrow{n \to \infty} 0 \end{aligned} $$

with $C(p,t)$ a constant defined in Prop. 3.2. Therefore, $(\Psi_t(\bar{\mathbf{v}}^{N(n)}))_n$ converges in $L^2$, and so in $L^1$ by Cauchy-Schwartz, to 0.

We know that $\bar{\mathbf{v}}^{N(n)}$ converges weakly to $\tilde{\mathbf{v}}$ and $\Psi_t$ is continuous at $\tilde{\mathbf{v}}$ a.s. after (B1) and lemma B.1. We deduce from the continuous mapping theorem that $|\Psi_t(\bar{\mathbf{v}}^{N(n)})|$ converges weakly to $|\Psi_t(\tilde{\mathbf{v}})|$. Moreover, under hyp. (H16), from prop. 3.2.1 and using hypotheses on boundedness and regularity of the different functions involved in the definition of $\Psi_t$, it is straightforward to show that the sequence $(\Psi_t(\bar{\mathbf{v}}^{N(n)}))_n$ is bounded in $L^p$ for some $p \ge 2$ and so it is uniformly integrable. Then, we have (see [24], theorem 3.5) that, $\forall q \le p$, $\lim_{n \to \infty} \mathbb{E}[|\Psi_t(\bar{\mathbf{v}}^{N(n)})|^q] = \mathbb{E}[|\Psi_t(\tilde{\mathbf{v}})|^q]$. In particular, $0 = \lim_{n \to \infty} \mathbb{E}[|\Psi_t(\bar{\mathbf{v}}^{N(n)})|] = \mathbb{E}[|\Psi_t(\tilde{\mathbf{v}})|]$, which proves that $\Psi_t(\tilde{\mathbf{v}}) = 0$ a. s.

Note that $\tilde{\mathbf{v}}$ is continuous in time a.s., and therefore $\bar{\mathbf{v}}_0^N$ converges weakly to $\tilde{\mathbf{v}}_0$. But, under hyp. (H15), $\bar{\mathbf{v}}_0^N$ converges in distribution to $\xi_0$, so $\xi_0$ is indeed the initial condition for $\tilde{\mathbf{v}}$.

This is true for any $t \in [0, T]$ and test functions $f$, we deduce that $\tilde{\mathbf{v}}$ is almost surely a solution of (17), and the sequence of concentrations $(c^{N(n)})_n$ converges as well to the solution of (18). $\square$

## Appendix C. Proof of lemma 3.8

We proceed in two steps. In a first step, we prove that on a finite time horizon, positive solutions of (17) have finite mass. Let $(c_0, \xi_0) \in H^1(\mathcal{Z}_a)^{N_m} \times \mathcal{M}_{F^+}(\mathcal{X})$ be the initial condition. We suppose that there exists a solution $(c_t, \xi_t)_{t \ge 0} \in \mathcal{C}([0, T], H^1(\mathcal{Z}_a)^{N_m}) \times \mathcal{D}([0, T], \mathcal{M}_{F^+}(\mathcal{X}))$. Let us show that the mass of $\xi_t$ is bounded at finite time by a constant depending only of $\xi_0$ and $t$. From (17), for $t \le T$,

$$ \begin{aligned} \langle \xi_t, 1 \rangle &= \langle \xi_0, 1 \rangle + \sum_{k=1}^{N_e} \int_0^t \sum_{l=1}^{N_c} \int_{\mathcal{Z}} \langle \mu_k(z, l), 1 \rangle q_k^{\infty} q_k((z, l), \xi_s * D(z), c_s * \psi_a(z)) \xi_s(dz, dl) ds \\ &\le \langle \xi_0, 1 \rangle + \sum_{k=1}^{N_e} q_k^{\infty} \int_0^t \sum_{l=1}^{N_c} \int_{\mathcal{Z}} \xi_s(dz, dl) ds = \langle \xi_0, 1 \rangle + \sum_{k=1}^{N_e} q_k^{\infty} \int_0^t \langle \xi_s, 1 \rangle ds. \quad (C1) \end{aligned} $$

By Gronwall lemma, we conclude that $\langle \xi_t, 1 \rangle \le \langle \xi_0, 1 \rangle \exp(t \sum_{k=1}^{N_e} q_k^{\infty}) < +\infty$.
In a second step, we prove that a positive solution of (17-18) is unique. Let $(c_t^1, \xi_t^1)$ and $(c_t^2, \xi_t^2)$ be two solutions of (17-18) on $[0, T]$, with $\xi_t^i$ positive as in the previous step. We use the flat distance (or Kantorovitch-Rubinstein distance) $d_F$ on $\mathcal{M}_{F^+}(\mathcal{X})$ defined as

$$d_F(v, \xi) = \sup_{f \in \mathcal{L}ip_{1,1}(\mathcal{X}, \mathbb{R})} \langle v - \xi, f \rangle. \qquad (C2)$$

Let us prove that $||c_t^1 - c_t^2||_{L^2(\mathcal{X})^N_m} + d_F(\xi_t^1, \xi_t^2) = 0$.

From step 1, we can define $C(T, \xi^1 + \xi^2) := \sup_{s \le T} \langle \xi_s^1 + \xi_s^2, 1 \rangle < +\infty$. Moreover, $\xi_t^1, \xi_t^2$ are finite positive measures and, working from prop. 3.1.2 and setting $C := |\psi|_{\vee} \text{Leb}(\mathcal{X})^{1/2} K_1 e^{K_2 C(T, \xi_1 + \xi_2)T}$, we obtain that

$$||c_t^1 - c_t^2||_{L^2(\mathcal{X})^N_m} \le C \int_0^t d_F(\xi_s^1, \xi_s^2) ds. \qquad (C3)$$

We now bound $d_F(\xi_t^1, \xi_t^2)$ using formula (17). As done in ref. [12], the idea is to cancel in equation (17) the integral term accounting for the deterministic movement of cells between jumps. Let $f \in \mathcal{L}ip_{1,1} \cap \mathcal{C}^{1,\bullet}(\mathcal{X}, \mathbb{R})$. From lemma F.1, there exists a unique function $u_1$ from $[0, T] \times \mathcal{X}$ to $\mathbb{R}$ such that, $\forall l \in [N_c]$, $u_s^1(z, l)$ is solution of

$$\partial_s u_s^1(z, l) + \mathcal{F}(z, \xi_s^1 * F(z)) \nabla_z u_s^1(z, l) = 0 \quad \text{and} \quad u_t^1(x) = f(x). \qquad (C4)$$

We know that $u_s^1(z, l) = f(Z^1(s; t, z), l)$ where $Z^1(\cdot)$ are the characteristics associated to (F1) when $\xi_s = \xi_s^1$. Similarly, let $u_s^2(z, l) = f(Z^2(s; t, z), l)$, with $\xi_s = \xi_s^2$ in (F1). Moreover, still from lemma F.1, $u^i(\cdot, l) \in \mathcal{C}^1([0, T] \times \mathcal{X}, \mathbb{R})$.
We then apply formula (17), with $\xi_0^1 = \xi_0^2 = \xi_0$.

$$
\begin{align*}
|\langle \xi_t^1 - \xi_t^2, u \rangle| &= |\langle \xi_0, u_0^1 - u_0^2 \rangle \\
&\quad + \sum_{k=1}^{N_e} \int_0^t \sum_{l=1}^{N_c} \int_{\mathcal{Z}} \langle \mu_k(z,l), u_s^1 \rangle q_k((z,l), \xi_s^1 * D(z), c_s^1 * \psi(z)) \xi_s^1 (dz, dl) ds \\
&\quad - \sum_{k=1}^{N_e} \int_0^t \sum_{l=1}^{N_c} \int_{\mathcal{Z}} \langle \mu_k(z,l), u_s^2 \rangle q_k((z,l), \xi_s^2 * D(z), c_s^2 * \psi(z)) \xi_s^2 (dz, dl) ds \\
&= |\langle \xi_0, u_0^1 - u_0^2 \rangle| \\
&\quad + \underbrace{\sum_{k=1}^{N_e} \int_0^t \sum_{l=1}^{N_c} \int_{\mathcal{Z}} \langle \mu_k(z,l), u_s^1 - u_s^2 \rangle q_k((z,l), \xi_s^1 * D(z), c_s^1 * \psi(z)) \xi_s^1 (dz, dl) ds}_{(B)} \\
&\quad + \underbrace{\sum_{k=1}^{N_e} \int_0^t \sum_{l=1}^{N_c} \int_{\mathcal{Z}} \langle \mu_k(z,l), u_s^2 \rangle (q_k((z,l), \xi_s^1 * D(z), c_s^1 * \psi(z)) - q_k((z,l), \xi_s^2 * D(z), c_s^1 * \psi(z))) \xi_s^1 (dz, dl) ds}_{(C)} \\
&\quad + \underbrace{\sum_{k=1}^{N_e} \int_0^t \sum_{l=1}^{N_c} \int_{\mathcal{Z}} \langle \mu_k(z,l), u_s^2 \rangle (q_k((z,l), \xi_s^2 * D(z), c_s^1 * \psi(z)) - q_k((z,l), \xi_s^2 * D(z), c_s^2 * \psi(z))) \xi_s^1 (dz, dl) ds}_{(D)} \\
&\quad + \underbrace{\sum_{k=1}^{N_e} \int_0^t \sum_{l=1}^{N_c} \int_{\mathcal{Z}} \langle \mu_k(z,l), u_s^2 \rangle q_k((z,l), \xi_s^2 * D(z), c_s^2 * \psi(z)) (\xi_s^1 - \xi_s^2) (dz, dl) ds|}_{(E)}
\end{align*}
$$

Before proceeding, we recall the following simple properties that we will use extensively.
Let $f$ be Lipschitz on $\mathcal{X}$, then $|\frac{f}{|f|_\infty}| \le 1$, $|\frac{f}{|f|_\infty}|_{\text{Lip}} \le 1$ and, for $v_1, v_2 \in \mathcal{M}_{F^+}(\mathcal{X})$

$$ \langle v_1 - v_2, f \rangle \le |f|_\vee \times d_F(v_1, v_2). $$

We first derive a bound for (C). Remember that $\langle \mu_k(z, l), u_s^2 \rangle$ is the sum of at most two punctual evaluations of $u_s^2$. So $|\langle \mu_k(z, l), u_s^2 \rangle|_\infty \le 2$. We supposed (hyp. (H5)) that $q_k$ is Lipschitz, therefore there exists a constant $K(q_k)$ such that

$$ (C) \le 2 \sum_{k=1}^{N_e} \int_{0}^{t} \sum_{l=1}^{N_c} \int_{\mathcal{Z}} K(q_k) |\xi_s^1 * D(z) - \xi_s^2 * D(z)| \xi_s^1 (dz, dl) ds. $$

Under hyp. (H6), $|D_{|\mathcal{Z}|}|_\vee < \infty$ on the compact set $\mathcal{Z} := \{z-y : (z,y) \in \mathcal{Z}^2\}$. Therefore, for all $z \in \mathcal{Z}$, $|\xi_s^1 * D(z) - \xi_s^2 * D(z)| \le |D_{|\mathcal{Z}|}|_\vee d_F(\xi_s^1, \xi_s^2)$. Let $C_{(C)} := 2 \sum_{k=1}^{N_e} K(q_k)C(t, \xi_1 + \xi_2) |D_{|\mathcal{Z}|}|_\vee$. Then

$$ (C) \le C_{(C)} \int_{0}^{t} d_{F}(\xi_{s}^{1}, \xi_{s}^{2}) ds. $$
For (D):

$$ (D) \leq 2 \sum_{k=1}^{N_e} K(q_k) \int_{0}^{t} \sum_{l=1}^{N_c} \int_{\mathcal{Z}} |(c_s^1 - c_s^2) * \psi(z)| \xi_s^1 (dz, dl) ds \\ \leq C_{(D)} \int_{0}^{t} ||c_s^1 - c_s^2||_{L^2(\mathcal{Z}_a)^{N_m}} ds $$

with $C_{(D)} := 2 \sum_{k=1}^{N_e} K(q_k) ||\psi||_{L^2(\mathbb{R}^d)} C(t, \xi_1 + \xi_2)$, using that $q_k$ is Lipschitz and $|(c_s^1 - c_s^2) * \psi(z)| \leq ||c_s^1 - c_s^2||_{L^2(\mathcal{Z}_a)^{N_m}} ||\psi||_{L^2(\mathbb{R}^d)}$.

We move on to term (E). For $s \leq t$, let $f_s: (z, l) \mapsto \langle \mu_k(z, l), u_s^2 \rangle q_k(z, l, \xi_s^2 * D(z), c_s^2 * \psi(z))$ be defined on $\mathcal{X}$. Let us show that $f_s$ is uniformly bounded and uniformly Lipschitz for $s \leq t$, and therefore that we have

$$ (E) \leq C_{(E)} \int_{0}^{t} d_{F}(\xi_{s}^{1}, \xi_{s}^{2}) ds, $$

with $C_{(E)} := \sup_{s \leq t} |f_s|_{\vee} < +\infty$. On the one hand, $q_k$ and $u_s^2$ are bounded, so $\sup_{s \leq t} |f_s|_{\infty} < +\infty$.

On the other hand, $\forall l \in [N_c]$, $z \mapsto q_k(z, l, \xi_s^2 * D(z), c_s^2 * \psi(z))$ is Lipschitz uniformly for $s \leq t$. Indeed, $z \mapsto \xi_s^2 * D(z)$ is Lipschitz uniformly for $s \leq t$ under hyp. (H6). It is the case of $c_s^2 * \psi(z)$ as well because $\sup_{s \leq t} ||c_s^2||_{L^2(\mathcal{Z}_a)^{N_m}} < +\infty$ from Prop. 3.1 and $\psi(z)$ is Lipschitz on $\mathcal{Z}$ under hyp. (H7). Eventually, $\forall l \in [N_c]$, $q_k(\cdot, l, \dots)$ is Lipschitz under hyp. (H5) so $z \mapsto q_k((z, l), \xi_s^2 * D(z), c_s^2 * \psi(z))$ is Lipschitz uniformly for $s \leq t$.

Moreover, for all $l$, $z \mapsto u_s^2(z, l)$ is Lipschitz uniformly in time as it is continuously derivable, so $u_s^2$ is Lipschitz on $\mathcal{X}$ uniformly in time. Then, when the jump type $k$ is a differentiation or an extrusion, $\langle \mu_k(z, l), u_s^2 \rangle$ is also Lipschitz uniformly in time. If $k$ is a division, it is the case as well because $\lambda(z)$ is Lipschitz according to hyp. (H2). Eventually, the product of two Lipschitz bounded functions is Lipschitz. Therefore, $f_s$ is Lipschitz uniformly for $s \leq t$.

In order to bound the remaining terms (A) and (B), we want an upper bound on $|u_s^1(z) - u_s^2(z)|$ depending of $d_F(\xi_s^1, \xi_s^2)$. We have

$$ |u_s^1(z) - u_s^2(z)| \leq |Z^1(t; s, z) - Z^2(t; s, z)| $$

since $u \in Lip_{1,1}(\mathcal{X}, \mathbb{R})$. To obtain this inequality is what motivates us to use of the flat metric $d_F$. Using proposition F.2, there exists a constant $C(t)$ such that $(A) \leq C_{(A)} \int_0^t d_F(\xi_s^1, \xi_s^2) ds$ with $C_{(A)} := \langle \xi_0, 1 \rangle C(t)$. We can also tackle term (B) : $(B) \leq C_{(B)} \int_0^t d_F(\xi_s^1, \xi_s^2) ds$ with $C_{(B)} := \sum_{k=1}^{N_e} 2q_k^\infty C(t, \xi^1 + \xi^2 C(t))t$.

Therefore, we showed that there exists a constant $C$ such that, $\forall f \in Lip_{1,1}(\mathcal{X}, \mathbb{R}) \cap \mathcal{C}^{1,\bullet}(\mathcal{X}, \mathbb{R})$,

$$ |\langle \xi_t^1 - \xi_t^2, f \rangle| \leq C \int_0^t ||c_s^1 - c_s^2||_{L^2(\mathcal{Z}_a)^{N_m}} + d_F(\xi_s^1, \xi_s^2) ds. $$

As $\mathcal{X}$ is compact, according to Stone-Weierstrass theorem, the set of functions $\mathcal{C}^{1,\bullet}(\mathcal{X}, \mathbb{R})$ is dense in $\mathcal{C}(\mathcal{X}, \mathbb{R})$ for the uniform norm. So $\sup_{f \in Lip_{1,1} \cap \mathcal{C}^{1,\bullet}} |\langle \xi_t^1 - \xi_t^2, f \rangle| = d_F(\xi_t^1, \xi_t^2)$ and

$$ d_F(\xi_t^1, \xi_t^2) \leq C \int_0^t ||c_s^1 - c_s^2||_{L^2(\mathcal{Z}_a)^{N_m}} + d_F(\xi_s^1, \xi_s^2) ds. \quad (C5) $$
Using (C3), we obtain that there exists a constant C such that

$$ ||c_t^1 - c_t^2||_{L^2(\mathcal{X}_a)^{N_m}} + d_F(\xi_t^1, \xi_t^2) \le C \int_0^t ||c_s^1 - c_s^2||_{L^2(\mathcal{X}_a)^{N_m}} + d_F(\xi_s^1, \xi_s^2) ds. \quad (C6) $$

From Gronwall lemma, $||c_t^1 - c_t^2||_{L^2(\mathcal{X}_a)^{N_m}} + d_F(\xi_t^1, \xi_t^2) = 0$ which proves the uniqueness of a positive solution.

## Appendix D. Metrics on $\mathcal{M}_{F^+}(\mathcal{X})$, properties of $\mathcal{D}([0,T], (\mathcal{M}_{F^+}(\mathcal{X}), d_{PR}))$

We remind the reader that $(\mathcal{X}, d_{\mathcal{X}})$, with $d_{\mathcal{X}}$ defined at (1), is a complete separable metric space (Polish space). For a closed set $F$ of $\mathcal{X}$ and $\varepsilon \ge 0$, let $F_{\varepsilon} = \{x \in \mathcal{X} : d_{\mathcal{X}}(x, F) < \varepsilon\}$. The Prokhorov distance $d_{PR}$ on $\mathcal{M}_{F^+}(\mathcal{X})$ is defined by (Appendix A2.5 of [25])

$$ d_{PR}(\mu, \nu) = \inf\{\varepsilon \ge 0: \begin{array}{l} \text{for all closed } F \subset \mathcal{X}, \mu(F) \le \nu(F_{\varepsilon}) + \varepsilon \\ \text{and } \nu(F) \le \mu(F_{\varepsilon}) + \varepsilon. \end{array} \quad (D1) $$

A useful property of the Prokhorov distance is that it metrizes the weak convergence topology on $\mathcal{M}_{F^+}(\mathcal{X})$ (Corollary A2.5.II in [25]). Moreover, if $(X, d)$ is Polish, the space $(\mathcal{M}_{F^+}(\mathcal{X}), d_{PR})$ is Polish too (Proposition A.2.5.III of [25]). Once a metric is defined on $\mathcal{M}_{F^+}(\mathcal{X})$, we define the Skorokhod metric (see section 3.5 of [9]) on the space $\mathcal{D}([0, T], (\mathcal{M}_{F^+}(\mathcal{X}), d_{PR}))$ of càdlàg functions of $[0, T]$ to $\mathcal{M}_{F^+}(\mathcal{X})$. The space $\mathcal{D}([0, T], (\mathcal{M}_{F^+}(\mathcal{X}), d_{PR}))$ equipped with the Skorokhod metric inherits the Polish space property from $(\mathcal{M}_{F^+}(\mathcal{X}), d_{PR})$. Polish spaces are well suited to study weak convergence of random variable.

On $\mathcal{M}_F(\mathcal{X})$, we also define the total variation norm

$$ |\mu|_{tv} = \mu^{+}(\mathcal{X}) + \mu^{-}(\mathcal{X}) = \sup_{\substack{f \text{ measurable,} \\ ||f||_{\infty} \le 1}} \langle \mu, f \rangle \quad (D2) $$

where the couple $\mu^+$, $\mu^-$ is the Hahn-Jordan decomposition of $\mu$.

**Proposition D.1.** Let $v \in \mathcal{M}_F(\mathcal{X})$,

$$ |v|_{tv} = \sup_{\substack{f \text{ continuous,} \\ ||f||_{\infty} \le 1}} \langle \mu, f \rangle. \quad (D3) $$

*Proof.* $v$ is finite $L^{\infty}(\mathcal{X}, |v|) \subset L^1(\mathcal{X}, |v|)$, and continuous functions are dense in $L^1(\mathcal{X}, |v|)$. $\square$

We now remind (and prove some) useful properties of these metrics.

**Proposition D.2** (Theorem 1.2 p.96 in [9]). Let $\mu$ and $v$ be two probability measures on $\mathcal{X}$. Let $\Pi(\mu, v)$ be the set of probability measures on $\mathcal{X} \times \mathcal{X}$ with marginal distributions $\mu$ and $v$. Then

$$ d_{PR}(\mu, v) = \inf_{\gamma \in \Pi(\mu, v)} \inf\{\varepsilon > 0 : \gamma(\{x, y \in \mathcal{X}^2, d(x, y) \ge \varepsilon\}) \le \varepsilon\}. \quad (D4) $$
**Proposition D.3.** Let $(v^n)_{n \in \mathbb{N}}$ be a sequence which converges to $v$ for the Skorokhod metric in $\mathcal{D}([0, T], \mathcal{M}_{F^+}(\mathcal{X}))$. We suppose that $v \in \mathcal{C}([0, T], \mathcal{M}_{F^+}(\mathcal{X}))$. Then, $(v^n)_n$ converges to $v$ uniformly in time i.e. $\sup_{0 \le t \le T} d_{PR}(v_t^n, v_t) \to 0$.

**Proof.** See Appendix: Skorohod topology of [13]. □

**Proposition D.4.** Let $(v^n)_{n \in \mathbb{N}}$ be a sequence which converges to $v$ in $\mathcal{D}([0, T], \mathcal{M}_{F^+}(\mathcal{X}))$. Then, $(v_t^n)_n$ converges to $v_t$ for almost every $t \in [0, T]$.

**Proof.** $v$ has at most a countable set of time discontinuities, written $A$ (lemma 5.1 in [9]). Let $\lambda_n$ be change of time functions as in proposition 5.3 of [9]. Then, $\forall t_0 \in [0, T] \setminus A$,

$$d_{PR}(v_{t_0}^n, v_{t_0}) \le d_{PR}(v_{t_0}^n, v_{\lambda_n(t_0)}) + d_{PR}(v_{\lambda_n(t_0)}, v_{t_0}) \to 0.$$

Indeed, the first term of the right hand side goes to 0 by definition of the convergence in Skorohod distance, and the second as well as $v_t$ is continuous at $t_0$. □

**Proposition D.5.** The total variation distance is greater than the Prokhorov distance on $\mathcal{M}_{F^+}(\mathcal{X})$.

**Proof.** For all $\mu \in \mathcal{M}_F(\mathcal{X})$, $\sup_{B \in \mathcal{B}(\mathcal{X})} |\mu(B)| \le |\mu|_{\text{tv}}$. Let $F$ be a closed set of $\mathcal{X}$, $\mu$ and $v \in \mathcal{M}_{F^+}(\mathcal{X})$,

$$\begin{aligned}
\mu(F) &= \mu(F) - v(F) + v(F) \le \sup_{B \in \mathcal{B}(\mathcal{X})} \{|\mu(B) - v(B)|\} + v(F) \\
&\le |\mu - v|_{\text{tv}} + v(F|_{|\mu - v|_{\text{tv}}})
\end{aligned}$$

where we used the fact that $v$ is a positive measure. We can switch $\mu$ and $v$ in those expressions and we obtain that $d_{PR}(\mu, v) \le |\mu - v|_{\text{tv}}$. □

**Proposition D.6.** Let $\alpha > 0$ and $v$, $\mu \in \mathcal{M}_{F^+}(\mathcal{X})$.

$$d_{PR}(\alpha\mu, \alpha v) \le \max(\alpha, 1) d_{PR}(\mu, v).$$

**Proof.** Let $F$ be a closed set of $\mathcal{X}$. By definition,

$$\begin{aligned}
\alpha v(F) &\le \alpha \mu(F_{d_{PR}(\mu, v)}) + \alpha d_{PR}(\mu, v) \\
\Rightarrow \quad &\alpha v(F) \le \alpha \mu(F_{\max(\alpha, 1)d_{PR}(\mu, v)}) + \max(\alpha, 1)d_{PR}(\mu, v).
\end{aligned}$$

**Proposition D.7.** Let $v \in \mathcal{M}_{F^+}(\mathcal{X})$ and $(v^n)_{n \in \mathbb{N}}$ a sequence of $\mathcal{M}_{F^+}(\mathcal{X})$ such that $d_{PR}(v, v^n) \to 0$. Let $G: \mathbb{R}^d \to \mathbb{R}$ Lipschitz bounded. Then,

$$\sup_{z \in \mathcal{X}} |(v - v^n) * G(z)| \to 0.$$

**Proof.** We remind the reader that if $d_{PR}(v^n, v) \to 0$ then $|v^n|_{\text{tv}} \to |v|_{\text{tv}}$, with $|v|_{\text{tv}} = v(\mathcal{X})$ when $v \in \mathcal{M}_{F^+}(\mathcal{X})$. Indeed, if $d_{PR}(v, v^n) \le \varepsilon$, from the definition (D1), $|v(\mathcal{X}) - v^n(\mathcal{X})| \le \varepsilon$.
We start by treating the case when $\nu$ is the null measure. Then, $\forall z \in \mathcal{Z}, |\nu^n * G(z)| \le |G|_\infty |\nu^n|_{tv} \xrightarrow{n \to \infty} 0$. Now, we assume that $|\nu|_{tv} > 0$. As $|\nu^n|_{tv} \to |\nu|_{tv}$, we also suppose that $|\nu^n|_{tv} > 0 \forall n \in \mathbb{N}$. We want to use proposition D.2 p.23 which is valid for probability measures, whereas here we may have $|\nu^n|_{tv} \neq |\nu|_{tv}$. But

$$d_{PR}\left(\frac{\nu^n}{|\nu^n|_{tv}}, \frac{\nu}{|\nu|_{tv}}\right) \le d_{PR}\left(\frac{\nu^n}{|\nu^n|_{tv}}, \frac{\nu^n}{|\nu|_{tv}}\right) + d_{PR}\left(\frac{\nu^n}{|\nu|_{tv}}, \frac{\nu}{|\nu|_{tv}}\right) \\ \xrightarrow{n \to 0} 0.$$

Indeed, the first term of the right hand side goes to 0 as

$$d_{PR}\left(\frac{\nu^n}{|\nu^n|_{tv}}, \frac{\nu^n}{|\nu|_{tv}}\right) \le \left|1 - \frac{|\nu^n|_{tv}}{|\nu|_{tv}}\right|$$

by proposition D.5. The second term of the right hand side goes to 0 by proposition D.6. Let $\varepsilon > 0$, from what we just showed, for $n$ large enough,

$$d_{PR}\left(\frac{\nu}{|\nu|_{tv}}, \frac{\nu^n}{|\nu^n|_{tv}}\right) \le \varepsilon \quad \text{and} \quad \left|1 - \frac{|\nu^n|_{tv}}{|\nu|_{tv}}\right| \le \varepsilon.$$

Let $G^+(z)$ and $G^-(z)$ be two positive functions such that $G(z) = G^+(z) - G^-(z)$. Then,

$$|(\nu - \nu^n) * G^+(z)| \le |\nu|_{tv} \left( \left| \left( \frac{\nu}{|\nu|_{tv}} - \frac{\nu^n}{|\nu^n|_{tv}} \right) * G^+(z) \right| + \left| 1 - \frac{|\nu^n|_{tv}}{|\nu|_{tv}} \right| \times \left| \frac{\nu^n}{|\nu^n|_{tv}} * G^+(z) \right| \right) \\ \le |\nu|_{tv} \times \left| \left( \frac{\nu}{|\nu|_{tv}} - \frac{\nu^n}{|\nu^n|_{tv}} \right) * G^+(z) \right| + \varepsilon |\nu|_{tv} |G^+|_\infty$$

Let $\hat{v}^n := \frac{\nu^n}{|\nu^n|_{tv}}$ and $\hat{v} := \frac{\nu}{|\nu|_{tv}}$, which are probability measures on $\mathcal{X}$. From proposition D.2 p.23, $\forall \delta > \varepsilon$, there exists a probability measure on $\mathcal{X} \times \mathcal{X}$ $\gamma \in \Pi(\hat{v}^n, \hat{v})$ and a set $E \supset \{(x, y) \in \mathcal{X}^2 : d(x, y) \ge \delta\}$ such that $\gamma(E) \le \delta$. Then

$$(\hat{v}^n - \hat{v}) * G^+(z) = \sum_{l=1}^{N_c} \int_{\mathcal{Z}} \sum_{k \in [N_c]} \int_{\mathcal{Z}} G^+(z-y) - G^+(z-u) \gamma(dy, dl, du, dk) \\ \le 2\delta |G^+|_\infty + \int_{\mathcal{X} \setminus E} G^+(z-y) - G^+(z-u) \gamma(dy, dl, du, dk).$$

We can moreover assume that $\varepsilon < \delta < \text{diam}(\mathcal{X})$ and use the fact that $G^+$ is Lipschitz in space to bound the integral on the domain $\mathcal{X} \setminus E$ with $\delta |G^+|_{\text{Lip}}$. And this is true for any $\delta > \varepsilon$ so

$$(\hat{v}^n - \hat{v}) * G^+(z) \le \varepsilon |\nu|_{tv} (2 |G^+|_\infty + |G^+|_{\text{Lip}}). \quad (D5)$$

The same kind of upper bound can be found for $|(\nu^n - \nu) * G^-(z)|$. And so, $\forall z \in \mathcal{Z}$,

$$|(\nu^n - \nu) * G(z)| \le 2\varepsilon |\nu|_{tv} (3|G|_\infty + |G|_{\text{Lip}}).$$

This concludes the proof as $\varepsilon \to 0$. $\square$
Appendix E. Reaction-diffusion equations

In this section we recall some classical results on linear reaction-diffusion equations with non-homogeneous boundary conditions. The main result is proposition E.1 which states existence and uniqueness of solution for the reaction-diffusion problem (E4). We use this result to obtain the existence and uniqueness of solutions to the problem governing the evolution of chemical concentrations in our model (proposition 3.1 in the main text).

Let $U \subset \mathbb{R}^d$ be a bounded open set of regularity $\mathcal{C}^\infty$. We denote $\partial U$ the boundary of $U$ and consider the homogeneous Robin boundary condition

$$ \partial_n u|_{\partial U} + bu|_{\partial U} = 0 \qquad (E1) $$

where $b$ is a function in $C^\infty(\bar{U})$ and $\partial_n$ is the normal derivative on $\partial U$. It is a standard result (see [23, 26]) that the unbounded operator $A = \sigma\Delta$ whose domain $\mathcal{D}(A)$ is the closure in $H^2(U)$ of the set of function in $C^2(\bar{U})$ which satisfy (E1) generates a strongly continuous semi-group on $L^2(U)$, denoted $S$. Moreover, this semi-group is positive and analytical.

Let us consider the non-homogeneous elliptic problem defined on $U$:

$$ -\lambda c + Ac = 0 \qquad (E2) $$

$$ \partial_n c|_{\partial U} + bc|_{\partial U} = f \qquad (E3) $$

where $\lambda \ge 0$ is a real number greater than the largest spectral value of $A$. For all $\alpha \ge 0$, the problem (E2) has a unique solution $c$ in $H^{2\alpha}(U)$ when $f \in H^{2\alpha-\frac{3}{2}}(\partial U)$, denoted $c = D(f)$. Operator $D$ is bounded from $H^{2\alpha-\frac{3}{2}}(\partial U)$ to $H^{2\alpha}(U)$. These standard results can be found in [27] or [23].

Using notations taken from [28], we define the following functional spaces, for all $\alpha \ge 0$ and for $L^2(U) = H^0(U)$,

$$ W_{L^2(U)}^{2\alpha, [\alpha]} = \left\{ f \in L^2([0, T], H^{2\alpha}(U)); \forall j \in \{1, \dots, [\alpha]\} \quad \frac{\partial^j f}{\partial t^j} \in L^2([0, T], H^{2\alpha-2j}(U)) \right\} $$

where $[\cdot]$ is the ceiling function. $W_{L^2(U)}^{2\alpha, [\alpha]}$ is a Hilbert space endowed with the norm

$$ \left( \|f\|_{L^2([0,T], H^{2\alpha}(U))}^2 + \sum_{j=1}^{[\alpha]} \| \frac{\partial^j f}{\partial t^j} \|_{L^2([0,T], H^{2\alpha-2j}(U))}^2 \right)^{\frac{1}{2}}. $$

Note that with these notations, $W_{L^2(U)}^{0,0} = L^2([0,T] \times U)$. We define in the same way the space $W_{L^2(\partial U)}^{2\alpha, [\alpha]}$. From [27] theorem 3.1 p. 19, we know that:

* if $f \in W_{L^2(\partial U)}^{\frac{3}{2},1}$, then $f \in L^2([0,T], H^{\frac{3}{2}}(\partial U))$ and $\frac{\partial f}{\partial t} \in L^2([0,T], H^{-\frac{1}{2}}(\partial U))$, so $f \in \mathcal{C}([0,T], H^{\frac{1}{2}}(\partial U))$ and the injection is continuous.

* Similarly, if $c \in W_{L^2(U)}^{2,1}$ then $c \in \mathcal{C}([0,T], H^1(U))$ and the injection is continuous.

With these definitions and properties, we have the following result
**Proposition E.1.** If $b \in \mathcal{C}^{\infty}(\bar{U})$, $c_0 \in H^1(U)$, $f \in W_{L^2(\partial U)}^{\frac{1}{2},1}$, and $u \in W_{L^2(U)}^{0,0}$, then

$$c(.,t) = S(t)c_0 + (\lambda I - \sigma \Delta)(\int_0^t S(t-s)D(f)(.,s)ds) + \int_0^t S(t-s)u(.,s)ds$$

is the unique weak solution in $W_{L^2(U)}^{2,1}$ of the following evolution problem

$$\begin{aligned}
\partial_t c - \sigma \Delta c &= u, \\
\partial_n c|_{[0,T] \times \partial U} + bc|_{[0,T] \times \partial U} &= f \\
c(.,0) &= c_0
\end{aligned} \qquad (E4)$$

This solution belongs to $W_{L^2(U)}^{2,1} \cap \mathcal{C}([0,T], H^1(U))$. Moreover, the operator $(c_0, f, u) \mapsto c$ is continuous from $(H^1(U) \times W_{L^2(\partial U)}^{\frac{1}{2},1} \times W_{L^2(U)}^{0,0})$ to $W_{L^2(U)}^{2,1}$.

**Proof.** It is a consequence of theorem 6.1 from [28] for problems with mixed (Robin) boundary conditions in the case $\alpha = 1/4$ and $k = 0$. $\square$

It should be noted the same type of results hold for Dirichlet boundary conditions with $f \in W_{L^2(\partial U)}^{\frac{3}{2},1}$ if the compatibility condition $c_{0|\partial U} = f(0)$ is satisfied (see[28]). However, it is difficult to state general results when $U$ is non smooth or in the case of mixed Dirichlet-Neumann boundary conditions, where the analysis of each specific case is required.

**Appendix F. Technical results for the proof of theorem 3.4**

**Lemma F.1.** Let $(c_s, \xi_s)_{s \ge 0}$ be a solution of (17-18) on $[0,T]$, let $t \in [0,T]$ and $f(z) \in \mathcal{C}^1(\mathcal{L}, \mathbb{R})$. The transport equation on $[0,T]$

$$\left\{
\begin{array}{l}
\partial_s u(s,z) + \mathcal{F}(z, \xi_s * F(z)) \nabla_z u(s,z) = 0, \\
u(t,z) = f(z),
\end{array}
\right.
\qquad (F1)$$

admits a unique solution constant along characteristics $Z(s;t,z)$:

$$u(s,z) = f(Z(t;s,z)), \qquad (F2)$$

where characteristics $Z(s;t,z)$ are solutions of

$$\left\{
\begin{array}{l}
d_s Z(s;t,z) = \mathcal{F}(Z(s;t,z), \xi_s * F(Z(s;t,z))), \\
Z(t;z) = z \quad \text{for } z \in \mathcal{L}.
\end{array}
\right.
\qquad (F3)$$

Moreover, the solution $u$ belongs to $\mathcal{C}^1(|0,T| \times \mathcal{L}, \mathbb{R})$.

**Proof.** This is an application of theorem 6.1 in Ref. [29], which relies on the continuous derivability of the flow associated to the ODE (F3). To prove the latter, one must show that the vector field $(t,z) \mapsto \mathcal{F}(z, \xi_t * F(z))$ is continuously derivable on $[0,T] \times \mathcal{L}$, which is a
standard consequence of the $\mathscr{C}^1$ regularity assumption (H8). The detailed proof can be found in A.3.2 p.120 of ref. [30]. $\square$

Let $(c_t^1, \xi_t^1)$ and $(c_t^2, \xi_t^2)$ be two solutions of (17-18) on $[0, T]$ with $\xi_t^i$ positive, and let $Z^1(s; t, z)$ and $Z^2(s; t, z)$ be the associated characteristics defined in (F3). Then, the following proposition is true.

**Proposition F.2.** Let $z \in \mathscr{L}$ and $0 \le s \le t$,

$$|Z^1(t; s, z) - Z^2(t; s, z)| \le C(t) \int_s^t d_F(\xi_{u1}^1, \xi_{u2}^2) du \quad (F4)$$

with $C(t)$ a constant which depends on $t$.

**Proof.** The proof can be found in property 2.9, section A.4 p.125 of ref. [30]. It relies on Gronwall lemma and regularity hypotheses made on the different functions involved. $\square$