domain: math sdv

text:
Léo Darrigade
email: leo.darrigade@inria.fr

Simon Labarthe
email: simon.labarthe@inrae.fr

Béatrice Laroche
email: beatrice.laroche@inrae.fr

Deterministic limit of a PDMP model of epithelial tissue interacting with diffusing chemicals and application to the intestinal crypt

Keywords: biological tissue model, individual-based model (IBM), large population asymptotics AMS Subject Classification: 60J85, 92D25, 92C15

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.

Introduction

Biological tissues are meso-scale biological structures between micro-scale individual cells and macroscale organs. Tissues are characterized by interacting cells with similar physiology embedded into an extracellular matrix. Thus, tissues have their own regulation mechanisms including cell metabolism, cell-cell communications, cell movement and mechanical interactions [START_REF] Meier-Schellersheim | Multiscale modeling for biologists[END_REF]. In the framework of cellpopulation studies at the tissue scale, mathematical models of biological tissues are particularly useful for integrating data at different scales or for testing hypothesis with high-throughput in silico experiments [START_REF] Alexander | Seven challenges in the multiscale modeling of multicellular tissues[END_REF]. A large set of tissue models have then been developed in different contexts as diverse as surgical robotics [START_REF] Famaey | Soft tissue modelling for applications in virtual surgery and surgical robotics[END_REF], embryogenesis [START_REF] Ben Amar | Physics of growing biological tissues: the complex cross-talk between cell activity, growth and resistance[END_REF][START_REF] Rosenbauer | Modeling of Wnt-mediated tissue patterning in vertebrate embryogenesis[END_REF] or mathematical oncology [START_REF] Russell C Rockne | The 2019 mathematical oncology roadmap[END_REF].

Usually, tissue models are either individual-based models (IBM) or continuous population models. IBMs treat cells as distinct entities, aggregating a large number of individual cell models. These models are interconnected through cell-to-cell interactions, facilitating the emergence of patterns at the tissue scale [START_REF] James | Comparing individual-based approaches to modelling the selforganization of multicellular tissues[END_REF]. IBMs offer several key advantages. Primarily, they align well with biological observations frequently made at the cellular level, such as through microscopy or single-cell experiments, making them suitable for integrating data and knowledge. Additionally, they are particularly well adapted to objectoriented programming, allowing to leverage high-performance and parallel computing. Furthermore, IBMs efficiently incorporate stochastic events at the individual level, allowing for the exploration of inter-individual variability within populations. However, the computational and memory requirements for IBM scale at least linearly with the number of cell models, posing challenges for large cell numbers that approach real tissue scales. In such cases, the synthetic data generated by IBMs can be as complex to analyze as actual experimental data [START_REF] James | Comparing individual-based approaches to modelling the selforganization of multicellular tissues[END_REF].

On the other hand, 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). These models, however, lack important properties of the IBMs. They are most often fully deterministic (although stochastic events can technically be incorporated) whereas stochasticity is an important aspect of cellular behaviour [START_REF] Losick | Stochasticity and Cell Fate[END_REF]. Additionally, they do not facilitate the modelling of individual trajectories hindering the comparison with individual-based data. The main advantage of continuous population models is their ability to analyze population-wide traits independent of the actual cell count, allowing scalability for large populations. Moreover, their relative simplicity compared to large scale IBMs make them more amenable to mathematical analysis.

Tissue models can integrate both formalism. For instance, IBMs often include smaller-scale elements such as nutrients through a continuous model coupled to the IBM [START_REF] Aaron | Incorporating chemical signalling factors into cell-based models of growing epithelial tissues[END_REF][START_REF] James | Comparing individual-based approaches to modelling the selforganization of multicellular tissues[END_REF]. Furthermore, continuous population models can be directly derived from IBMs using large population approximations [START_REF] Stewart | Markov Processes: Characterization and Convergence[END_REF]. 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 asymptotic model derived from the 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. This qualitative behavior may then be extrapolated to the IBM for large enough populations.

The goal of this paper is twofold: we aim to propose a versatile stochastic individual-based model of epithelial tissues, that depicts a population of cells interacting with diffusing chemicals, alongside a rigorous deterministic approximation. We want the IBM general enough to be applicable to a broad range of tissues, encapsulating essential aspects of biological tissues turnover. The model is structured by cell types and spatial dimensions, enabling the representation of population distributions along these axes. Additionally, it accounts for cell fate regulation by diffusing extracellular molecules and the associated spatial variations in chemical concentrations. We therefore model : 1) spatial movement of cells driven by mechanical interactions with neighbouring cells, 2) stochastic cell division, extrusion and differentiation regulated by chemical concentrations, cell location and local cell density, 3) spatial diffusion of chemicals, 4) consumption or production of chemicals by cell metabolism.

Although these cellular functions are broadly applicable, our primary focus is on the epithelium of the intestinal colon crypt, a tissue with unique characteristics that make it an ideal test case for our model. First, the gut epithelium is one of the tissue with the highest renewal rate in the human body: the cell turn over is in the range of a few days and includes strong dynamical processes such as differentiation and cell migration [START_REF] Barker | Adult intestinal stem cells: critical drivers of epithelial homeostasis and regeneration[END_REF]. Next, a lot of different cell-cell interactions occur, including mechanical repression of division (contact inhibition), contact signaling [START_REF] Batlle | β-Catenin and TCF mediate cell positioning in the intestinal epithelium by controlling the expression of EphB/ephrinB[END_REF] and longer scale interactions with diffusing molecules [START_REF] Wang | Bioengineered systems and designer matrices that recapitulate the intestinal stem cell niche[END_REF]. Furthermore, complex interactions with the cell environment take place, in particular the uptake of gut-microbes-derived metabolites [START_REF] Kaiko | The Colonic Crypt Protects Stem Cells from Microbiota-Derived Metabolites[END_REF]. Models of the gut crypt already exist, both in the IBMs [DNO13; [START_REF] Buske | On the biomechanics of stem cell niche formation in the gut -Modelling growing organoids[END_REF][START_REF] Thalheim | Linking stem cell function and growth pattern of intestinal organoids[END_REF] and continuous population model [START_REF] Matthew D Johnston | Mathematical modeling of cell population dynamics in the colonic crypt and in colorectal cancer[END_REF][START_REF] Philip | Comparing a discrete and continuum model of the intestinal crypt[END_REF] frameworks. However, these IBMs lack a rigorous mathematical foundation, particularly in their stochastic components, limiting the scope for mathematical asymptotic analysis. We note a previous derivation of a population model from IBMs [START_REF] Philip | Comparing a discrete and continuum model of the intestinal crypt[END_REF] but the derivation is purely formal. Our paper addresses these gaps by defining an IBM within the PDMP (Piecewise Deterministic Markov Process) framework and proving convergence to the population model using functional analysis tools. This mathematical framework was first introduced in theoretical ecology and evolution [START_REF] Fournier | A microscopic probabilistic description of a locally regulated population and macroscopic approximations[END_REF][START_REF] Champagnat | Individual-based probabilistic models of adaptive evolution and various scaling approximations[END_REF][START_REF] Viet | Stochastic particle models for problems of adaptive evolution and for the approximations of statistical solutions[END_REF]. It was then adapted to other biological context with less emphasis on evolutionary dynamics such as the chemostat [START_REF] Campillo | Weak Convergence of a Mass-Structured Individual-Based Model[END_REF] or the ovarian follicle [START_REF] Michel | Large population and size scale limit of a stochastic particle model[END_REF]. Our work extends previous studies by accounting for a greater diversity of biological phenomena. It provides a strong theoretical grounding for the models used in the community for mathematical oncology, morphogenesis and gut microbiota studies [START_REF] Darrigade | A PDMP model of the epithelial cell turn-over in the intestinal crypt including microbiota-derived regulations[END_REF].

Our individual based PDMP model (PDMP-IBM) represents the population of cells as a point measure and concentrations of molecules as smooth densities. In this class of IBM continuous processes modeled by deterministic equations (like ODEs) are punctuated by stochastic events: the random process determines the time occurrence and the magnitude of stochastic jumps in the state space, with individual trajectories driven by the continuous deterministic model between jumps [START_REF] Viet | Large population limit and time behaviour of a stochastic particle model describing an age-structured population[END_REF][START_REF] Campillo | Weak Convergence of a Mass-Structured Individual-Based Model[END_REF][START_REF] Benaım | Qualitative properties of certain piecewise deterministic Markov processes[END_REF]. Our model's deterministic part includes cell-to-cell interactions and the reaction-diffusion of chemicals, while the stochastic part models cell division, extrusion and differentiation. Custom notations used throughout the paper are defined in section 2. The PDMP is introduced in section 3, with its deterministic part presented in section 3.1 and its stochastic part in 3.2. The well-posedness and stability of the process are shown in 4. We then prove in section 5 weak convergence of this stochastic cell-based model to a deterministic model of self-aggregation type, under the hypothesis of large population. This part is inspired by previous works [START_REF] Fournier | A microscopic probabilistic description of a locally regulated population and macroscopic approximations[END_REF][START_REF] Viet | Stochastic particle models for problems of adaptive evolution and for the approximations of statistical solutions[END_REF][START_REF] Campillo | Weak Convergence of a Mass-Structured Individual-Based Model[END_REF] with adaptations to take into account mechanical interactions between cells and the strong coupling of the cell population with reaction-diffusion PDEs. Finally, section 6 assesses the convergence numerically, using a model that simulates epithelial turn-over in an intestinal crypt.

Preliminary notations and definitions

We consider a finite set T of cell types, and a spatial domain Z which is a compact subset of R d , with d ≤ 3. Cells are characterised by their position and their type x = (z, l) ∈ Z ×T . We model a population of cells ν as a point measure on X = Z × T :

ν = n i=1 δ x i (dx) (1)
where δ x (dx) is the Dirac measure on x and n = ⟨ν, 1⟩ is the total number of cells in the population. We also model the spatial and temporal evolution of c, the vector of concentrations of N m ∈ N * chemicals on Z.

Topology and metric on Z and X

We endow Z ⊂ R d with the topology associated to the euclidean norm and denote B(Z) the corresponding Borel σ-algebra, and endow T with the topology generated by P(T ), the set of its parts. The product space X = Z × T is equipped with the product topology for which an associated metric d X is defined as

d X : X × X → R + ((z 1 , l 1 ), (z 2 , l 2 )) → |z 1 -z 2 | if l 1 = l 2 , β otherwise, (2)
with β > diam(Z) and | • | the euclidean norm on Z. Then the product σ-algebra B(X ) = B(Z) ⊗ P(T ) is the Borel σ-algebra on (X , d X ). Moreover, we sometimes use Ẑ := {z -y : (z, y) ∈ Z 2 }.

Functions defined on X

Let f be a function defined on X . For all l in T , we denote f l : z → f (z, l) the restriction of f to Z × {l}.

Then f is:

• measurable (resp. continuous) iff, for all l ∈ T , f l is measurable (resp. continuous) with respect to the Borel σ-algebra (resp. the euclidian norm) on Z.

• f ∈ C k,• (X , R) iff, for all l ∈ T , f l is continuous of class k in Z (f l ∈ C k ), meaning that it is k times differentiable on Z, with k-th derivative continuous on Z. When k ≥ 1, we denote ∇ z f (x) ∈ R d
the gradient operator with respect to the space variable z.

• f ∈ C k1,k2,• (R + × X , R) iff, for all l ∈ T , f l (t, z) is continuous of class k 1 with respect to t ∈ R +
and continuous of class k 2 with respect to z ∈ Z.

For f : X → R, we write

||f || ∞ = sup x∈X |f (x)| and ||f || Lip = sup (x,y)∈X ,x̸ =y |f (x) -f (y)| d X (x, y) .
We have the following elementary proposition.

Proposition 1. f : X → R is Lipschitz on X iff, ∀l ∈ T , f l is Lipschitz on Z. Moreover, ||f || Lip ≤ max max l∈T ||f l || Lip(Z) , 2||f || ∞ β with ||f l || Lip(Z) is the Lipschitz constant of f l on Z and β is defined in (2).
For f Lipschitz on X , we define

||f || ∨ = max(||f || ∞ , ||f || Lip ) (3)
and introduce the space of functions bounded by 1 with Lipschitz constant lower than 1

Lip 1,1 (X , R) = {f : X → R : ||f || ∨ ≤ 1}.

Measures on X

We denote M F (X ) (M F + (X )) the space of finite (resp. finite positive) measures on (X , B(X )), and M P (X ) the subset of finite positive point measures:

M P (X ) = n i=1 δ x i ; n ∈ N, x i = (z i , l i ) ∈ Z × T , with 0 i=1 δ x i (x) the null measure by convention.
Let ν be a measure of M F (X ), then, for all l ∈ T , we define ν l ∈ M F Z as

ν l (Z) = ν(Z × {l})
with Z any Borel set of Z. Then:

ν(dz, dl) = k∈T δ k (dl)ν k (dz) (4)
as for any set X of B(X ), X = ∪ l∈T Z l × {l} with Z l a Borel set of Z.

For a measurable function f : X → R and a measure ν ∈ M F (X ), we define

⟨ν, f ⟩ := X f (x)ν(dx) = l∈T Z f (z, l)ν(dz, l). (5) If ν ∈ M P (X ), that is to say ν = n i=1 δ x i for some n ∈ N, ⟨ν, f ⟩ = n i=1 f (x i ) = n i=1 f li (z i ).
In particular, ⟨ν, 1⟩ = n.

Convolution

For a measurable function D : R d × T → R and ν ∈ M F (X ), we set for

x = (z, l) ∈ X ν * D(x) := l∈T Z D(z -y, l)ν(dy, l). (6) If ν ∈ M P (X ), ν * D(x) = ⟨ν,1⟩ i=1 D(z -z i , l i ).

Definition of the process

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

Deterministic movement of cells and reaction-diffusion of chemicals

Movement of cells

For each cell i with mass m i << 1, the position of its center z i (t) evolves continuously in the spatial domain Z, in a non-inertial regime:

¨¨¨m i d 2 z i dt 2 (t) inertia + dz i dt (t) friction = F( z i (t), ν t * F (z i (t)) )
other forces

(7) with F : R d → R d and F : Z × R d → R d .
The first variable z i (t) in F is introduced to take into account local modulations (if any) in the total force (apart from friction). 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

n j=1,j̸ =i F (z i -z j ) = n j=1 F (z i -z j ) = ν * F (z i ). ( 8
)
We assume that Hypothesis 2.

1. F ∈ C 1 (Z × R d , R d ), F ∈ C 1 (R d , R d ), 2. F (0) = 0, 3. F(z, •) = 0 ∀z ∈ ∂Z.
We have the following standard result.

Proposition 3 (Cauchy-Lipschitz). Thanks to hypothesis (3) and regularity assumptions on F and F, the system of coupled equations ( 7) has a unique solution on R + , taking values in Z n , for any initial condition in Z n .

Remark 4. Proposition 3 holds under the weaker hypothesis that F and F are locally Lipschitz, the C 1 regularity will be used later in section 5.

Remark 5. Assumption (2) deserves an explanation. It means that two cells can occupy the same location without exerting a force on each other. Such situation occurs when the space dimension is reduced due to symmetry assumptions, in which case cells are allowed to form vertical or spherical layers corresponding to the same altitude or radial distance, see for instance [START_REF] Darrigade | A PDMP model of the epithelial cell turn-over in the intestinal crypt including microbiota-derived regulations[END_REF]. More generally, we assume that when they get very close, cells exert a repulsive force on each other. However, when the distance between cells is too small, this correspond to a high compression and cells are damaged, so that they cannot exert a force on each other any more. This usually comes with a high probability of extrusion (expulsion of damaged cells), as reminded in the description of cell fate events below.

Remark 6. We could take into account an effect of cell types in mechanical interaction by using for the force F (z i -z j , l i , l j ) with F : R d × T × T → R d , and l i , l j are respectively the cell types of cell i and cell j. By extending hypotheses done on F to F (•, l i , l j ), all of the following results are still valid.

Reaction-diffusion of chemicals

Let a be the typical diameter of a cell and B R d (z, r) be the open ball of R d centered on z of radius r. We define

Z a := z∈Z B R d z, 1 2 a
the spatial domain we use for modelling concentrations of chemicals. It is an enlarged version of Z which contains not only the position of the center of cells, as Z, but also the whole cells area (or volume), which is appropriate to model surfacic absorption of chemicals. We assume the following:

Hypothesis 7. Z a is an open set of R d of regularity C ∞ , with boundary Γ.
Let c(z) be the vector of concentrations of the N m chemicals located at z. Each cell type l ∈ T can degrade or produce chemical m, 1 ≤ m ≤ N m , at speed γ ml (c(z)). A cell i located at z i can catalyze a reaction only if it is physically close to z. Therefore, we weight the reaction speed γ ml (c(z)) with a term ψ a (z -z i ) where Hypothesis 8. ψ a is a positive function with radial symmetry, Lipschitz on R d with support contained in B R d (0, 1 2 a) and its integral is 1. On the boundary of Z a , concentrations are influenced by the external media through quantities f m which depend of time and space. Then, concentrations c on Z a evolve according to the following system of coupled reaction-diffusion equations

∀1 ≤ m ≤ N m , ∂ t c m -σ∆c m = l∈T γ ml (c) ν l t * ψ a (z), (9)
∂ n c m,|[0,T ]×Γ + b m c m,|[0,T ]×Γ = f m , c m (•, 0) = c 0m ,
where ∂ n is the normal derivative on Γ.

For an open set U ⊂ R d and α ≥ 0, let

W 2α,⌈α⌉ L 2 (U ) := f ∈ L 2 ([0, T ], H 2α (U )); ∀j ∈ {1, • • • , ⌈α⌉} ∂ j f ∂t j ∈ L 2 ([0, T ], H 2α-2j (U ))
with ⌈•⌉ the ceiling function.

We assume that Hypothesis 9. For all

1 ≤ m ≤ N m , 1. b m is in C ∞ (Z a ),
2. ∀l ∈ T , γ ml is bounded and uniformly Lipschitz from R Nm

+ to R, 3. the initial condition c 0m is in H 1 (Z a ), 4. f m is in W 1 2 ,1 L 2 (Γ)
. With the notations introduced in section B, we denote S m the analytic semi-group S on L 2 (Z a ) generated by σ∆ with boundary conditions of chemical m, D m the related lifting operator (D in the appendix) and λ m the associated scalar (λ in the appendix). Moreover, for ν ∈ M F (X ), let ||ν|| T V be the total variation norm defined in equation (63) of the appendix. Then the following proposition is verified.

Proposition 10. Under hypothesis 9, for all ν measurable on [0, T ] with values in M F + (X ) endowed with the weak convergence topology and the associated σ-algebra, and such that sup t∈[0,T ] ||ν t || T V < +∞, the non-linear evolution problem (9) has a unique solution c in W 2,1

L 2 (Za) Nm ∩ C([0, T ], H 1 (Z a ) Nm ), verifying for 1 ≤ m ≤ N m : c m (., t) = S m (t)c 0m + (λ m I -∆)( t 0 S m (t -s)D m (f m )(., s)ds) (10) + t 0 S m (t -s) l∈T γ ml (c(s)) ν l s * ψ a (•)ds Proof. Let us define the function G on [0, T ] × L 2 (Z a ) Nm as G(t, c) := l∈T γ 1l (c) ν l t * ψ a (•), • • • , l∈T γ Nml (c) ν l t * ψ a (•) ∈ L ∞ ([0, T ], L 2 (Z a ) Nm ).
G is Lipschitz with respect to c on L 2 (Z a ) Nm , uniformly for t ∈ [0, T ]. Indeed, as γ ml are Lipschitz under hypothesis 2,

||γ|| Lip = sup m,l ||γ ml || Lip < ∞. ( 11
) Let c 1 , c 2 in L 2 (Z a ) Nm , then sup t∈[0,T ] ∥G(t, c 1 ) -G(t, c 2 )∥ L 2 (Za) Nm ≤ Nm m=1 l∈T sup t∈[0,T ] ||ν l t || T V ∥γ ml (c 1 ) -γ ml (c 2 )∥ L 2 (Za) ≤ sup t∈[0,T ] ||ν t || T V × N m |T |K∥c 1 -c 2 ∥ L 2 (Za) Nm .
With a standard fixed point argument (see [CZ12; Ama93]) on the system (10) and the semi-group of L 2 (Z a ) Nm defined by S = diag(S 1 , . . . , S Nm ), we prove the existence of a weak solution c in

C([0, T ], L 2 (Z a ) Nm ). Then u(z, t) = l∈T γ ml (c(z, t))ν l t * ψ a (z) ∈ L ∞ ([0, T ], L 2 (Z a )) ⊂ W 0,0 L 2 (Za) .
Indeed, γ ml are bounded lipschitz on R Nm + under hypothesis 2, so

γ ml (c(z, t)) ∈ C([0, T ], L 2 (Z a )) as c ∈ C([0, T ], L 2 (Z a ) Nm ) according to the previous step. Moreover, ν l t * ψ a (z) ∈ L ∞ ([0, T ], L 2 (Z a )) because ψ a is bounded and sup t∈[0,T ] ||ν t || T V < ∞.
We obtain the adequate regularity result by observing that c m is a solution of a problem of type (73) with c 0m ∈ H 1 (Z a ), f m verifying the regularity conditions stated in proposition 35 and u ∈ W 0,0

L 2 (Za) . Therefore c m ∈ W 2,1 L 2 (Za) ∩ C([0, T ], H 1 (Z a )).
The following lemma will be used later. It bounds the distance between two populations of cells ν 1 and ν 2 , by a quantity depending on the distance between the two corresponding solutions of (9). Lemma 11. Let ν 1 , ν 2 be two measurable functions from [0, T ] to M F + (X ) (endowed with the topology of weak convergence and the associated Borel σ-algebra) such that

sup t∈[0,T ] ||ν 1 t || T V ≤ B < +∞ and sup t∈[0,T ] ||ν 2 t || T V ≤ B < +∞.
Let c 1 et c 2 be the associated solutions of (9), with c 0 and f verifying the regularity hypothesis 9. Then, ∀t ∈ [0, T ],

∥c 1 (•, t) -c 2 (•, t)∥ L 2 (Za) Nm ≤ N m e Ct γ ∞ t 0 l∈T ||ν 1l s * ψ a (•) -ν 2l s * ψ a (•)|| L 2 (Za) ds, ( 12
)
where C = N m ||ψ a || ∞ ||γ|| Lip B and γ ∞ are positive constants. We can deduce that

∥c 1 (•, t) -c 2 (•, t)∥ L 2 (Za) Nm ≤ N m γ ∞ ||ψ a || ∞ |Z a | 1/2 |T |e Ct t 0 d F (ν 1 s , ν 2 s ) ds, ( 13
)
where d F is the flat metric (see A) defined on M F + (X ) by

d F (ν, ξ) = sup f ∈Lip1,1(X ,R) ⟨ν -ξ, f ⟩.
Proof. By formula (10),

∥c 1 m (•, t) -c 2 m (•, t)∥ L 2 (Za) = || t 0 S m (t -s) l∈T γ ml (c 1 (•, s))ν 1l s * ψ a (•) -γ ml (c 2 (•, s))ν 2l s * ψ a (•) ds|| L 2 (Za) = || t 0 S m (t -s) l∈T γ ml (c 1 (•, s)) ν 1l s * ψ a (•) -ν 2l s * ψ a (•) -γ ml (c 2 (•, s)) -γ ml (c 1 (•, s)) ν 2l s * ψ a (•) ds|| L 2 (Za) ≤ t 0 l∈T ||γ ml (c 1 (•, s)) ν 1l s * ψ a (•) -ν 2l s * ψ a (•) || L 2 (Za) + || γ ml (c 2 (•, s)) -γ ml (c 1 (•, s)) ν 2l s * ψ a (•)|| L 2 (Za) ds
because S m is a contraction semi-group (see B). Moreover, there exists a constant γ ∞ such that ∀m, l,

|γ ml (•)| ≤ γ ∞ < ∞. Therefore, ∥c 1 m (•, t) -c 2 m (•, t)∥ L 2 (Za) ≤γ ∞ t 0 l∈T ||ν 1l s * ψ a (•) -ν 2l s * ψ a (•)|| L 2 (Za) ds +||ψ a || ∞ ||γ|| Lip sup [0,T ] ||ν 2 s || T V t 0 ∥c 1 (•, s) -c 2 (•, s)∥ L 2 (Za) Nm ds.
We set C = N m ||ψ a || ∞ ||γ|| Lip B and deduce from the previous inequality that

∥c 1 (•, t) -c 2 (•, t)∥ L 2 (Za) Nm ≤N m γ ∞ t 0 l∈T ||ν 1l s * ψ a (•) -ν 2l s * ψ a (•)|| L 2 (Za) ds +C t 0 ∥c 1 (•, s) -c 2 (•, s)∥ L 2 (Za) Nm ds By Gronwall lemma, ∀t ∈ [0, T ], ∥c 1 (•, t) -c 2 (•, t)∥ L 2 (Za) Nm ≤ N m e Ct γ ∞ t 0 l∈T ||ν 1l s * ψ a (•) -ν 2l s * ψ a (•)∥ L 2 (Za) ds, ( 14
)
which is the first result of the lemma. We obtain the second as for all 0 ≤ s ≤ t,

||ν 1l s * ψ a (•) -ν 2l s * ψ a (•)|| L 2 (Za) ≤ ||ψ a || ∞ |Z a | 1/2 d F (ν 1 s , ν 2 s )
where |Z a | is the Lebesgue measure of Z a .

Stochastic model of cell fate

Cells can divide symmetrically, differentiate or die. These phenomena are modeled as independent random events, indexed by 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

ν t = ν t -+ µ k (x i (t -))
where

µ k : X → M F + (X ), (15)
and

• if k is an extrusion: µ k (x) = -δ x , • if k is a differenciation of a cell of type l into a type l * : µ k (z, l) = -δ (z,l) + δ (z,l * ) ,
• if k is a division: we introduce a new daughter cell at a distance λ(z) from mother cell that is to say µ k (z, l) = δ (z+λ(z),l) . We assume the following

Hypothesis 12. -λ : Z → R d is continuous and such that, ∀z ∈ Z, z + λ(z) ∈ Z, -λ is locally Lipschitz from Z to R d .
Remark 13. We conjecture that our results could be extended to a stochastic division shift λ, as done in [START_REF] Fournier | A microscopic probabilistic description of a locally regulated population and macroscopic approximations[END_REF][START_REF] Viet | Stochastic particle models for problems of adaptive evolution and for the approximations of statistical solutions[END_REF][START_REF] Campillo | Weak Convergence of a Mass-Structured Individual-Based Model[END_REF].

For a cell at x ∈ X , in a population ν ∈ M F + (X ), the individual jump rate associated to jump k ∈ E takes the form

q ∞ k × q k (x, ν * D(z), c * ψ a (z)) (16)
with

D : R d → R + , q k : X × R × R + → R + ,
ψ a the function representing the cell volume defined in hypothesis 8, and

q ∞ k < +∞ are constants. Moreover, c * ψ a (z) = Za c 1 (y)ψ a (z -y)dy, • • • , Za c Nm (y)ψ a (z -y)dy
is the vector of each coordinate of c convoluted with ψ a . Therefore, we model a cell population where jump rates depend of the cell state x, the local density of the population through the term ν * D(z), and concentrations of chemicals near the cell through the term c * ψ a (z). For instance, as mentioned earlier, a high local density could promote extrusion events as observed experimentally [START_REF] George T Eisenhoffer | Crowding induces live cell extrusion to maintain homeostatic cell numbers in epithelia[END_REF].

We assume the following Hypothesis 14.

• For all k ∈ E,

q k ≤ 1, • For all k ∈ E, for all l ∈ T , (z, y, c) → q k ((z, l), y, c) is Lipschitz on Z × R × R. We denote K k,l be a Lipschitz constant and set C(q k ) = max l∈T K k,l . • D ∈ C 1 (R d , R + )

Rescaled process

We introduce a sequence of processes (c N t , ν N t ) N ∈N * similar to the model described in sections 3.1-3.2 with parameters F N , D N , γ N

• such that

F N = F N , D N = D N , γ N • = γ • N , ( 17
)
leaving the other model parameters unchanged. Note that N = 1 corresponds to the process described in sections 3.1-3.2. The rescaled process is defined as

νN t := 1 N ν N t .
Remark 15. Obviously, the global existence and uniqueness properties of cell spatial trajectories and concentrations of chemicals as stated in propositions 3 and 10 readily extended to

(c N t , ν N t ) N ∈N * , for all N ∈ N * . We consider independent Poisson random measure on [0, 1] × N * × R + , one for each type of jump k ∈ E, denoted N k (dθ, di, dt), with intensity n k (dθ, di, ds) = q ∞ k dθdt k∈N * δ k (di).
From hypothesis 17, for all z ∈ Z,

ν N t * F N (z) = ν N t * F (z) N = νN t * F (z), ν N t * D N (z) = νN t * D(z),
and, for all (m, l)

∈ {1, • • • , N m } × T , γ N ml (c N ) ν lN t * ψ a (z) = γ ml (c N ) νlN t * ψ a (z).
Provided trajectories exist on R + (which will be proved in proposition 17 below), the process

(c N t , νN t ) t ∈ C([0, T ] , H 1 (Z a ) Nm ) × D([0, T ] , M F + (X )) is solution of, ∀f ∈ C 1,1,• (R + × X , R), ⟨ν N t , f t ⟩ =⟨ν N 0 , f 0 ⟩ + t 0 ⟨ν N s , F(•, νN s * F ) ∇ z f s + ∂ s f s ⟩ ds + k∈E t 0 N * 1 0 ⟨ν N s -+ µ k N (x N i s -), f s ⟩ -⟨ν N s -, f s ⟩ × 1 {i≤⟨ν s -,1⟩} (i)1 {θ≤q k (x N i s -,ν N s - * D(z N i s -),c N s * ψa(z N i s -))} (θ)N k (dθ, di, ds), (18)
and, for chemical m,

1 ≤ m ≤ N m , ∂ t c N m -σ∆c N m = l∈T γ ml (c N ) νlN t * ψ a (z), ( 19
)
∂ n c N |[0,T ]×Γ + b m c N |[0,T ]×Γ = f m , c N m (•, 0) = c 0m .
Remark 16. The motivation for employing test functions f t depending on both space and time comes from our proof of the convergence theorem 21, more specifically in lemma 26, second step.

Without stochastic jumps (or between jumps), the deterministic flow of the process is defined on R + for all N ∈ N * (see remark 15). For (c N t , νN t ) to have trajectories on 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 shown in the next proposition.

Proposition 17.

1. (Control of the population size) Let t ≥ 0 and p

∈ N * . If sup N ∈N * E ⟨ν N 0 , 1⟩ p < +∞, then sup N ∈N * E sup s≤t ⟨ν N s , 1⟩ p ≤ C(p, t) < +∞ ( 20
)
with C(p, t) a constant that depends of p and t.

(Stability)

For all N ∈ N * , the process (c N t , νN t ) solution of (18-19) is stable and trajectories exist on R + .

Proof. (1) : Let p, N ∈ N * . For all n ∈ N, we define the stopping-time

τ N n = inf{t ≥ 0, ⟨ν N t , 1⟩ = n}. The process (c N t∧τ N n , νN t∧τ N n
) t is well-defined as the cell population is bounded by n and, by property 15, c N exists.

Note that the determinist evolution does not change population size, and jumps modify it at most by one. Let t ≥ 0, applying Itô formula (theorem 5.

1 chapter II in [IW89]) to ⟨ν N t∧τ N n , 1⟩ p , we get ⟨ν N t∧τ N n , 1⟩ p ≤ ⟨ν N 0 , 1⟩ p + k∈E t∧τ N n 0 N * 1 0 ⟨ν N s -, 1⟩ + 1 N p -⟨ν N s -, 1⟩ p 1 {i≤⟨ν N s -,1⟩} (i) 1 {θ≤q k (•)} (θ) N k (dθ, di, ds).
The integral is increasing with time so we can switch to sup s≤t∧τ N n on the left hand-side. For positive x,

(x + 1 N ) p -x p ≤ C(p) N p ((N x) p-1 + 1) ≤ C(p) N (x p-1 + 1) with C(p) = p-1 k=1 p k . Therefore, sup s≤t∧τn ⟨ν N s , 1⟩ p ≤ ⟨ν N 0 , 1⟩ p + k∈E t∧τ N n 0 N * 1 0 C(p) N ⟨ν N s -, 1⟩ p-1 + 1 1 {i≤⟨ν N s -,1⟩} (i)
1 {θ≤q k (•)} (θ)N k (dθ, di, ds).
Taking expectation:

E sup s≤t∧τn ⟨ν N s , 1⟩ p ≤ E ⟨ν N 0 , 1⟩ p + k∈E E   t∧τ N n 0 ⟨ν N s -,1⟩ i=1 C(p) N ⟨ν N s -, 1⟩ p-1 + 1 q ∞ k q k (•) ds   ≤ E ⟨ν N 0 , 1⟩ p + k∈E E t 0 q ∞ k C(p) ⟨ν N s∧τ N n , 1⟩ p + ⟨ν N s∧τ N n , 1⟩ ds ≤ E ⟨ν N 0 , 1⟩ p + 2C(p) k∈E q ∞ k t 0 E sup u≤s∧τ N n ⟨ν N u , 1⟩ p ds.
as for p ≥ 1, n + n p ≤ 2n p . Using Gronwall lemma and our hypotheses, we can find a bound which does not depend of n nor N :

E sup s≤t∧τ N n ⟨ν N s , 1⟩ p ≤ E ⟨ν N 0 , 1⟩ p exp 2C(p) k∈E q ∞ k t = C(p, t). ( 21
)
Therefore τ n → ∞ a.s. Then, by Fatou lemma:

E sup s≤t ⟨ν N s , 1⟩ p = E lim inf n→∞ sup s≤t∧τ N n ⟨ν N s , 1⟩ p ≤ lim inf n→∞ E sup s≤t∧τ N n ⟨ν N s , 1⟩ p ≤ C(p, t).
(2) : A direct consequence of (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.

We now establish a martingale property which will be used in our convergence study. To do so, we introduce for each jump type k ∈ E the compensated random Poisson measure:

Ñk (dθ, di, ds) = N k (dθ, di, ds) -n k (dθ, di, ds). Proposition 18. We suppose that E ⟨ν N 0 , 1⟩ 2 < ∞. For all f ∈ C 1,1,• (R + × X , R) and t ≥ 0, ⟨ν N t , f t ⟩ = ⟨ν N 0 , f 0 ⟩ + t 0 ⟨ν N s , F(•, νN s * F )∇ z f s + ∂ s f s ⟩ds (22) + k∈E q ∞ k t 0 X ⟨µ k (x), f s ⟩q k (x, νN s * D(z), c N s * ψ a (z)) νN s (dx)ds + Z f,N t where Z f,N t is the square integrable martingale Z f,N t = k∈E t 0 N * 1 0 ⟨µ k (x N i s -), 1 N f s ⟩1 {i≤⟨ν N s -,1⟩} (i)1 {θ≤q k (•)} (θ) Ñk (dθ, di, ds), ( 23
)
with quadratic variation

⟨Z f,N ⟩ t = 1 N k∈E q ∞ k t 0 X ⟨µ k (x), f s ⟩ 2 q k (x, νN s * D(z), c N s * ψ a (z))ν N s (dx)ds. ( 24
)
Proof. The integrands in (23) are left-continuous and therefore predictable. From proposition 17, for all t ≥ 0,

E k∈E t 0 N * 1 0 ⟨µ k (x N i s -), 1 N f s ⟩ 1 {i≤⟨ν N s -,1⟩} (i)1 {θ≤q k (•)} (θ)n k (dθ, di, ds) = k∈E t 0 X E q ∞ k ⟨µ k (x), 1 N f s ⟩ q k (x, νN s * D(z), c N s * ψ a (z)) νN s (dx)ds < +∞
and

E k∈E t 0 N * 1 0 ⟨µ k (x N i s -), 1 N f s ⟩ 2 1 {i≤⟨ν N s -,1⟩} (i)1 {θ≤q k (•)} (θ)n k (dθ, di, ds) < +∞.
We can then apply the classical result given p.62 of [START_REF] Ikeda | Stochastic differential equations and diffusion processes[END_REF].

Remark 19. Note that results presented in this section are valid when the rate functions q k are measurable. The stronger hypothesis 14 will prove necessary in section 5

5 Deterministic limit of the rescaled process by a large population approximation.

Let T be a positive time. Let D([0, T ] , M F + (X )) be the space of càdlàg functions of [0, T ] to M F + (X ). We now prove the weak convergence of the rescaled process sequence (c N t , νN t ) N ∈N * defined in section 4 towards a deterministic process in the space

C([0, T ] , H 1 (Z a ) Nm ) × D([0, T ] , M F + (X )).
We endow D([0, T ] , M F + (X )) with the Skorokhod metric, using the Prokhorov distance d P R on M F + (X ). The topology associated to the Prokhorov metric is the topology of the weak convergence on M F + (X ). For more details on these points, see A.

We suppose that hypotheses 2, 7, 8, 9, 12, 14 are verified. Additionally, we assume that :

Hypothesis 20.

• the sequence νN 0 = 1 N ν N 0 , seen as a sequence of random variables on M F + (X ), converges weakly to ξ 0 ∈ M F + (X ) where ξ 0 can be random.

• There exists an integer p ≥ 2 such that sup N ∈N * E ⟨ν N 0 , 1⟩ p < +∞. We can now state the main result of the paper.

Theorem 21. Under all the hypotheses mentioned above, the process

(c N t , νN t ) t≥0 converges weakly on C([0, T ] , H 1 (Z a ) Nm ) × D([0, T ] , M F + (X )) to the solution of ⟨ξ t , f t ⟩ = ⟨ξ 0 , f 0 ⟩ + t 0 ⟨ξ s , F(•, ξ s * F )∇ z f s + ∂ t f s ⟩ds + k∈E t 0 l∈T Z ⟨µ k (z, l), f s ⟩q ∞ k q k ((z, l), ξ s * D(z), c s * ψ a (z))ξ s (dz, l)ds (25) for all f ∈ C 1,1,• (R + × X , R), and ∀1 ≤ m ≤ N m , ∂ t c m -σ∆c m = l∈T γ ml (c) ξ l t * ψ a (z), ( 26
)
∂ n c m,|[0,T ]×Γ + b m c m,|[0,T ]×Γ = f m , c m (•, 0) = c 0m ,
The proof of theorem 21 is established through three successive steps:

• the sequence of laws of (ν N ) N ∈N * , considered as probability measures on D([0, T ] , M F + (X )), is tight (lemma 22).

• All accumulation points of (ν N ) N ∈N * are almost surely solutions of (25), and the corresponding accumulation points in the concentration space are solutions of (26) (lemma 25).

• Uniqueness of the solution of (25-26) in lemma 26.

Lemma 22. The sequence of laws of (ν N ) N ∈N * , considered as probability measures on D([0, T ] , M F + (X )), is tight.

Proof. Based on a theorem from Roelly (Thm. 2.1 in [START_REF] Roelly-Coppoletta | A criterion of convergence of measure-valued processes: application to measure branching processes[END_REF]), the tightness of (ν N ) N ∈N * is equivalent to the tightness of laws of (⟨ν N , f ⟩) N ∈N * on D([0, T ] , R), for all f in a dense set of C(X , R) endowed with the topology of uniform convergence. Here, we choose f ∈ C 1,• (X , R).

A sufficient condition for this property is given by the Aldous-Rebolledo criteria (Cor. 2.3.3 in [START_REF] Joffe | Weak Convergence of Sequences of Semimartingales with Applications to Multitype Branching Processes[END_REF]) according to which we have to establish that 1) for t ≥ 0, the sequence (⟨ν N t , f ⟩) N ∈N * is tight on R, 2) considering the semi-martingale decomposition defined at Prop. 18

⟨ν N t , f ⟩ = ⟨ν N 0 , f ⟩ + A f,N t + Z f,N t with A f,N t
a finite variation process and Z f,N t a martingale, then for t ≥ 0, ϵ > 0, η > 0, there exists δ > 0 and N 0 ∈ N * such that for all sequence of stopping times (τ

N ) N ∈N * (τ N is a F N -stopping time) with τ N ≤ t we have sup N ≥N0 sup θ∈[0,δ] P (|A f,N τ N +θ -A f,N τ N | ≥ η) ≤ ϵ, (27)
sup

N ≥N0 sup θ∈[0,δ] P (|⟨Z f,N ⟩ τ N +θ -⟨Z f,N ⟩ τ N | ≥ η) ≤ ϵ. ( 28
)
Proof of 1): Let f ∈ C 1,• (X , R) and t ≥ 0. Under Hyp. 20, from Prop. 17.1 p.10, there exists an integer p ≥ 2 and a constant C(p, t) such that sup N ∈N * E sup s≤t ⟨ν N s , 1⟩ p ≤ C(p, t). Then, for every integer q ≤ p,

sup N ∈N * E sup s≤t ⟨ν N s , 1⟩ q ≤ C(p, t) + 1. ( 29
)
From Markov inequality, for any N ,

P (⟨ν N t , f ⟩ ≥ K) ≤ 1 K ||f || ∞ E ⟨ν N t , 1⟩ ≤ 1 K ||f || ∞ (C(p, t) + 1).
Therefore, the sequence (⟨ν

N t , f ⟩) N ∈N * is tight. Proof of 2): Let t, η, ϵ, δ ≥ 0, and θ ≤ δ, τ N < t. From proposition 18, E |A f,N τ N +θ -A f,N τ N | ≤ E | τ N +θ τ N ⟨ν N s , F(•, νN s * F )∇ z f ⟩ds| + k∈E E | τ N +θ τ N X q k (x, νN s * D(z), c N s * ψ a (z))⟨µ k (x), f ⟩ νN s (dx)ds| A consequence of hypothesis 2 is that there exists a constant C(F) such that, ∀(z, y 1 , y 2 ) ∈ Z × R d × R d , |F(z, y 1 ) -F(z, y 2 )| ≤ C(F)|y 1 -y 2 |,
and so

|F(z, y 1 )| ≤ C(F) × |y 1 | + F(z, 0) ≤ C(F) × |y 1 | + ||F(•, 0)|| ∞ < +∞ as z → F(z, 0) is continuous on Z compact. Then, |⟨ν N s , F(•, νN s * F )∇ z f ⟩| ≤ ||∇ z f || ∞ ⟨ν N s , |F(•, νN s * F )|⟩ ≤||∇ z f || ∞ ⟨ν N s , C(F)|ν N s * F (•)| + ||F(•, 0)|| ∞ ⟩ ≤||∇ z f || ∞ (C(F)||F | Ẑ || ∞ ⟨ν N s , 1⟩ 2 + ||F(•, 0)|| ∞ ⟨ν N s , 1⟩) as, ∀ z ∈ Z, |ν N s * F (z)| ≤ ||F | Ẑ || ∞ ⟨ν N s , 1⟩. Therefore, E |A f,N τ N +θ -A f,N τ N | ≤ ||∇ z f || ∞ C(F)||F | Ẑ || ∞ E τ N +θ τ N ⟨ν N s , 1⟩ 2 ds + (||∇ z f || ∞ ||F(•, 0)|| ∞ + k∈E 2q ∞ k ||f || ∞ )E τ N +θ τ N ⟨ν N s , 1⟩ds .
We have that

τ N +θ τ N ⟨ν N s , 1⟩ 2 ds ≤ δ sup s≤t ⟨ν N s , 1⟩ 2
. So, applying proposition 17.1 under Hyp. 20, we get that

E |A f,N τ N +θ -A f,N τ N | ≤ δ × ||∇ z f || ∞ C(F)||F | Ẑ || ∞ (C(p, t + δ) + 1) + δ × (||∇ z f || ∞ ||F(•, 0)|| ∞ + k∈E 2q ∞ k ||f || ∞ )(C(p, t + δ) + 1).
As δ can be made arbitrary small and the upper bound is uniform in N and θ, we get (27) by Markov inequality.

Following the same reasoning, we get that

E |⟨Z N ⟩ τ N +θ -⟨Z N ⟩ τ N | ≤ δ N × (C(p, t + δ) + 1) k∈T 2q ∞ k ||f 2 || ∞
and we deduce (28), which ends the proof of the lemma.

We now define two mappings. The first one is

Φ : D([0, T ], M F + (X )) → W 2,1 L 2 (Za) Nm ∩ C([0, T ], H 1 (Z a ) Nm ) ν → Φ(ν) = c, ( 30
)
which associates to ν the solution of the reaction-diffusion equation (9), with fixed c 0 and boundary conditions satisfying hypothesis 9. The second one is defined for a fixed t

∈ [0, T ], f ∈ C 1,1,• ([0, T ]×X , R) by ∀ν ∈ D([0, T ], M F + (X )), Ψ t (ν) = ⟨ν t , f t ⟩ -⟨ν 0 , f 0 ⟩ - t 0 ⟨ν s , F(•, ν s * F (•))∇ z f s + ∂ s f s ⟩ds - k∈E t 0 X ⟨µ k (x), f s ⟩q k (x, ν s * D(z), Φ(ν) s * ψ a (z))ν s (dx)ds. ( 31
)
Let us establish two technical lemmas.

Lemma 23. Φ is continuous from D([0, T ], M F + (X )) to C([0, T ], H 1 (Z a ) Nm ).
Proof. Let (ν N ) N be a sequence which converges to a given ξ in D([0, T ], M F + (X )). We denote c N = Φ(ν N ) and c = Φ(ξ) ∈ C([0, T ], H 1 (Z a ) Nm ) the solutions of the reaction-diffusion equation (9), for fixed c 0 and boundary conditions . Let

G N = ( l∈T γ 1l (c N ) ν N l * ψ a (•), • • • , l∈T γ Nml (c N ) ν N l * ψ a (•)) ∈ W 0,0 L 2 (Za) Nm and G = ( l∈T γ 1l (c) ξ l * ψ a (•), • • • , l∈T γ Nml (c) ξ l * ψ a (•)) ∈ W 0,0 L 2 (Za)

Nm

.

By the continuous injection of W 2,1

L 2 (Za) Nm in C([0, T ], H 1 (Z a ) Nm
) and the proposition 35 applied to

u = G N -G, we get that sup t∈[0,T ] ∥c N (t) -c(t)∥ (L 2 (Za)) Nm ≤ sup t∈[0,T ] ∥c N (t) -c(t)∥ H 1 (Za) Nm ≤ β 1 ∥c N -c∥ W 2,1 L 2 (Za ) Nm (32) and ∥c N -c∥ W 2,1 L 2 (Za ) Nm ≤ β 2 ∥G N -G∥ W 0,0 L 2 (Za ) Nm . ( 33
)
Then, using the same estimates as in Prop. 10 and lemma 11,

∥G N -G∥ W 0,0 L 2 (Za ) Nm ≤ Nm m=1 T 0 ∥ l∈T γ ml (c N (t)) ν N l t * ψ a -γ ml (c(t)) ξ l * ψ a ∥ 2 L 2 (Za) dt 1/2 ≤ N m γ ∞ l∈T T 0 ∥(ν lN t -ξ l t ) * ψ a ∥ L 2 (Za) dt + ||ψ a || ∞ ||γ|| Lip C(T, ξ) T 0 ∥c N (t) -c(t)∥ (L 2 (Za)) Nm . ( 34
)
with, for all ν ∈ D([0, T ], 33), (34) and applying Gronwall lemma, we obtain that

M F + (X )), C(T, ν) = sup [0,T ] ||ν t || T V < ∞. Combining (32), (
C := N m β 1 β 2 ||ψ a || ∞ ||γ|| Lip C(T, ξ) is a positive constant independent of N such that sup t∈[0,T ] ∥c N (t) -c(t)∥ H 1 (Za) Nm ≤ N m β 1 β 2 γ ∞ e CT l∈T T 0 ∥(ν lN t -ξ l t ) * ψ a ∥ L 2 (Za) dt. ( 35
)
From proposition 29, for almost every t ∈ [0, T ],

d P R (ν N t , ξ t ) → n 0. So, for almost every t, d P R (ν lN t , ξ l t ) → n 0.
Then, from proposition 34 (to be perfectly rigorous: replace Z by Z a in the statement of this proposition), for almost every t,

sup z∈Za |(ν lN t -ξ l t ) * ψ a (z)| → n 0 and so ∥(ν lN t -ξ l t ) * ψ a ∥ L 2 (Za) → n 0.
From the dominated convergence theorem, we conclude that the right hand side of (35) goes to 0 when N → ∞, which shows the desired continuity result.

From lemma 22 and Prokhorov theorem, we deduce that we can find a subsequence (ν N (n) ) n of (ν N ) N ∈N * which converges weakly in D([0, T ] , M F + (X )) to a limit denoted ν. Thanks to lemma 23 and the continuous mapping theorem, it is actually the couple (ν N (n) , c N (n) ) which converges weakly in the product space to (ν, c), c being the solution of (9) associated with ν.

Before proving that the limit ν is solution of (25), we need a second technical lemma.

Lemma 24. Let ν ∈ C([0, T ], M F + (X )), then for any sequence (ν n ) n∈N of D([0, T ], M F + (X )) such that lim n→∞ ν n = ν, we have that lim n→∞ Ψ t (ν n ) = Ψ t (ν). Proof. Let C(T, ν) = sup [0,T ] ⟨ν t , 1⟩. Under the hypothesis that lim n→∞ ν n = ν with ν ∈ C([0, T ], M F + (X )),
we know that the convergence is uniform in time (Prop. 28). So, for all ϵ > 0, there exists N ϵ ∈ N such that for all n ≥ N ϵ , sup

t∈[0,T ] d P R (ν n t , ν t ) ≤ ϵ (36)
As a consequence of the definition of d P R (see ( 62)), we have

sup t∈[0,T ] ⟨ν n t , 1⟩ ≤ sup t∈[0,T ] ⟨ν t , 1⟩ + ϵ ≤ C(T, ν) + ϵ (37)
Therefore, we restrict ourselves without loss of generality to sequences such that, for all n ∈ N,

sup t∈[0,T ] ⟨ν n t , 1⟩ ≤ 2C(T, ν). ( 38
)
Moreover, uniform convergence of (ν n ) n∈N to ν implies that for all t ≤ T , the projection at t π t : ν → ν t from D([0, T ] , M F + (X )) to M F + (X ) at every element ν continuous in time. Prokhorov distance metrizes the weak topology on M F + (X ), that is to say that for all continuous bounded function f

∈ C b (X , R), the mapping ν → ⟨ν, f ⟩ is continuous on M F + (X ). So ν → ⟨π t (ν), f ⟩ is continuous at each element of C([0, T ] , M F + (X )).
We can now show the continuity of Ψ t at ν. We denote c n = Φ(ν n ) and c = Φ(ν), then

Ψ t (ν n ) -Ψ t (ν) = ⟨ν n t -ν t , f t ⟩ -⟨ν n 0 -ν 0 , f 0 ⟩ - t 0 ⟨ν n s -ν s , ∂ s f s ⟩ds (An) - t 0 ⟨ν n s , F(•, ν n s * F (•))∇ z f s ⟩ds - t 0 ⟨ν s , F(•, ν s * F (•))∇ z f s ⟩ds (Bn) - k∈E t 0 X ⟨µ k (x), fs⟩q k (x, ν n s * D(z), c n s * ψa(z))ν n s (dx)ds - X ⟨µ k (x), fs⟩q k (x, νs * D(z), cs * ψa(z))νs(dx) ds (Cn)
.

According to the foregoing, regularity of f and compactness of X , we have

⟨ν n t -ν t , f t ⟩ → n→+∞ 0, ⟨ν n 0 -ν 0 , f 0 ⟩ → n→+∞ 0 and s → ⟨ν n s -ν s , ∂ s f s ⟩ converges pointwise to 0. Moreover, |⟨ν n s -ν s , ∂ s f s ⟩| ≤ sup (s,x) {∂ s f s (x)} × (⟨ν n s , 1⟩ + ⟨ν s , 1⟩) ≤ 3C(T, ν) sup (s,x) {∂ s f s (x)} < +∞.
By dominated convergence, t 0 ⟨ν n s -ν s , ∂ s f s ⟩ds → n→+∞ 0. Therefore, lim n→+∞ A n = 0. We have

B n = B 1 n + B 2 n with B 1 n = t 0 ⟨ν n s -ν s , F(•, ν s * F (•))∇ z f s ⟩ds and B 2 n = t 0 ⟨ν n s , (F(•, ν n s * F (•)) -F(•, ν s * F (•)))∇ z f s ⟩ds. F(•, ν s * F (•))∇ z f s is continuous bounded as X is compact, so the integrand in B 1 n converges pointwise to 0. As for A n , by the dominated convergence theorem, B 1 n -→ n→+∞ 0. Moreover, F is Lipschitz according to Hyp. 2.1, therefore B 2 n ≤ t 0 ⟨ν n s , ||F|| Lip × |ν n s * F (•) -ν s * F (•)| × |∇ z f s |⟩ds.
As F is Lipschitz and bounded on the compact set Ẑ under Hyp. 2.1, we can apply proposition 34:

∀t ∈ [0, T ]: sup z∈Z |(ν n t -ν t ) * F (z)| → n→∞ 0. So B 2 n -→ n→+∞ 0.
Finally, let

C n = C 1 n + C 2 n with C 1 n = k∈E t 0 X ⟨µ k (x), f s ⟩q k (x, ν n s * D(z), c n s * ψ a (z)) (ν n s (dx) -ν s (dx))ds
and

C 2 n = k∈E t 0 X ⟨µ k (x), f s ⟩ (q k (x, ν s * D(z), c s * ψ a (z)) -q k (x, ν n s * D(z), c n s * ψ a (z))) ν s (dx)ds
The term C 1 n → n→+∞ 0 as q k are continuous bounded according to Hyp. 14. To tackle the term C 2 n , we observe that for all t ∈ [0, T ] and all z ∈ Z,

|(c n t -c t ) * ψ a (z)| ≤ ∥ψ a ∥ L 2 ∥c n t -c t ∥ L 2 (Za) Nm ≤ ∥ψ a ∥ L 2 sup t∈[0,T ] ∥Φ(ν n ) t -Φ(ν) t ∥ H 1 (Za) Nm -→ n→+∞ 0
where we used lemma 23. Due to Hyp. 8 on ψ a , we use proposition 34 and get

sup z∈Z |(ν n t -ν t ) * D(z)| → n→∞ 0.
Finally, using that q k is Lipschitz with respect to its second and third variable (Hyp. 14), proposition 34 for the term involving D (Lipschitz and bounded on Ẑ according to Hyp. 14) and the dominated convergence theorem, we get that C 2 n → n 0. This ends the proof of the continuity of Ψ t at ν.

Lemma 25. ν is almost surely a solution of (25).

Proof. We observe that ν is solution of (25) if and only if Ψ t (ν) = 0 a.s. for all t and f . We therefore prove the latter. We first show that

ν ∈ C([0, T ], M F + (X )) a.s. (39)
ν is right continuous by definition. Let us prove the left continuity. We use total variation norm on M F + (X ) defined by (see A p.29)

||µ -ν|| T V = sup {⟨µ -ν, f ⟩; f ∈ C(X ), ||f || ∞ ≤ 1} . Let f ∈ C(X ), ||f || ∞ ≤ 1.
For all n, random jumps of νN(n) happen one at a time a.s. so, for all t ≤ T ,

|⟨ν N (n) t , f ⟩ -⟨ν N (n) t - , f ⟩| ≤ 2 ||f || ∞ N (n) ≤ 2 N (n) a.s. ⇒ sup 0≤t≤T ||ν N (n) t - νN(n) t - || T V ≤ 2 N (n) a.s..
Total variation distance dominates over Prokhorov distance on M F + (X ) (30 p.30). Therefore,

lim n→∞ sup 0≤t≤T d P R (ν N (n) t , νN(n) t - ) = 0 a.s.
and, by theorem 10.2 p.148 in [START_REF] Stewart | Markov Processes: Characterization and Convergence[END_REF],

sup 0≤t≤T d P R (ν t , νt -) = 0 a.s.
which proves (39). We now observe that Ψ t (ν N ) = Z f,N t with Z f,N t the martingale defined in Prop. 18. Then, under Hyp. 20:

E Ψ t (ν N (n) ) 2 = E ⟨Z f,N (n) ⟩ t ≤ 1 N (n) ||f || 2 ∞ t(C(p, t) + 1) k∈E q ∞ k -→ n→∞ 0 ( 40
)
with C(p, t) a constant defined in Prop. 17. Therefore, (Ψ t (ν N (n) )) n converges in L 2 , and so in L 1 by Cauchy-Schwartz, to 0.

Given νN(n) converges weakly to ν and Ψ t is continuous at ν, we get from lemma 24 and the continuous mapping theorem that |Ψ t (ν N (n) )| converges weakly to |Ψ t (ν)|. Moreover, under Hyp. 20 and from Prop. 17.1, the sequence (Ψ t (ν N (n) )) n is bounded in L p with p ≥ 2 and so uniformly integrable. Then, we have (see theorem 3.5 in [START_REF] Billingsley | Convergence of probability measures[END_REF]) that, ∀q ≤ p,

lim n→∞ E |Ψ t (ν N (n) )| q = E [|Ψ t (ν)| q ] .
In particular,

0 = lim n→∞ E |Ψ t (ν N (n) )| = E [|Ψ t (ν|] ⇔ Ψ t (ν) = 0 as.
Note that ν is continuous in time a.s., and therefore νN 0 converges weakly to ν0 . But, under Hyp. 20, νN 0 converges in distribution to ξ 0 so ξ 0 and ν0 have the same law. This is true for any t ≥ 0, we deduce that ν is solution of (25), and the sequence of concentrations (c N (n) ) n converges as well to the solution of (26

) in C([0, T ] , H 1 (Z a ) Nm ).
So far, we showed that all accumulation points of (ν N ) N ∈N * are positive measures, solutions of (25), and related accumulation points in the concentrations' space are solutions of (26) in C([0, T ] , H 1 (Z a ) Nm ). To complete the proof of theorem 21, we tackle uniqueness of the solution.

Lemma 26. If (25-26) has a solution (ξ, c) such that ξ ∈ C([0, T ], M F + (X )) then this solution is unique.

Proof. We proceed in two steps.

Step 1 : on a finite time horizon, positive solutions of (25) have finite mass. Let (c 0 , ξ 0 ) ∈ H 1 (Z a ) Nm × M F + (X ) be a deterministic initial condition. We suppose that there exists a solution (c t , ξ

t ) t≥0 ∈ C([0, T ] , H 1 (Z a ) Nm ) × D([0, T ] , M F + (X )
). Let show that the mass of ξ t is bounded at finite time by a constant depending only of ξ 0 and t. From (25), for t ≤ T ,

⟨ξ t , 1⟩ = ⟨ξ 0 , 1⟩ + k∈E t 0 l∈T Z ⟨µ k (z, l), 1⟩q ∞ k q k ((z, l), ξ s * D(z), c s * ψ a (z))ξ s (dz, l)ds ≤ ⟨ξ 0 , 1⟩ + k∈E q ∞ k t 0 l∈T Z ξ s (dz, l)ds = ⟨ξ 0 , 1⟩ + k∈E q ∞ k t 0 ⟨ξ s , 1⟩ds. ( 41
)
By Gronwall lemma, we see that

⟨ξ t , 1⟩ ≤ ⟨ξ 0 , 1⟩ exp(t k∈E q ∞ k ) < +∞,
which is the desired result.

Step 2 : a positive solution of (25-26) is unique. Let (c 1 t , ξ 1 t ) and (c 2 t , ξ 2 t ) be two solutions of (25-26) on [0, T ], with ξ i t positive as in the previous step. We use the flat metric d F on M F + (X ) defined as (more details are given in A)

d F (ν, ξ) = sup f ∈Lip1,1(X ,R) ⟨ν -ξ, f ⟩,
and derive a bound for

||c 1 t -c 2 t || L 2 + d F (ξ 1 t , ξ 2 t ), for t ∈ [0, T ]. From step 1, we can define C(t, ξ 1 + ξ 2 ) := sup s≤t ⟨ξ 1 s + ξ 2 s , 1⟩ < +∞.
Moreover, ξ 1 t , ξ 2 t are finite positive measures so, from lemma 11, there exists a constant C such that

||c 1 t -c 2 t || L 2 ≤ C t 0 d F (ξ 1 s , ξ 2 s ) ds. ( 42
)
We now bound d F (ξ 1 t , ξ 2 t ) using formula (25). As done in [START_REF] Viet | Stochastic particle models for problems of adaptive evolution and for the approximations of statistical solutions[END_REF], the idea is to cancel in equation (25) the integral term accounting for the deterministic movement of cells between jumps. Let u ∈ Lip 1,1 ∩ C 1,• (X , R). From lemma 38, there exists a unique function

u 1 : [0, T ] × X → R (s, z, l) → u 1 s (z, l)
such that, ∀l ∈ T , u 1 s (z, l) is solution of (74), with ξ s = ξ 1 s and the condition u 1 t (z, l) = u(z, l), that is to say ∂ s u 1 s (z, l) + F(z, ξ 1 s * F (z))∇ z u 1 s (z, l) = 0 and u 1 t (x) = u(x). We also know that u 1 s (z, l) = u(Z 1 (s; t, z), l) where Z 1 (•) are the characteristics associated to (74) when ξ s = ξ 1 s . Similarly, let u 2 s (z, l) = u(Z 2 (s; t, z), l), with ξ s = ξ 2 s in (74). Moreover, still from lemma 38,

u i • (•, l) ∈ C 1 ([0, T ] × Z, R).
We then apply formula (25), with ξ 1 0 = ξ 2 0 = ξ 0 .

|⟨ξ 1 t -ξ 2 t , u⟩| = |⟨ξ 0 , u 1 0 -u 2 0 ⟩ + k∈E t 0 l∈T Z ⟨µ k (z, l), u 1 s ⟩q k ((z, l), ξ 1 s * D(z), c 1 s * ψ a (z))ξ 1 s (dz, l)ds - k∈E t 0 l∈T Z ⟨µ k (z, l), u 2 s ⟩q k ((z, l), ξ 2 s * D(z), c 2 s * ψ a (z))ξ 2 s (dz, l)ds = | ⟨ξ 0 , u 1 0 -u 2 0 ⟩ (A) + k∈E t 0 l∈T Z ⟨µ k (z, l), u 1 s -u 2 s ⟩q k ((z, l), ξ 1 s * D(z), c 1 s * ψ a (z))ξ 1 s (dz, l)ds (B) + k∈E t 0 l∈T Z ⟨µ k (z, l), u 2 s ⟩(q k (z, l), ξ 1 s * D(z), c 1 s * ψ a (z)) -q k ((z, l), ξ 2 s * D(z), c 1 s * ψ a (z)))ξ 1 s (dz, l)ds (C) + k∈E t 0 l∈T Z ⟨µ k (z, l), u 2 s ⟩(q k ((z, l), ξ 2 s * D(z), c 1 s * ψ a (z)) -q k ((z, l), ξ 2 s * D(z), c 2 s * ψ a (z)))ξ 1 s (dz, l)ds (D) + k∈E t 0 l∈T Z ⟨µ k (z, l), u 2 s ⟩q k ((z, l), ξ 2 s * D(z), c 2 s * ψ a (z))(ξ 1 s -ξ 2 s )(dz, l)ds (E) |.
Before proceeding, we recall the following simple properties that we will use extensively. Let f be Lipschitz on X , then

|| f ||f ||∨ || ∞ ≤ 1, || f ||f ||∨ || Lip ≤ 1 and, for ν 1 , ν 2 ∈ M F + (X ) ⟨ν 1 -ν 2 , f ⟩ ≤ ||f || ∨ × d F (ν 1 , ν 2 ).
We first derive a bound for (C). Remember that ⟨µ k (z, l), u 2 s ⟩ is the sum of at most two punctual evaluations of u 2 s . So ||⟨µ k (z, l), u 2 s ⟩|| ∞ ≤ 2. We supposed (Hyp. 14) that q k is Lipschitz with respect to its second variable, therefore there exists a constant C(q k ) such that

(C) ≤ 2 k∈E t 0 l∈T Z C(q k )|ξ 1 s * D(z) -ξ 2 s * D(z)|ξ 1 s (dz, l)ds.
D is locally Lipschitz in z under Hyp. 14 on the compact set Ẑ := {z -y : (z, y) ∈ Z 2 }. Therefore, for all z ∈ Z, |ξ

1 s * D(z) -ξ 2 s * D(z)| ≤ ||D a| Ẑ || ∨ d F (ξ 1 s , ξ 2 s ). Let C (C) := 2 k∈E C(q k )C(t, ξ 1 + ξ 2 )||D a| Ẑ || ∨ . Then (C) ≤ C (C) t 0 d F (ξ 1 s , ξ 2 s )ds.
We find an upper bound for the term (D) using again the fact that q k is Lipschitz with respect to its third variable:

(D) ≤ 2 k∈E C(q k ) t 0 l∈T Z |(c 1 s -c 2 s ) * ψ a (z)| ξ 1 s (dz, l)ds ≤ 2 k∈E C(q k )||ψ a || L 2 C(t, ξ 1 + ξ 2 ) t 0 ||c 1 s -c 2 s || L 2 (Za) Nm ds ≤ C (D) t 0 ||c 1 s -c 2 s || L 2 (Za) Nm ds setting C (D) := 2 k∈E C(q k )||ψ a || L 2 (R d ) C(t, ξ 1 + ξ 2 ).
We move on to term (E). For s ≤ t, let f s : (z, l) → ⟨µ k (z, l), u 2 s ⟩q k (z, l, ξ 2 s * D(z), c 2 s * ψ a (z)) be defined on X . Let us show that f s is uniformly bounded and uniformly Lipschitz for s ≤ t.

On the one hand, q k and u 2 are bounded, so sup s≤t ||f s || ∞ < +∞.

On the other hand, ∀l ∈ T , z → q k (z, l, ξ

(•, l, • • • ) is Lipschitz under Hyp. 14 so z → q k ((z, l), ξ 2 s * D(z), c 2 s * ψ a (z)) is Lipschitz uniformly for s ≤ t from Prop. 1.
Moreover, for all l, z → u 2 s (z, l) is Lipschitz uniformly in time as it is continuously derivable, so u 2 s is Lipschitz on X uniformly in time. Then, when the jump type k is a differentiation or an extrusion, ⟨µ k (z, l), u 2 s ⟩ is also Lipschitz uniformly in time. If k is a division, it is the case as well because λ(z) is Lipschitz according to Hyp. 12. Eventually, the product of two Lipschitz bounded functions is Lipschitz. Therefore, f s is Lipschitz uniformly for

s ≤ t. Let C (E) := sup s≤t ||f s || ∨ < +∞. Then, (E) ≤ C (E) t 0 d F (ξ 1 s , ξ 2 s ) ds.
In order to bound the remaining terms (A) and (B), we want an upper bound on

|u 1 s (z) -u 2 s (z)| depending of d F (ξ 1 s , ξ 2 s ). We have |u 1 s (z) -u 2 s (z)| = |u(Z 1 (t; s, z)) -u(Z 2 (t; s, z))| ≤ ||u|| ∨ |Z 1 (t; s, z) -Z 2 (t; s, z)| ≤ |Z 1 (t; s, z) -Z 2 (t; s, z)| since u ∈ Lip 1,1 (X , R).
Using proposition 39, there exists a constant C(t) such that

(A) = l∈T Z u(Z 1 (t; 0, z), l) -u(Z 2 (t; 0, z), l) ξ 0 (dz, l) ≤ l∈T Z |Z 1 (t; 0, z) -Z 2 (t; 0, z)| ξ 0 (dz, l) ≤ ⟨ξ 0 , 1⟩C(t) t 0 d F (ξ 1 s , ξ 2 s ) ds ≤ C (A) t 0 d F (ξ 1 s , ξ 2 s ) ds
with C (A) := ⟨ξ 0 , 1⟩C(t). Note that the first inequality comes from the assumption that u is 1-Lipschitz, which motivates the use of the flat metric d F .

With proposition 39 and as u is 1-Lipschitz, we can also tackle term (B):

(B) ≤ k∈E 2q ∞ k C(t, ξ 1 + ξ 2 )C(t) t 0 t s d F (ξ 1 r , ξ 2 r )drds ≤ k∈E 2q ∞ k C(t, ξ 1 + ξ 2 )C(t)t t 0 d F (ξ 1 s , ξ 2 s ) ds ≤ C (B) t 0 d F (ξ 1 s , ξ 2 s ) ds with C (B) := k∈E 2q ∞ k C(t, ξ 1 + ξ 2 C(t)t. Therefore, we showed that there exists a constant C such that, ∀u ∈ Lip 1,1 (X , R) ∩ C 1,• (X , R), |⟨ξ 1 t -ξ 2 t , u⟩| ≤ C t 0 ||c 1 s -c 2 s || L 2 (Z) + d F (ξ 1 s , ξ 2 s ) ds.
As X is compact, according to Stone-Weierstrass theorem, the set of functions C

1,• (X , R) is dense in C(X , R) for the uniform norm. So sup u∈Lip1,1∩C 1,• |⟨ξ 1 t -ξ 2 t , u⟩| = sup u∈Lip1,1 |⟨ξ 1 t -ξ 2 t , u⟩| = d F (ξ 1 t , ξ 2 t )
and

d F (ξ 1 t , ξ 2 t ) ≤ C t 0 d F (ξ 1 s , ξ 2 s ) + ||c 1 s -c 2 s || L 2 (Z) ds. ( 43
)
Using (42), we obtain that there exists a constant C such that

||c 1 t -c 2 t || L 2 + d F (ξ 1 t , ξ 2 t ) ≤ C t 0 ||c 1 s -c 2 s || L 2 (Z) + d F (ξ 1 s , ξ 2 s ) ds. ( 44
)
From Gronwall lemma,

||c 1 t -c 2 t || L 2 + d F (ξ 1 t , ξ 2 t )
= 0 which proves the uniqueness of a positive solution.

This shows that the sequence (ν N ) N ∈N * is tight with a unique accumulation point, which is solution of (25), so it is convergent and the sequence of concentrations converges to the corresponding solution of (26) in C([0, T ] , H 1 (Z a ) Nm ). This ends the proof of theorem 21.

6 Numerical convergence of the PDMP towards a deterministic limit

Description of the PDMP model

We illustrate our result on a simple model inspired from a more complex intestinal crypt model previously introduced [START_REF] Darrigade | A PDMP model of the epithelial cell turn-over in the intestinal crypt including microbiota-derived regulations[END_REF]. We first set the precise numerical experiments that are conducted, considering two populations of different types and stochastic regulation processes including cell division, extrusion and differentiation. Then, we introduce the corresponding deterministic limit in large population for this model. We consider a one dimensional spatial model on Z = [0, z max ], 0 < z max < ∞, with two different type of cells: progenitor cell (pc) and enterocyte cells (ent) that is to say T = {pc, ent}, and two metabolites: butyrate and oxygen (N m = 2). We denote c b (z, t) and c o (z, t) the respective concentration of butyrate and oxygen. The size of the cells is modelled through a parameter a > 0.

Mechanical interactions for cell migration. We set

dz i dt (t) = F( z i (t), ν t * F a (z i (t)) ) := ϕ(z i (t)) ν t * F a (z i (t)) ( 45
)
with

F a (z) = k a f a (z) ||f a || ∞ where f a (z) =    z a exp - |z| a 1- |z| a if |z| ≤ a, 0 otherwise. ( 46
)
k a is a scaling parameter driving the maximal intensity of the force. F a is plotted in Fig. 1a with

k a = 1. ϕ(z) is a positive continuously derivable function, such that 0 ≤ ϕ(z) ≤ 1, ϕ(0) = ϕ(z max ) = 0 and ϕ(z) = 1 if z ∈ [ϵ, z max -ϵ].
The function ϕ models the vanishing of interaction forces near the boundaries, induced by geometrical constraints (see [START_REF] Darrigade | A PDMP model of the epithelial cell turn-over in the intestinal crypt including microbiota-derived regulations[END_REF] for more details). The exact expression of ϕ is given in equation ( 78).

Reaction-diffusion of chemicals.

The oxygen and butyrate concentrations c o and c b follow the reaction-diffusion equation ( 9), leading to the system:

∂ t c o -σ∆c o = l∈{pc,ent} γ ml (c o , c b ) ν l t * ψ a (z), ( 47
)
∂ t c b -σ∆c b = l∈{pc,ent} γ ml (c o , c b ) ν l t * ψ a (z). ( 48
)
Progenitor cells do not use oxygen and butyrate for their metablism, while enterocytes catalyze them with the following chemical reaction 1 butyrate + 4 O 2 -→ energy + byproducts.

Therefore we set for the reaction speed in equation ( 47

γ m,l (c o , c b ) = -s m γ ∞ c 4 o c b c 4 o c b + K 5 (49)
where s m is the stoechiometric coefficient of metabolite m. γ ∞ is the maximal reaction speed of the reaction and K > 0 is a constant. We define ψ a (z) as follow. Let

D a (z) := d a D z a , where D(z) := exp -z 2 1-z 2 if 0 ≤ |z| ≤ 1 0 otherwise, ( 50
)
and d a is a scaling parameter (see Fig. 1a for a plot with d a = 1). We then derive a normalized version of the density kernel D a by writing

ψ a (z) := D a/2 (z) a/2 -a/2 D a/2 (u)du . ( 51
)
Then, for all z ∈ [0, z max ], Za ψ a (z -y)dy = 1. We see with this kernel definition that the extended domain Z a has been introduced for compatibility reasons between the reaction-diffusion domain and the epithelial cell domain [0, z max ]. This system of equations is supplemented with boundary conditions modeling the chemical fluxes at the crypt extremities. For butyrate:

c b = c b,lum at z = z max + a 2 ∂ z c b = 0 at z = -a 2 (52)
and for oxygen:

∂ z c o = 0 at z = z max + a 2 c o = c o,bot at z = -a 2 (53)
Constant parameters c i,lum , c i,bot represent the chemical concentrations in the gut lumen or at the bottom of the crypt, respectively reflecting microbial or host physiology.

Stochastic jumps modelling cell fate. We consider three types of stochastic events. Progenitor cells can divide symmetrically or differentiate into enterocytes, and enterocytes can die. As explained in section 3.2, these stochastic events are regulated in various ways, which are summarized in table 1:

• the rate of progenitor cell division q pc,pc (•) is a decreasing function of local butyrate concentration, height in the crypt (z ∈ [0, z max ]) and local cell density. Upon division, the spatial shift of the daughter cell λ is

λ(z) =      a 2 if z ≤ z max /2 -a, -a 2 if z ≥ z max /2 + a, a 2 -a z-(zmax/2-a) 2a otherwise,
• The rate of progenitor cell differentiation q pc,ent (•) is positively correlated to height in the crypt.

• The rate of enterocyte death q ent,∅ (•) is positively correlated to height in the crypt and local cell density.

More precisely, we define the local density of the population d as

d(ν t , z) = d(z)ν t * D a (z) (54)
where D a is defined in Eq.( 50) and

d(z) = a -a D a (u) du min(a,zmax-z) max(-a,-z) D a (u) du (55)
is a rescaling function ensuring that the density kernel sums to 1 even near the boundaries. As a << z max , 1 ≤ d(z) ≤ 2. We also define the local concentration of butyrate c b as

c b (c b , z) = c b * ψ a (z) (56)
where ψ a (z) is introduced in Eq. (51). Eventually, for a cell at x := (z, l) ∈ X , in a population ν ∈ M F + (X ) and for butyrate and oxygen concentration c := (c o , c b ), we define the individual jump rate of progenitor cell division by

q pc,pc (x, ν * D(z), c * ψ a (z)) = 0 if not prog. cell 1{l = pc} × inhibition by z (1 -R(z, K pc,pc [z], κ pc,pc [z])) × inhibition by density (1 -R(d(ν t , z), K pc,pc [dens], κ pc,pc [dens]) × (1 -R(c b (c b , z), K pc,pc [but], κ pc,pc [but])) inhibition by butyrate (57)
with K pc,pc [•] and κ pc,pc [•] strictly positive constants and R(y, K, κ) is the piecewise polynomial function (plotted on Fig. 1b for several values of K, and κ) defined by

R(y, K, κ) :=      0 if y ≤ K -κ, -1 4κ 3 y 3 + 3K 4κ 3 y 2 -3K 2 -3κ 2 4l 3 y + K 3 +2κ 3 -3Kκ 2 4κ 3 if K -κ ≤ y ≤ K + κ, 1 if K + κ ≥ y.
Note that R connects continuously the values 0 and 1 for K -κ < z < K + κ and has value 1/2 in y = K.

We see that the jump rate depends on a product of three regulation terms involving the function R that corresponds to the three biological regulation of progenitor cell division given in table 1. Using the same construction, we define functions q pc,ent and q ent,∅ based on biological regulations given in table 1. For a positive regulation, we use a term R(•) instead of 1 -R(•). And no regulation correspond to a constant term equal to 1.

The values of parameters used for simulation are given in table 2.

Regulation pathway (j)

Jump type Jump index (k ∈ E) butyrate [but] position [z] density [dens]

Progenitor cell division pc, pc

---

Progenitor cell differentiation pc, ent

+ + ∅ Enterocyte death ent, ∅ ∅ + +
Table 1: Enumeration of stochastic jump types and their corresponding regulation. For the jump k and a given regulation pathway j, a negative regulation (-) is represented by the term

(1 -R(y, K k [j], κ k [j]
)), a positive regulation (+) by the term R(y, K k [j], κ k [j]) and a neutral regulation (∅) by a constant value equal to 1, where y = z, y = d(ν t , z) or y = c b (c b , z) for respectively position, density or butyrate regulation. The values of K k [j] and κ k [j] can be found in table 2.

Initial conditions

Initial conditions for the position of the cells in the stochastic process are sampled in distributions ρ pc and ρ ent defined in eqs. ( 79)-( 80), in D with a Monte-Carlo sampling, hence guarantying the weak convergence of the rescaled initial condition towards ρ pc and ρ ent (see next section). Furthermore, initial conditions for chemicals are defined : c b (0, z) = c b0 (z) and c o (0, z) = c o0 (z), for smooth c b0 and c o0 functions defined in eq. ( 81).

Rescaled parameters. To illustrate the convergence of the stochastic model, we used the following rescaled parameters:

F N a = F a N , D N a = D a N , γ N = γ N .
where N is a multiplicative factor applied to the progenitor and enterocyte cell populations included in the initial condition.

Deterministic model obtained as a limit of the PDMP: SA equation

Since the PDMP model definition complies with the hypothesis 2, 7, 8, 9, 12, 14 and 20, according to Th. 21 and formula 25, the deterministic limit ξ

(SA) t
of the PDMP model is the weak solution of the following system, here written in the strong form:

∂ t ξ (SA) t + ∇ z F(•, ξ (SA) t * F a )ξ (SA) t = k∈E l∈T Z µ k (z, l)q ∞ k q k ((z, l), ξ (SA) t * D a (z), c t * ψ a (z))ξ (SA) t (dz, l). ( 58
)
This equation can be interpreted as a non-linear advection-reaction equation, that will be termed selfaggregation (SA) equation in the sequel. Note that the non linear transport term ∇ z (F(•, ξ t * F a )ξ t ) comes from the cell interaction model of the deterministic flow of the PDMP model. The reaction term is a deterministic approximation of the stochastic terms of the PDMP.

For the metabolites m ∈ {but, oxy}, the same PDE equations as in the PDMP model are kept, except for the reaction term that now includes the density ξ

(SA) t ∂ t c (SA) m -σ∆c (SA) m = l∈T γ ml (c (SA) ) ξ (SA) t l * ψ a (z). ( 59
)
This model is supplemented by a deterministic version of the boundary conditions

F(•, ξ (SA) t * F )ξ (SA) t (0, t) = 0 c (SA) m (z, t) = f m (z, t) for z ∈ Γ, ξ (SA) 0 (•, 0) = (ρ pc , ρ ent ) c (SA) m (•, 0) = c 0m ,

Numerical implementation

Simulations are performed with dedicated Python scripts available here https://forgemia.inra.fr/ leo.darrigade/cvg. To solve the PDMP, we use a standard Gillespie algorithm with an acceptationreject method for stochastic events, while the deterministic flows are solved with an explicit Euler method for the cell movement, and with a semi-implicit Euler scheme with operator splitting for the reactiondiffusion equations on butyrate and oxygen. To solve the deterministic limit, we use a semi-implicit scheme for the non-linear diffusion term and a Runge-Kutta4 scheme for time integration. Scripts use NumPy [START_REF] Harris | Array programming with NumPy[END_REF], Scipy [START_REF] Virtanen | SciPy 1.0: Fundamental Algorithms for Scientific Computing in Python[END_REF] sparse algebra and Numba [START_REF] Siu Kwan Lam | Numba: A LLVM-based Python JIT Compiler[END_REF] packages for additional code acceleration.

Numerical assessment of the convergence

As the deterministic limit of cell densities ξ (SA) t is actually continuous in time (see proof of lemma 25), the weak convergence of νN t to ξ (SA) t at a fixed time t is true for almost every t. Then, for a continuous and bounded function f : X → R, ⟨ν N t , f ⟩ converges weakly to ⟨ξ

(SA) t , f ⟩.
From uniform integrability of the sequence (⟨ν N t , f ⟩) N , as in the proof of lemma 25, we get that

E ⟨ν N t , f ⟩ → N ⟨ξ (SA) t , f ⟩. (60)
We now illustrate this convergence using some specific functions f .

Total populations

First, taking f = 1, the quantity E ⟨ν N t , f ⟩ and ⟨ξ

(SA) t
, f ⟩ represent the total cell population at time t in the PDMP and the SA equation, respectively. Fig. 2 displays the evolution of the total populations in the PDMP (mean and quantiles 2.5% and 97.5% for n r = 30 replicates) when the population scaling N = 60. We can observe that the progenitor cells of the SA equation is in the range of the PDMP outputs, while the total number of enterocytes is over-estimated by the SA deterministic limit by 10% at t = 29h. This discrepancy could be due to the effect of death events occuring in the upper part of the crypt near z max . This results in lower cell populations in this region, for which the convergence towards the deterministic limit is very slow.

The inacurracies in total enterocyte populations are propagated to the evaluation of the total number of cells (progenitor and enterocyte cells).

Convergence for growing N

Let us first introduce, for (z a , z 2 ) ∈ (0, z max ) 2 and a scaling parameter α > 0 the kernel

φ(z 1 , z 2 |α) = exp( -1 1-( |z 1 -z 2 | α ) 2 ) if |z 1 -z 2 | < α 0 otherwise
and, for z ∈ (0, z max , we define the function ϕ z indexed by z so that ϕ z (x) := φ(z, x) for any x ∈ (0, z max ).

The functions ϕ z allow to test the PDMP and SA solutions in different zones of the crypt. We finally define the residual for cell populations with

R N (t, z) := 1 2 l∈{pc,ent} E ⟨ν N l,t -ξ (SA) l,t , ϕ z ⟩ (61)
where the true expectation E [•] is approximated by the mean over n r replicates of the PDMP, and νN l,t = 1 N ν N l,t is the rescaled solution of the PDMP for the rescaling parameter N and the population l ∈ {pc, ent} at time t > 0.

According to the convergence theorem, we must have for any t > 0 and z ∈ (0, z max )

R N (z, t) → N 0
We display the curves N → R N (t, z) (N ∈ {15, 20, 25, 30, 35, 40, 45, 50, 55, 60}) for several t (8 h, 16 h, 24 h, 29 h) and several z (34 µm, 42 µm, 59 µm, 68 µm), either grouped by t (Figure 3) or by z (Figure 4). In all curves, we observe a consistent reduction in residuals between the PDMP and SA equation as N increases, despite some local fluctuations. The convergence rate is relatively slow : an average calculation shows that, with a fourfold increase in the population scaling parameter (from N = 15 to N = 60), residuals decrease by an average factor of 0.86.

We note that at a given time, this decrease follows a common pattern for every z ∈ {34, 42, 59, 68}, since the curves are nearly parallel in Figure 3. However, convergence seems to be faster for higher z, since red curves (z = 68) and green curves (z = 59) have the steepest decrease in the four panels of Figure 3. At a given z (Figure 4, the decrease of residuals tends to be more stochastic than for the curves grouped by times (Figure 3). The convergence rate is lower at the end of the simulation, for t = 29 h, but no general trend can be deciphered for earlier times.

The residual R N (z, t) is a metric that only focuses on cell populations. We now look at the convergence of the PDMP towards the SA equation for the metabolite concentrations. In Figure 5 computed with the deterministic SA equation. We can see that the concentration maps are very similar whatever the value of N . The SA equation renders very well the butyrate distribution, which has the main impact on the cell population dynamics through butyrate regulation of growth rates (see eq.( 57)). Comparatively, discrepancies are higher for oxygen distribution, specially at the top of the crypt (z = 100 µm), probably for the same reasons as above: the very slow convergence towards the asymptotic cell density in this region, with an underestimation of the enterocytes, is reflected in the underestimation of the oxygen consumption. Note that the Dirichlet condition at the top of the crypt renders this effect invisible for the butyrate concentration. The deterministic self-aggregation limit is plotted with thick dots.

Conclusive discussion

In this study, we introduced a generic model of epithelial tissues interacting with diffusive compounds, in a rigorous PDMP framework. We determined the deterministic limit of the PDMP model, which is a self-aggregation like equation. We then proved the convergence of the IBM model towards this limit under a large population assumption. This convergence is finally checked numerically by solving both the deterministic limit model and the PDMP model for several population scaling parameters. This work is of special interest for modelers of epithelial tissues. It effectively links a PDMP IBM and its deterministic, continuous population counterpart, combining the strength of both modeling frameworks for epithelial tissues modeling. Micro-scale data like microscopy imaging can be used to calibrate the IBM model, while the deterministic limit allows upscaling this calibration to the continuous population model. Conversely, continuous population model parameters inferred from macro-scale data like metabolomics can be down-scaled to the IBM. Following the same idea, a PDMP IBM model can be used to model accurately intra-specific distribution in cell populations, but the deterministic limit would be more appropriate for intensive numerical exploration such as sensitivity analysis. Such an approach has been previously used to analyse a PDMP model of the cell turn-over in the intestinal crypt [START_REF] Darrigade | A PDMP model of the epithelial cell turn-over in the intestinal crypt including microbiota-derived regulations[END_REF].

The main achievement of this work is the rigorous derivation of the deterministic limit and the proof of the convergence of the PDMP model towards this limit. Although the general functional analysis framework is well established, in particular in the theoretical ecology and evolution community [FM04; CFM07; Tra06], technical issues had to be overcome in the specific context of epithelial tissues. Notably, dealing with cell-to-cell interactions introduced non-linearities in the deterministic limit, and complex regulatory processes in epithelial cell turnover added technical difficulties in proofs.

The main limitation of our study is the lack of convergence rate estimates. We could check numerically that in our test case, this convergence rate is relatively slow. Consequently, we have no guarantee that the number of cells observed in real tissues is in the convergence regime. Most probably, at the crypt scale, cell numbers might be too low for stochastic events to approximate the "average process" accurately, suggesting the importance to retain stochasticity at this scale to model population distribution around the mean. However, at the tissue scale, with numerous crypts, the large population assumption could be more valid, making the deterministic limit accurate. In conclusion, modelers should not avoid a discussion of the number of cells in the tissue of interest before choosing between a PDMP or a continuous population model, bearing in mind that connections can be made between these formalisms.

A Metrics on M

F + (X ), properties of D([0, T ] , (M F + (X ), d P R ))
We remind the reader that (X , d), with d defined at (2) is a complete separable metric space (Polish space). For a closed set F of X and ϵ ≥ 0, let F ϵ = {x ∈ X : d(x, F ) < ϵ}. The Prokhorov distance d P R on M F + (X ) is defined by (Appendix A2.5 of [DVon])

d P R (µ, ν) = inf {ϵ ≥ 0 : for all closed F ⊂ X , µ(F ) ≤ ν(F ϵ ) + ϵ and ν(F ) ≤ µ(F ϵ ) + ϵ} . (62)
The first important property of the Prokhorov distance is that it metrizes the weak convergence topology on M F + (X ) (Corollary A2.5.II in [DVon]) and this is why we use it. Another useful property is that, if (X, d) is Polish, the space (M F + (X ), d P R ) is Polish too (Proposition A.2.5.III of [DVon]). On M F (X ), we also define the total variation norm:

||µ|| T V = µ + (X ) + µ -(X ) = sup f measurable,||f ||∞≤1 X f dµ ( 63
)
where the couple µ + , µ -is the Hahn-Jordan decomposition of µ.

Eventually, we will use the flat distance (or Kantorovtich-Rubinstein distance) defined on

M F + (X ) as d F (ν, ξ) = sup f ∈Lip1,1(X ,R) ⟨ν -ξ, f ⟩. (64)
We now remind (and prove some) useful properties on these metrics.

Proposition 27 (Theorem 1.2 p.96 in [START_REF] Stewart | Markov Processes: Characterization and Convergence[END_REF]). Let µ and ν be two probability measures on X . Let Π(µ, ν) be the set of probability measures on X × X with marginal distributions µ and ν. Then

d P R (µ, ν) = inf γ∈Π(µ,ν) inf ϵ > 0 : γ{(x, y) ∈ X 2 : d(x, y) ≥ ϵ} ≤ ϵ . ( 65
)
Once a metric is defined on M F + (X ), we can define the Skorokhod metric (see [START_REF] Billingsley | Convergence of probability measures[END_REF]) on the space D([0, T ] , (M F + (X ), d P R )) of càdlàg functions of [0, T ] to M F + (X ). The space D([0, T ] , (M F + (X ), d P R )) equipped with the Skorokhod metric inherits the Polish space property from (M F + (X ), d P R ). Polish spaces are well suited to study weak convergence of random variable.

Proposition 28. Let (ν n ) n∈N be a sequence which converges to ν in D([0, T ], M F + (X )). We suppose that ν ∈ C([0, T ], M F + (X )). Then, (ν n ) n converges to ν uniformly in time:

sup 0≤t≤T d P R (ν n t , ν t ) → n 0.
Proof. From the definition of the Skorokhod metric (Prop. But, when n → ∞, the first term of the right hand side goes to 0 from the definition of convergence in the Skorokhod space, and the second term goes also to 0 from the uniform continuity in time of ν on [0, T ].

Proposition 29. Let (ν n ) n∈N be a sequence which converges to ν in D([0, T ], M F + (X )). Then, (ν n t ) n converges to ν t for almost every t ∈ [0, T ], that is to say

d P R (ν n t , ν t ) → n 0 t -a.e..
Proof. ν has at most a countable set of time discontinuities, written A (lemma 5.1 in [START_REF] Stewart | Markov Processes: Characterization and Convergence[END_REF]). Then, by definition of the Skorokhod metric (see Prop. 28), ∀t 0 ∈ [0, T ] \ A,

d P R (ν n t0 , ν t0 ) ≤ d P R (ν n t0 , ν λn(t0) ) + d P R (ν λn(t0) , ν t0 ) → n 0.
Indeed, the first term of the right hand side goes to 0 by definition of the convergence in Skorokhod distance, and the second as well as ν t is continuous at t 0 .

Proposition 30. The total variation distance is greater than the Prokhorov distance on M F + (X ).

Proof. For all µ ∈ M F (X ), we admit that (and this can be simply proven using the Hahn-Jordan decomposition of µ):

||µ|| T V ≤ 2 sup B∈B(X ) {|µ(B)|} ≤ 2||µ|| T V , (66)
with B(X ) the Borel sets. Let F be a closed set of X , µ and ν ∈ M F + (X ),

µ(F ) = µ(F ) -ν(F ) + ν(F ) ≤ sup B∈B(X ) {|µ(B) -ν(B)|} + ν(F ) ≤ 1 2 ||µ -ν|| T V + ν(F 1 2 ||µ-ν|| T V )
where we used the fact that ν is a positive measure. We can switch µ and ν in those expressions and we obtain that

d P R (µ, ν) ≤ 1 2 ||µ -ν|| T V . Proposition 31. Let ν ∈ M F (X ), ||ν|| T V = sup {⟨ν, f ⟩; f ∈ C(X ), ||f || ∞ ≤ 1} . ( 67
)
The first term of the right hand side goes to 0 by proposition 33. The second term of the right hand side goes to 0 by proposition 32. Let ϵ such that 0 < ϵ. From what was shown above, there exists

N ϵ such that, if n ≥ N ϵ , d P R ν ||ν|| T V , ν n ||ν n || T V ≤ ϵ and |1 - ||ν n || T V ||ν|| T V | ≤ ϵ. Let n ≥ N ϵ . Let G + (z) and G -(z) be two positive functions such that G(z) = G + (z) -G -(z). |(ν n -ν) * G(z)| ≤ |(ν n -ν) * G + (z)| + |(ν n -ν) * G -(z)|.
Then,

|(ν -ν n ) * G + (z)| = ||ν|| T V × ν ||ν|| T V - ||ν n || T V ||ν|| T V × ν n ||ν n || T V * G + (z) ≤ ||ν|| T V × ν ||ν|| T V - ν n ||ν n || T V * G + (z) + ||ν|| T V × (1 - ||ν n || T V ||ν|| T V ) × ν n ||ν n || T V * G + (z) ≤ ||ν|| T V × ν ||ν|| T V - ν n ||ν n || T V * G + (z) + ||ν|| T V ϵ||G + || ∞ Let νn := ν n ||ν n || T V and ν := ν ||ν|| T V
, which are probability measures on X . From proposition 27 p.29, ∀δ such that ϵ < δ, there exists a probability measure on X × X γ ∈ Π(ν n , ν) and a set E ⊃ {(x, y) ∈ X 2 : d(x, y) ≥ δ} such that γ(E) ≤ δ. Then We can moreover assume that ϵ < δ < diam(Z). Then, for all x 1 = (z 1 , l 1 ), x 2 = (z 2 , l 2 ) ∈ X 2 , if l 1 ̸ = l 2 then (x 1 , x 2 ) ∈ E. Thus, we can use the fact that G + is Lipschitz in space to find an upper bound for the last integral on the domain X \ E: This concludes the proof as ϵ → 0.

B 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 35 which states existence and uniqueness of solution -λc + Ac = 0 (71)

∂ n c |∂U + bc |∂U = f (72)
where λ ≥ 0 is a real number greater than the largest spectral value of A. For all α ≥ 0, the problem (71) has a unique solution c in H 2α (U ) when f ∈ H 2α-3 2 (∂U ), denoted c = D(f ). Operator D is bounded from H 2α-3 2 (∂U ) to H 2α (U ). These standard results can be found in [START_REF] Louis | Non-homogeneous boundary value problems and applications[END_REF] or [START_REF] Amann | Nonhomogeneous linear and quasilinear elliptic and parabolic boundary value problems[END_REF] . Using notations taken from [START_REF] Lasiecka | Unified theory for abstract parabolic boundary problems-a semigroup approach[END_REF] , we define the following functional spaces, for all α ≥ 0 and for L 2 (U ) = H 0 (U ),

W 2α,⌈α⌉ L 2 (U ) = f ∈ L 2 ([0, T ], H 2α (U )); ∀j ∈ {1, • • • , ⌈α⌉} ∂ j f ∂t j ∈ L 2 ([0, T ], H 2α-2j (U ))
where ⌈•⌉ is the ceiling function. W 2α,⌈α⌉ L 2 (U ) is a Hilbert space endowed with the norm

  ∥f ∥ 2 L 2 ([0,T ],H 2α (U )) + ⌈α⌉ j=1 ∥ ∂ j ∂t j f ∥ 2 L 2 ([0,T ],H 2α-2j (U ))   1 2
.

Note that with these notations, W 0,0 L 2 (U ) = L 2 ([0, T ] × U ). We define in the same way the space W 2α,⌈α⌉ L 2 (∂U ) . From [LM72] theorem 3.1 p. 19, we know that: This solution belongs to W 2,1 L 2 (U ) ∩ C([0, T ], H 1 (U )). Moreover, the operator (c 0 , f, u) → c is continuous from

• if f ∈ W
H 1 (U ) × W 1 2 ,1 L 2 (∂U ) × W 0,0 L 2 (U ) to W 2,1 L 2 (U ) .
Proof. It is a consequence of theorem 6.1 from [START_REF] Lasiecka | Unified theory for abstract parabolic boundary problems-a semigroup approach[END_REF] for problems with mixed (Robin) boundary conditions in the case α = 1/4 and k = 0.

Remark 36. The same type of results hold for Dirichlet boundary conditions with f ∈ W 3 2 ,1 L 2 (∂U ) if the compatibility condition c 0|∂U = f (0) is satisfied (see [START_REF] Lasiecka | Unified theory for abstract parabolic boundary problems-a semigroup approach[END_REF]).

Remark 37. 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.

Initial conditions Let

where c b,lum and c o,bot are the Dirichlet boundary conditions for respectively the butyrate and the oxygen (see Table 2).

  (a) Interaction kernels F a (with k a = 1, blue) and D a (with d a = 1, orange). The function F a describes repulsive forces, while D a renders density dependant interactions in jump rates. (b) Regulation function y → R(y, K, κ) for various values of K and κ. The parameter K tunes the location of the transition and κ its slope.

Figure 1 :

 1 Figure 1: Interactions and regulations. Shapes of interaction kernels and regulation functions used in the model

  )-(48), γ m,l (c o , c b ) = 0 if l = pc, and, if l = ent,

  , the concentration distributions z → c N o and z → c N b obtained with the PDMP for different scaling parameters N ∈ {15, 20, 50, 60} are displayed at t = 29 h together with the functions z → c

Figure 2 :

 2 Figure2: Total populations over time : The total population of progenitor cell and enterocytes are plotted over time for the PDMP (with population scaling = 60) and the SA equation. The averaged over 30 PDMP repetitions are plotted (continuous line) together with quantiles 2.5% and 97.5% (thin dots). The deterministic self-aggregation limit is plotted with thick dots.

Figure 4 :Figure 5 :

 45 Figure 4: Evolution of the residual R N (z, t) for increasing N at different t and z, grouped by z: we display for different times (8h, 16h, 24h, 29h) and different z (34 µm, upper left; 42 µm, upper right; 59 µm, lower left; 68 µm lower right) the evolution of the residual for growing N (N ∈ {15, 20, 25, 30, 35, 40, 45, 50, 55, 60}), grouped by z. To visualize all the graphs in a common plot, we normalized the residuals by the value of the residual for N = 15: we then plot the functions N → R N (z, t)/R 15 (z, t)

  (ν nν) * G + (z) = l∈T Z G + (z -y) (ν nν)(dy, dl) = l∈T Z G + (z -y) νn (dy, dl) -k∈T Z G + (z -u) ν(du, dk) = l∈T Z k∈T Z G + (z -y) -G + (z -u) γ(dy, dl, du, dk) ≤ 2δ||G + || ∞ + X \E ||G + || Lip δ γ(dy, dl, du, dk) ≤ δ(2||G + || ∞ + ||G + || Lip ).

  l∈T Z k∈T Z G + (z -y) -G + (z -u) γ(dy, dl, du, dk) ≤ 2δ||G + || ∞ + X \E ||G + || Lip δ γ(dy, dl, du, dk) ≤ δ(2||G + || ∞ + ||G + || Lip ).As this is true for any δ > ϵ, we get that(ν nν) * G + (z) ≤ ϵ(2||G + || ∞ + ||G + || Lip ). (69)The same kind of upper bound can be found for|(ν n -ν) * G -(z)|. And so, ∀z ∈ Z, |(ν n -ν) * G(z)| ≤ 2ϵ||ν|| T V (3||G|| ∞ + ||G|| Lip ).

  for the reaction-diffusion problem (73). 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 10 in the main text).Let U ⊂ R d be a bounded open set of regularity C ∞ . We denote ∂U the boundary of U and consider the homogeneous Robin boundary condition∂ n u |∂U + bu |∂U = 0 (70)where b is a function in C ∞ ( Ū ) and ∂ n is the normal derivative on ∂U . It is a standard result (see[START_REF] Pazy | Semigroups of linear operators and applications to partial differential equations[END_REF][START_REF] Amann | Nonhomogeneous linear and quasilinear elliptic and parabolic boundary value problems[END_REF]) that the unbounded operator A = σ∆ whose domain D(A) is the closure in H 2 (U ) of the set of function in C 2 ( Ū ) which satisfy (70) 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 :

  ∂U ) , then f ∈ L 2 ([0, T ], H 3 2 (∂U )) and ∂f ∂t ∈ L 2 ([0, T ], H -1 2 (∂U )), so f ∈ C([0, T ], H 1 2 (∂U )) and the injection is continuous.• Similarly, if c ∈ W 2,1 L 2 (U ) then c ∈ C([0, T ], H 1 (U )) and the injection is continuous. With these definitions and properties, we have the following resultProposition 35. If b ∈ C ∞ ( Ū ), c 0 ∈ H 1 (U ), f ∈ W 1 2 ,1 L 2 (∂U ) , and u ∈ W 0,0 L 2 (U ) , then c(., t) = S(t)c 0 + (λI -σ∆)( t 0 S(t -s)D(f )(., s)ds) + t 0 S(t -s)u(., s)ds is the unique weak solution in W 2,1 L 2 (U ) of the following evolution problem ∂ t c -σ∆c = u, (73)∂ n c |[0,T ]×∂U + bc |[0,T ]×∂U = f c(., 0) = c 0

  simulations of the PDMP, we sample randomly the position of 55 × N progenitor cells in ρ pc and 100 × N enterocytes in ρ ent .The initial conditions c b0 and c o0 are defined withc b0 (z) = z z max c b,lum c o0 (z) = z max -z z max c o,bot

  2 s * D(z), c 2 s * ψ a (z)) is Lipschitz uniformly for s ≤ t. Indeed, z → ξ 2 s * D(z) is Lipschitz uniformly for s ≤ t under Hyp. 14. It is the case of c 2 s * ψ a (z) as well because sup s≤t ||c 2 s || L 2 (Za) Nm < +∞ from Prop. 10 and ψ a (z) is Lipschitz on Ẑ under Hyp. 8. Eventually, ∀l ∈ T , q k

Evolution of the residual R

			t= 8h				t= 16h
		1.00			1.00
		0.98			0.98
	rescaled rel. err.	0.90 0.92 0.94 0.96			0.92 0.94 0.96
					0.90
		0.88
		0.86			0.88
					0.86				z=34
		1.00	t= 24h		1.025		t= 29h		z=42 z=59 z=68
		0.95			1.000
	rescaled rel. err.	0.85 0.90			0.900 0.925 0.950 0.975
		0.80			0.875
					0.850
		20	30 scaling parameter Ns 40	50	60	20	30 scaling parameter Ns 40	50	60
	Figure 3:

N (z, t)

for increasing N at different t and z, grouped by times: we

  display for different times (8h, upper left; 16h, upper right; 24h, lower left and 29h, lower right) and different z (34 µm,42 µm,59 µm,68 µm) the evolution of the residual for growing N (N ∈ {15, 20, 25, 30, 35, 40, 45, 50, 55, 60}), grouped by times. To visualize all the graphs in a common plot, we normalized the residuals by the value of the residual for N = 15: we then plot the functions N → R N (z, t)/R 15 (z, t).

				z=:34 m				z=:42 m
		1.025			1.025
		1.000			1.000
	rescaled rel. err.	0.925 0.950 0.975			0.900 0.925 0.950 0.975
		0.900			0.875
		0.875			0.850
				z=:59 m				z=:68 m	t=8h t=16h t=24h
			1.00			1.00			t=29h
	rescaled rel. err.	0.90 0.95			0.85 0.90 0.95
			0.85
						0.80
			0.80
			20	30 scaling parameter Ns 40	50	60	20	30 scaling parameter Ns 40	50	60

Table 2 :

 2 Parameters for Sect. 6.1. For references, see[START_REF] Darrigade | A PDMP model of the epithelial cell turn-over in the intestinal crypt including microbiota-derived regulations[END_REF]. max. rate of progenitor cell to enterocyteK pc,ent [z] = K pc,pc [z] µm space regulation of progenitor cell to diff. κ pc,ent [z] = κ pc,pc [z]µm space regulation of progenitor cell to diff. K pc,ent[but] = K pc,pc [but] µm butyrate regulation of progenitor cell to diff. κ pc,ent[but] = κ pc,pc [but] µm butyrate regulation of progenitor cell to diff.

		Value	Unit	Description
	z max	100	µm	height of the crypt
	r 0	12.5	µm	shape parameter of the crypt
	ϵ	= 0.1 × r 0	µm	shape parameter of the crypt
	a	10	µm	cell diameter
	k a	15	µm × h -1	maximal intensity of F
	d a	1	∅	maximal value of D
	q ∞ pc,pc	0.11	h -1	maximal rate of progenitor cell division
	K pc,pc [z]	20	µm	space regulation for progenitor cell div.
	κ pc,pc [z]	7.5	µm	space regulation for progenitor cell div.
	K pc,pc [dens]	20	cell	cell density regulation for progenitor cell div.
	κ pc,pc [dens]	6	cell	space regulation for progenitor cell div.
	K pc,pc [but]	2.5		butyrate regulation for progenitor cell div.
	κ pc,pc [but]	5		butyrate regulation for progenitor cell div.
	q ∞ pc,ent	0.125	h -1
	q ∞ ent,∅	0.17	h -1	maximal rate of diff. extrusion
	K ent,∅ [z]	90	µm	space regulation for diff. extrusion
	κ ent,∅ [z]	10	µm	space regulation for diff. extrusion
	K ent,∅ [dens]	20	cell	cell density regulation for diff. extrusion
	κ ent,∅ [dens] = κ pc,pc [dens]	cell	cell density regulation for diff. extrusion
	s o	-4	∅	stoechiometric coefficient for O 2
	s b	-1	∅	stoechiometric coefficient for butyrate
	σ o	3.6 × 10 4	µm 2 /h	diffusion coefficient for O 2
	σ but	= σ o	µm 2 /h	diffusion coefficient for butyrate
	γ ∞	1.6 × 10 2	10 -15 × mmol h×µm 2	max. reaction speed of β-oxydation
	K	1840 1/5	mM	affinity of β-oxydation
	c b,lum	5	mM	luminal concentration of butyrate
	c o,bot	10	a.u.

O 2 concentration at the bottom of the crypt

Proof. As ν is finite L ∞ (X , |ν|) ⊂ L 1 (X , |ν|), and continuous functions are dense in L 1 (X , |ν|).

Proposition 32. Let α ∈ R + , ν ∈ M F + (X ),

Proof. Suppose that 0 ≤ α ≤ 1. Let ϵ ≥ 0, F a closed set of X . Then αν(F ) ≤ ν(F ϵ ) + ϵ as ν is a positive measure. Moreover

Using what we have shown above:

This completes the proof.

Proposition 33. Let α > 0 and ν, µ ∈ M F + (X ).

Proof. Let F be a closed set of X .

The second inequality implies that

From the definition of the Prokhorov metric, this concludes the proof.

Proof. We remind the reader that if

We start by treating the case when ν is the null measure. Then, ∀z ∈ Z,

Let us suppose now that ||ν|| T V > 0. As ||ν n || T V → ||ν|| T V , we also suppose that ||ν n || T V > 0 ∀n ∈ N. We want to use proposition 27 p.29 which is valid for probability measures. We have that

Indeed:

C Technical results for the proof of theorem 21

Lemma 38. Under hypotheses 2, let (c s , ξ s ) s≥0 be a solution of (25-26) on [0, T ], let t ∈ [0, T ] and u * (z) ∈ C 1 (Z, R). The transport equation on [0, T ]

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

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

Proof. This is an application of theorem 6.1 in [START_REF] Perthame | Transport Equations in Biology[END_REF], which relies on the continuous derivability of the flow associated to the ODE (76). To prove the latter, one must show that the vector field (t, z) → F(z, ξ t * F (z)) is continuously derivable on [0, T ]×Z, which is a standard consequence of the C 1 regularity assumption (1) in hypothesis 2. The detailed proof can be found in A.3 p.119 of [START_REF] Darrigade | Modélisation du dialogue hôte-microbiote au voisinage de l'épithélium de l'intestin distal[END_REF].

Let (c 1 t , ξ 1 t ) and (c 2 t , ξ 2 t ) be two solutions of (25-26) on [0, T ] with ξ i t positive, and let Z 1 (s; t, z) and Z 2 (s; t, z) be the associated characteristics defined in (76). Then, the following proposition is true.

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

Proof. Proof can be found in A.4 p.125 of [START_REF] Darrigade | Modélisation du dialogue hôte-microbiote au voisinage de l'épithélium de l'intestin distal[END_REF].

D Exact expression and parameters values for the numerical implementation

Shape of the domain The function ϕ appearing in equation ( 45) is defined as

where we assume that r 0 is small compared to z max , and ϵ is small compared to r 0 .