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

Marie Haghebaert
email: marie.haghebaert@inrae.fr

Claire Cherbuy
email: claire.cherbuy@inrae.fr

Simon Labarthe
email: simon.labarthe@inrae.fr

Béatrice Laroche

Beatrice Laroche
email: beatrice.laroche@inrae.fr

B Laroche

A PDMP model of the epithelial cell turn-over in the intestinal crypt including microbiota-derived regulations

Keywords: Piecewise Deterministic Markov Processes, Porous Media Equation, Crypt model, Host-microbiota interactions

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

1.1 Colonic crypts: structure, regulations and interactions with the microbiota The mammalian colon harbours a complex microbial community called the colonic microbiota. This microbiota is in constant and reciprocal interactions with its host and plays a key role in its nutrition and immunity. At the forefront of these interactions is the intestinal epithelium. It is a monolayer of cells lining the inside of the gastrointestinal wall, wrapped with a mucus layer in direct contact with the gut content. The intestinal epithelium is very densely folded and its smallest fold is called the intestinal crypt. The intestinal crypt is also the simplest unit of the host-microbiota crosstalk. It is an intensively studied biological system from various standpoints: stem cell biology [START_REF] Barker | Adult intestinal stem cells: Critical drivers of epithelial homeostasis and regeneration[END_REF][START_REF] Sasaki | Reg4+ deep crypt secretory cells function as epithelial niche for Lgr5+ stem cells in colon[END_REF], cell differentiation [START_REF] Tetteh | Replacement of Lost Lgr5-Positive Stem Cells through Plasticity of Their Enterocyte-Lineage Daughters[END_REF][START_REF] Tóth | Early commitment and robust differentiation in colonic crypts[END_REF], biophysic [START_REF] Shyer | Villification: How the gut gets its villi[END_REF][START_REF] Serra | Self-organization and symmetry breaking in intestinal organoid development[END_REF][START_REF] Krndija | Active cell migration is critical for steady-state epithelial turnover in the gut[END_REF] or interactions with the microbiota [START_REF] Tomas | Early colonizing Escherichia coli elicits remodeling of rat colonic epithelium shifting toward a new homeostatic state[END_REF][START_REF] Naito | Lipopolysaccharide from crypt-specific core microbiota modulates the colonic epithelial proliferation-to-differentiation balance[END_REF].

In the colon, a crypt is roughly shaped as a test tube of size 200 × 50 µm and contains hundreds of cells [START_REF] Neumann | Gut commensal bacteria and regional Wnt gene expression in the proximal versus distal colon[END_REF][START_REF] Sunter | A comparison of cell proliferation at different sites within the large bowel of the mouse[END_REF]. The colonic crypt is organised in three zones along its main axis [START_REF] Barker | Adult intestinal stem cells: Critical drivers of epithelial homeostasis and regeneration[END_REF], as shown in the Fig. 1. The bottom of the crypt, the farthest zone from the microbiota, is the stem cell niche. In this specific microenvironment reside epithelial stem cells and deep crypt secretory (DCS) cells [START_REF] Sasaki | Reg4+ deep crypt secretory cells function as epithelial niche for Lgr5+ stem cells in colon[END_REF]. Stem cells produce all other types of epithelial cells through successive divisions and differentiation, while DCS cells are differentiated cells with a key role in maintaining a favourable environment for stem cells. Above the stem cell niche are found progenitor cells. Progenitors are still dividing but they have started their differentiation path. Finally, in the upper part of the crypt, transit-amplifying cells have terminally differentiated in mostly two types of cells: enterocytes and goblet cells. Enterocytes are involved in the uptake of nutrients and water from the gut lumen and goblet cells secrete a protective layer of mucus covering the epithelium.

A striking feature of the colonic crypt is cell migration. Epithelial cells are produced at the bottom of the crypt by division of stem cells and are removed from the epithelium by extrusion at the top of the crypt. Cells therefore migrate from the bottom to the top of the crypt and differentiate along the way.

The functioning of the crypt relies on numerous and complex regulations. Of interest are regulations relying on contacts between neighbouring cells and regulations by molecules whose concentrations change along the crypt main axis. Other noteworthy examples of such regulation in the crypt are the Notch pathway [START_REF] Sancho | Stem cell and progenitor fate in the mammalian intestine: Notch and lateral inhibition in homeostasis and disease[END_REF][START_REF] Sasaki | Reg4+ deep crypt secretory cells function as epithelial niche for Lgr5+ stem cells in colon[END_REF] and the EPH/ephrin signaling [START_REF] Batlle | β -catenin and TCF mediate cell positioning in the intestinal epithelium by controlling the expression of EphB/EphrinB[END_REF]. For several transcription factors (Wnt, BMP or Hedgehog), the concentration is a gradient between the bottom and the top of the crypt [START_REF] Wang | Bioengineered Systems and Designer Matrices That Recapitulate the Intestinal Stem Cell Niche[END_REF]. These transcription factors play an important part in crypt organization. For example, Wnt is secreted by pericryptic cells at the bottom of the crypt and is necessary to maintain the pool of stem cells [START_REF] Degirmenci | GLI1-expressing mesenchymal cells form the essential Wnt-secreting niche for colon stem cells[END_REF].

Epithelial cells interact with the colonic microbiota in multiple ways. A major interaction is the regulation of colon hypoxia, i.e. the depletion of oxygen in gut lumen, regulated by microbe-derived butyrate. Butyrate is a short-chain fatty acid produced at homeostasis by bacteria of the colonic microbiota by anaerobic fermentation of food fibers [START_REF] Litvak | Colonocyte metabolism shapes the gut microbiota[END_REF]. Dioxygen diffuses to the crypt from nearby blood vessels. At the bottom of the crypt, dividing host cells relie on anaerobic glycolysis (oxygen-free conversion of glucose to lactate) for their energy production, and consume little oxygen [START_REF] Fan | A bioassay to measure energy metabolism in mouse colonic crypts, organoids, and sorted stem cells[END_REF][START_REF] Lopez | Virulence factors enhance ¡i¿citrobacter rodentium¡/i¿ expansion through aerobic respiration[END_REF]. At the top of the crypt, differentiated epithelial cells consume butyrate and dioxygen through β -oxydation coupled to respiration to produce energy [START_REF] Litvak | Colonocyte metabolism shapes the gut microbiota[END_REF][START_REF] Kelly | Crosstalk between microbiota-derived short-chain fatty acids and intestinal epithelial HIF augments tissue barrier function[END_REF]. This metabolic activity of differentiated epithelial cells in the upper part of the crypt renders the colonic epithelium hypoxic [START_REF] Litvak | Colonocyte metabolism shapes the gut microbiota[END_REF]. Little oxygen diffuses from the gastrointestinal wall to the colonic lumen which helps to maintain anaerobiosis in the colonic lumen. This is very important as some beneficial microbial species are obligate anaerobes whereas pathogens (Salmonella enterica, Citrobacter rodentium) or dysbiotic bacteria of the phylum Proteobacteria can utilize oxygen for their growth [START_REF] Litvak | Colonocyte metabolism shapes the gut microbiota[END_REF][START_REF] Byndloss | The germ-organ theory of non-communicable diseases[END_REF]. Conversely, if microbe-derived butyrate drops off, differentiated epithelial cells switch their metabolism to anaerobic glycolysis, leading to oxygen release in the lumen which initiates a vicious circle harmful to butyrate producers and beneficial to dysbiotic bacteria. The idea is that the host keeps the colonic microbiota "on a leash" by controlling the availability of oxygen (and other electron acceptors) in the colonic lumen.

Butyrate also affects cellular processes beyond metabolism [START_REF] Martin-Gallausiaux | Scfa: mechanisms and functional importance in the gut[END_REF]. For example, butyrate has been reported to change the activity of important regulators of epithelial cell activity such as the nuclear receptor PPAR-γ [START_REF] Alex | Short-Chain Fatty Acids Stimulate Angiopoietin-Like 4 Synthesis in Human Colon Adenocarcinoma Cells by Activating Peroxisome Proliferator-Activated Receptor[END_REF], the heat shock protein Hsp25 [START_REF] Uchiyama | Butyrate and bioactive proteolytic form of Wnt-5a regulate colonic epithelial proliferation and spatial development[END_REF] or the cell-cycle regulator Foxo3 [START_REF] Kaiko | The Colonic Crypt Protects Stem Cells from Microbiota-Derived Metabolites[END_REF]. Due to these interactions, butyrate was shown to reduce the proliferation of stem and progenitor cells, even inducing their apoptosis at high concentration [START_REF] Kaiko | The Colonic Crypt Protects Stem Cells from Microbiota-Derived Metabolites[END_REF][START_REF] Uchiyama | Butyrate and bioactive proteolytic form of Wnt-5a regulate colonic epithelial proliferation and spatial development[END_REF], and to impact the differentiation processes of epithelial cells [START_REF] Pearce | Intestinal enteroids recapitulate the effects of short-chain fatty acids on the intestinal epithelium[END_REF]. The effect of butyrate on stem and progenitor cells activity is observed in vitro at luminal concentration of butyrate. Therefore, it is hypothesized that the consumption of butyrate by differentiated cells of the upper part of the crypt is sufficient to drop concentrations of butyrate at a low level in the stem cell niche [START_REF] Kaiko | The Colonic Crypt Protects Stem Cells from Microbiota-Derived Metabolites[END_REF]. However, the concentrations of butyrate or oxygen are hard to measure in vivo and are therefore most of the time missing of the picture. Our model is a step to fill this gap.

Existing models and achievement

There is a rich literature of intestinal crypt models, see [START_REF] Fletcher | Multiscale modelling of intestinal crypt organization and carcinogenesis[END_REF] and [START_REF] Almet | Modeling perspectives on the intestinal crypt, a canonical system for growth, mechanics, and remodeling[END_REF] for a very recent review and references therein. Most of the proposed models are stochastic individualbased models (IBM) at the scale of epithelial cells. There exists also macroscopic models formulated in terms of cells density [START_REF] Hannezo | Instabilities of Monolayered Epithelia: Shape and Structure of Villi and Crypts[END_REF][START_REF] Hannezo | Interplay of migratory and division forces as a generic mechanism for stem cell patterns[END_REF], notably coming from physics. The link between the microscopic scale of IBMs and the macroscopic scale of cell density models has rarely been explored [START_REF] Murray | Comparing a discrete and continuum model of the intestinal crypt[END_REF], most probably due to the complexity of IBMs and the lack of mathematical formalism. These crypt models are all largely motivated by issues arising from oncology such as clonal dynamics and stem cells competition in the crypt [START_REF] Lopez-Garcia | Intestinal stem cell replacement follows a pattern of neutral drift[END_REF][START_REF] Snippert | Intestinal crypt homeostasis results from neutral competition between symmetrically dividing Lgr5 stem cells[END_REF][START_REF] Thalheim | Linking stem cell function and growth pattern of intestinal organoids[END_REF] or morphogenesis [START_REF] Buske | On the biomechanics of stem cell niche formation in the gut -Modelling growing organoids[END_REF][START_REF] Langlands | Paneth Cell-Rich Regions Separated by a Cluster of Lgr5+ Cells Initiate Crypt Fission in the Intestinal Stem Cell Niche[END_REF][START_REF] Thalheim | Linking stem cell function and growth pattern of intestinal organoids[END_REF]. Some other models describe microbiota dynamics and digestive activity at the colonic scale, with however a very poor description of the epithelium [START_REF] Muñoz-Tamayo | Mathematical modelling of carbohydrate degradation by human colonic microbiota[END_REF][START_REF] Labarthe | A mathematical model to investigate the key drivers of the biogeography of the colon microbiota[END_REF][START_REF] Cremer | Effect of water flow and chemical environment on microbiota growth and composition in the human colon[END_REF][START_REF] Moorthy | A spatially continuous model of carbohydrate digestion and transport processes in the colon[END_REF].

In this work, our objective is to account for interactions with the microbiota in a healthy crypt model. To our knowledge, this is a new step in the modelling of intestinal crypts. In comparison with previous literature, we simplify crypt geometry and physical interactions between cells while retaining their essential features. Indeed these elements are not as important for understanding interactions with the microbiota as they are for other issues explored by crypt models, such as cell proliferation during oncogenesis or cryptogenesis. Instead, we put the emphasis on regulation processes directly or indirectly involved in crypt-microbiota cross-talk (see Sect. 2).

We propose an original model formulated as a piecewise deterministic Markov process (PDMP), which provides a rigorous mathematical framework. Our model couples an atomic measure that represents epithelial cells with spatially-structured concentrations of chemicals that mediate interactions between the crypt and the microbiota. The epithelial cells experience stochastic events governing their fate, and between two such events they move according to ODEs describing the mechanical constraints exerted on them, while reaction-diffusion PDEs model the chemical concentrations dynamics. Such a mixed ODE-PDE framework for a PDMP was already proposed e. g. in [START_REF] Buckwar | An exact stochastic hybrid model of excitable membranes including spatio-temporal evolution[END_REF] in other biological contexts and our work follows the same lines. This detailed mathematical setting allows in particular for the derivation of a deterministic limit of the PDMP, which is more suitable than the stochastic version for intensive simulation required for instance in sensitivity analysis. On the other hand, as outlined in [START_REF] Almet | Modeling perspectives on the intestinal crypt, a canonical system for growth, mechanics, and remodeling[END_REF], the PDMP model provides outputs directly comparable to biological data, and which are more difficult to access with the deterministic model such as migration speed.

The description of the model is done in Sect. 2, and its deterministic limit is formally derived in Sect. 3. In Sect. 4, we propose a stepwise validation of the IBM by decomposing it in elementary mechanisms that are studied sequentially and assessed by comparison with experimental values. Next, the global model is used to study the microbe-host interactions through butyrate and oxygen use. In Sect. 5, the deterministic limit of the PDMP model is implemented and explored by sensitivity analysis (Morris method). These models are finally used to explore a therapeutic strategy to accelerate recovery after infection by an enteric pathogen in Sect. 6.

A PDMP individual based-model of the crypt

We constructed an individual-based model of a colonic crypt at the scale of epithelial cells, interacting with butyrate or dioxygen which diffuse in cells environment. As shown in Fig. 1, this model accounts for 1) cell migration due to mechanical interactions between neighbouring cells, 2) cell division, extrusion and differentiation regulated by chemicals concentrations, cell density and spatial position, and 3) diffusion and consumption of chemicals, leading to chemical gradients in the crypt. This model is a PDMP, or hybrid system: it mixes discrete stochastic events modeling sudden changes in the cell population such as differentiation, division or extrusion, and a deterministic flow modeling continuous population dynamics between this events, such as migration. With these elementary mechanisms, we want to reproduce biomarkers of an healthy crypt, such as proportion and spatial distribution of cell populations, migration speed of epithelial cells or crypt turn-over duration.

Model state description

We assume that the crypt has cylinder symmetry around its vertical axis [0, z max ] and is characterized by a truncated test-tube like shape (see Fig. 1b), with radius z → r(z) defined as

f (z) =            r 0 z+ε r 0 2 -z+ε r 0 if z ≤ r 0 -ε r 0 if r 0 -ε < z < z max -r 0 + ε r 0 2 -z max -z+ε r 0 2 -z max -z+ε r 0 if z ≥ z max -r 0 + ε (1)
where we assume that r 0 is small compared to z max , and ε is small compared to r 0 . For modeling interactions between epithelial cells and microbiota activity, the most relevant component of cell position in the crypt is the position along the crypt vertical axis. Therefore, we consider a spatial model of dimension one where each cell is characterized by its position z ∈ [0, z max ] only, that is the projection of its actual spatial location on the crypt vertical axis. Consequently, several cells can be described by the same coordinate z, as several cells can be at the same height z in the 3D crypt.

A cell also has a type ℓ ∈ T , where T is a finite set comprising 5 different cell types: stem cell (sc), DCS cells (dcs), progenitor cells (pc), goblet cells (gc), enterocytes (ent).

Therefore, setting X = [0, z max ] × T , a cell is characterized by the vector be the set of positive finite atomic measures on X , with ∑ 0 i=1 δ x i (dx) the null measure by convention. A population of n epithelial cells at time t is represented by:

x = (z, ℓ) ∈ X . Let M P + (X ) = n ∑ i=1 δ x i (dx) : n ∈ N * ; ∀i, x i ∈ X (
ν t (dx) = n ∑ i=1 δ x i (t) (dx) ∈ M P + (X )
where x i (t) = (z i (t), ℓ i (t)) is the vector of position and type of cell i at time t. For each cell type ℓ ∈ T , we denote by ν ℓ (dz) ∈ M P + ([0, z max ]) the restriction of ν to a specific cell type (that is to say ν ℓ (Z) = ν(Z × {ℓ}) for any real borel set Z).

For a real-valued measurable function g on X , the integral of g against ν will be

⟨ν, g⟩ = X g(x)ν(dx) = n ∑ i=1 g(x i ).
For a real-valued measurable function h on R, we define the convolutions

ν * h(z) = ∑ ℓ∈T [0,z max ] h(z -y)ν(dy, ℓ) = n ∑ i=1 h(z -z i )
and, for ℓ ∈ T ,

ν ℓ * h(z) = [0,z max ] h(z -y)ν ℓ (dy) = n ∑ i=1 h(z -z i )1 {ℓ i =ℓ} .
Butyrate and dioxygen will be described by their respective concentrations c b (z,t) and c o (z,t) depending on space and time in suitable function spaces (see Sect. 2.2). In the following sections, we describe in details the behaviour of cells, as well as butyrate and dioxygen concentrations, by providing in Sect. 2.2 the deterministic flow of the PDMP and the stochastic events in Sect. 2.3.

Deterministic mechanisms of the PDMP model

The deterministic part of the PDMP models the cell movements and the reactiondiffusion of chemicals Mechanical interactions for cell migration. Cells are submitted to forces exerted by neighbouring cells and the epithelial basal membrane. We assume that for a given cell the components that are normal to the basal membrane are compensated by reaction and adhesion forces so that only the tangential forces are active (see Fig. 1c, right). These tangent components are modelled as repulsive forces. Moreover, as our model is 1D, we provide a phenomenological expression for the projection on the vertical axis of the total force exerted on cell i in a population of n cells as:

F i (t) = φ (z i (t)) n ∑ j=1 F a (z i (t) -z j (t)) = φ (z i (t))ν t * F a (z i (t)) (2)
where φ and F a are fixed functions that we describe now. In the 3D crypt, only neighbouring cells interact. This translates in our model into F a (z) := F( z a ) where F is odd, compactly supported in [-1, 1], and a is the typical diameter of one cell. Moreover, in the original 3D conformation, cells that are sufficiently close to each other exert a force mostly normal to the z-axis. This translates in our 1D model into a decreasing amplitude as z goes to zero and F a (0) = 0. Therefore, we used the following expression for the force

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. (3)
where k a is a scaling parameter driving the maximal intensity of the force. F a is plotted in Fig. 2a with k a = 1.

The function φ appearing in Eq.( 2) aims to take into account the impact of curvature at the top and bottom of the crypt: in the curved zones, the tangent interactions forces are exerted more radially (see Fig. 1c, left), reducing their impact on the longitudinal cell kinematics. Consequently a phenomenological correction is applied to the total force, which takes the form of a linear function of the crypt radius f (z), defined in (1), from 1 when f (z) = r 0 (i.e. no curvature) to 0 when z = z max or z = 0 (maximal curvature at both ends of the crypt). Namely, the correction is displayed in Fig. 1b and reads

φ (z) :=              f (z) r 0 - f (0) r 0 1- f (0) r 0 if z ≤ r 0 -ε 1 if r 0 -ε < z < z max -r 0 + ε f (zmax) r 0 - f (z) r 0 f (zmax) r 0 -1 if z ≥ z max -r 0 + ε (4)
Except for DCS cells (see below), the position of cell i evolves according to laws of classical mechanics taken in a non-inertial regime. This regime is known to correctly render the motion of individual cells in cell tissues where inertial terms can be neglected against other mechanisms [START_REF] Purcell | Life at low Reynolds number[END_REF]. The resulting equation of cell motion is

dz i dt (t) = F i (t) if ℓ i ̸ = dcs (5)
with F i (t) defined in [START_REF] Almet | Modeling perspectives on the intestinal crypt, a canonical system for growth, mechanics, and remodeling[END_REF]. Note that, as φ (0) = φ (z max ) = 0, the speed of cells at the boundaries z = 0 and z max is zero. DCS cells are assumed to be "inactive", that is to say motionless ( dz i dt (t) = 0), not dividing, not differentiating, not getting extruded. This is a simplification since it has been evidenced that the behaviour of DCS cells is quite complex and different from the rest of epithelial cells [START_REF] Sasaki | Reg4+ deep crypt secretory cells function as epithelial niche for Lgr5+ stem cells in colon[END_REF]. For example, they are thought to migrate downwards the crypt. Although inactive in the model, they exert mechanical forces on neighbouring cells and contribute to local cell density near the bottom of the crypt.

Reaction-diffusion of chemicals. As explained in section 1.1, butyrate is produced by the microbiota from fibers fermentation and O 2 is transported to the intestinal epithelium by the vascular system. Both molecules diffuse in the crypt and are consumed by goblet cells and enterocytes for energy production through β -oxidation. We summarize this metabolic pathway with the equation

1 butyrate + 4 O 2 -→ energy + byproducts. Concentration of chemical i ∈ {butyrate, O 2 } evolves according to the reaction/diffusion equation on ]0, z max [ a := ]-a/2, z max + a/2[ ∂ t c i -σ i ∂ zz c i = -s i γ ∞ β c 4 o c b c 4 o c b + K 5 β (ν ent t + ν gc t ) * ψ a (z) (6)
where σ i > 0 and s i are respectively the diffusion constant and the stoechiometric coefficient of metabolite i. γ ∞ β is the maximal reaction speed of β -oxydation, K β is a positive constant, and (ν ent t + ν gc t ) * ψ a (z) is the quantity of cell at z catalyzing the reaction. Indeed, we introduce the function

D a D a (z) := d a D( z a ), where D(z) := exp -z 2 1-z 2 if 0 ≤ |z| ≤ 1 0 otherwise, (7)
and d a is a scaling parameter (see Fig. 2a with

d a = 1). Note that D a is continuous, supported on ] -1, 1[ and has sup z D a (z) = D(0) = 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 . (8)
Then, for all z ∈ [0, z max ], ]0,z max [ a ψ a (zy)dy = 1. We see with this kernel definition that the extended domain ]0, z max [ 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 ( 9
)
and for oxygen:

∂ z c o = 0 at z = z max + a 2 c o = c o,bot at z = -a 2 ( 10
)
For butyrate, we set a no flux boundary condition at the bottom of the crypt as we consider that butyrate cannot diffuse across the epithelium (it can only be actively absorbed by epithelial cells). For oxygen, we set a no flux boundary condition at the top of the crypt as we assume that the concentration of oxygen at the top of the crypt and in the lumen zone near the crypt are very close. 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, and will be chosen depending on the biological context we want to model. Furthermore initial conditions 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. Deterministic flow. With classical arguments of system dynamics and PDE analysis (see Annex B), we can check for global existence of solution to the system (5-6), and define the associated deterministic flow for all t > 0.

A t : H 1 (]0, z max [ a ) 2 × M P + (X ) -→ H 1 (]0, z max [ a ) 2 × M P + (X ) (c, ν) -→ A t (c, ν)
where

A t (c, ν) = A 0 t (c, ν), Ãt (c, ν) with Ãt (c, ν) := ⟨ν,1⟩ ∑ i=1 δ (A i t (c,ν),ℓ i ) .
A t (c, ν) is the vector of solutions of the system (5-6) at time t for initial conditions (c, ν), with A 0 t (c, ν) the solutions (c b , c o ) of Eq. ( 6) and A i t (c, ν), for 1 ≤ i ≤ ⟨ν, 1⟩, the solution of Eq. [START_REF] Bansaye | Stochastic models for structured populations: scaling limits and long time behavior[END_REF]. Note that ⟨ν, 1⟩ is the integral of the atomic measure ν against the function 1 and gives the number n of cells in the population ν. This number is constant under the deterministic flow.

By a direct application of the fundamental theorem of calculus, we have the following equalities, holding for f

∈ C 1 (R + × [0, z max ] , R), s,t ∈ R + and 1 ≤ i ≤ n, f t (A i t (c, ν)) -f s (A i s (c, ν)) (11)
= t s φ (A i u (c, ν)) Ãu (c, ν) * F a (A i u (c, ν))∇ z f u (A i u (c, ν)) + ∂ u f u (A i u (ν)) du and ⟨ Ãt (c, ν), f t ⟩ -⟨ Ãs (c, ν), f s ⟩ = t s ⟨ Ãu (c, ν), φ (•) Ãu (c, ν) * F a (•)∇ z f u (•) + ∂ u f u (•)⟩ du. ( 12
)

Stochastic modeling of cell fate

The PDMP takes into account seven possible event types in cell population dynamics, modeled as stochastic punctual events, or jumps. They are indexed by k ∈ E the finite set of event types. These events are the symmetric division of stem and progenitor cells, the death by extrusion of enterocytes and goblet cells and the differentiation of stem cells into progenitors, as well as the differentiation of progenitors into enterocytes and goblet cells (see Fig. 1a for an overview of the events). Their construction is described below.

Stochastic jumps during cell fate. Jump functions are defined for all type of cells. When at T the cell i undergoes a jump of type k, the population is modified as

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

µ k : X → M P (X ), (13)
with M P (X ) = {∑ n i=1 ±δ x i (dx) : n ∈ N * ; ∀i, x i ∈ X } the set of atomic finite measures on X . More precisely,

• if k is a division, a new daughter cell is introduced at a distance λ (z) above (resp. under) the mother cell in the lower (resp. upper) part of the crypt, leading to

µ k (z, ℓ) = δ (z+λ (z),ℓ) with λ (z) =      a 2 if z ≤ z max /2 -a, -a 2 if z ≥ z max /2 + a, a 2 -a z-(z max /2-a) 2a otherwise, so that ∀z ∈ [0, z max ], z + λ (z) ∈ [0, z max ], • if k is an extrusion µ k (x) = -δ x , • if k is a differentiation of a cell of type ℓ into a type ℓ k then µ k (z, ℓ) = -δ (z,ℓ) + δ (z,ℓ k ) .
Jump rates. We model three biological regulations acting on events of cell fate. The first one depends on cell variables (position and cell type). The dependence on the position models the action of morphogenes, in particular Wnt, on the cell differentiation: we assume a fixed linearly decreasing gradient of Wnt in the crypt, so that the Wnt density can be modeled by the position z ∈]0, z max [ up to scaling parameters. The dependence on the cell type ensures that only the seven possible jump types occur, by setting the jump rate to zero when necessary.

The second one depends on the local cell density d defined by

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

d(z) = a -a D a (u) du min(a,z max -z) max(-a,-z) D a (u) du . ( 14
)
The function d(z) is then a normalization coefficient insuring that the kernel D a keeps the same integral along [0, z max ], included near the boundaries where the integral can be cut off, in order to account for the 1D approximation.

The third regulation depends on the local concentration of butyrate c b defined by

c b (z) = c b * ψ a (z)
where ψ a (z) is introduced in Eq. ( 8).

Then, the individual jump rate associated to each jump k ∈ E is defined as

q ∞ k q k (x, d(ν t , z), c b (z)) (15)
where q ∞ k is a positive constant equals to the maximal rate of jump and 0 ≤ q k (•) ≤ 1. To define the normalized jump rate q k , we introduce a generic regulation function R(y, K, κ) for y ≥ 0, K > 0, κ > 0, that will model the dependence of q k on the variable y (y = z for space regulation i.e. Wnt regulation, y = d(ν t , z) for density regulation or y = c b (z) for butyrate regulation) parametrized by jump and regulation dependant parameters K k [ j] and κ k [ j] where index k refers to the jump type and j to the regulation pathway (see equation [START_REF] Campolongo | An effective screening design for sensitivity analysis of large models[END_REF] for an explicit example). Namely, we define the piecewise polynomial function R(y, K, κ) (plotted on Fig. 2b for several values of K and κ) 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. (16
) Note that R connects continuously the values 0 and 1 for Kκ < z < K + κ and has value 1/2 in y = K. The normalized individual jump rate q k (x, d(ν t , z), c b (z)) associated to event k is a product over the regulations j of functions R( an activation or

(1 -R(•, K k [ j], κ k [ j]
)) for an inhibition. Below, we give the transition rate for division of stem cells:

q div,sc (z, ℓ, d(ν t , z), c b (z)) = 0 if not stem cell 1 {ℓ=sc} × regulation by Wnt (1 -R(z, K div,sc [z], κ div,sc [z])) × regulation by density (1 -R(d(ν t , z), K div,sc [dens], κ div,sc [dens]) × (1 -R(c b (z), K div,sc [but], κ div,sc [but]))
regulation by butyrate .

The form of the transition rate reflect different biological regulations based on the biological knowledge presented in Sect. 1.1. For other events and cell types, the form of the transition rate is explained in Table 1 and given explicitely in Appendix A.2. The parameter values q ∞ k , K k [ j] and κ k [ j] can be found in Table 11. Note that each jump type is not necessarily regulated by all three regulation pathways (for example, the rate of extrusion does not depend on butyrate concentration). When a regulation pathway is neutral for a jump, the corresponding term is set to the constant value 1. Remember as well that, as we supposed a stable linearly decreasing spatial gradient of Wnt, an inhibitory effect of Wnt corresponds in our model to an activation by the position z. Finally, the choice of using a product of regulation functions R(•) (and 1 -R(•)) instead of a sum corresponds to the biological hypothesis that the regulation pathways interact: if one of the regulation is close to zero, the overall jump rate is close to zero despite the contribution of others regulation terms.

Numerical implementation of the PDMP

We implemented the model in Python 3. The code sources can be found at https: //forgemia.inra.fr/decrypt/pdmp-model. To solve the PDMP, we use a standard Gillespie algorithm with an acceptation-reject method for stochastic events. The deterministic flow is solved with an explicit Euler method for the cell movement, and with a semi-implicit Euler scheme with operator splitting to compute the reactiondiffusion equations on butyrate and oxygen. Cells are sorted in a NumPy array according to their position, allowing to compute interactions only between groups of neighbouring cells. Numba package is used for further code acceleration.

3 Deterministic approximation of the PDMP.

To derive the deterministic approximation of the PDMP, we first set a functional framework in which the PDMP model is well-posed. Then, a formal approximation is derived in this setting.

Functional setting and well-posedness

With the definitions above, all the ingredients needed for the definition of a PDMP are now gathered (see e.g. [START_REF] Jacobsen | Point process theory and applications: marked point and piecewise deterministic processes[END_REF]):

the process state space G = H 1 (]0, z max [ a ) 2 × M P + (X ), which is a Polish space equipped with a distance defined as the product of the distance induced by the H 1 norm and the Prokhorov metric (see Proposition A.2.5.III [START_REF] Daley | An Introduction to the Theory of Point Processes[END_REF] and Proposition A.5 [START_REF] Darrigade | Modélisation du dialogue hôte-microbiote au voisinage de l'épithélium de l'intestin distal[END_REF]). a measurable and time-continuous deterministic flow A, also continuous with respect to its parameters. To prove the later, the only technical point is to show the continuity of the concentrations A 0 t (c, ν) with respect to ν and it can be deduce from Property 2.13, Chapter 4 of [START_REF] Darrigade | Modélisation du dialogue hôte-microbiote au voisinage de l'épithélium de l'intestin distal[END_REF].

a total jump intensity q : G → R + , defined by

q(c, ν) = ∑ k∈E ⟨ν,1⟩ ∑ i=1 q ∞ k q k (x i , d(ν t , z i ), c b (z i )) = ⟨ ∑ k∈E q ∞ k q k (., d(ν t , .), c b (.)), ν⟩ (18)
which only depends on the current process state, and is continuous in its arguments. Moreover, q is locally right-integrable in time along deterministic trajectories. Indeed, let (c, ν) ∈ G and h > 0, then

t+h t q(A s-t (c, ν))ds ≤ h|E | max k (q ∞ k )⟨ν, 1⟩ < +∞
since the cell population size only changes during stochastic jumps and does not vary during deterministic trajectories.

a Markov kernel r defined on G and taking values in the probability measures on G with

r(c, ν) = ∑ k∈E ⟨ν,1⟩ ∑ i=1 q ∞ k q k (x i , d(ν, z i ), c b (z i )) q(c, ν) δ (c,ν+µ k (x i )) ,
where q is the total jump intensity defined in [START_REF] Cherbuy | Effect of germfree state on the capacities of isolated rat colonocytes to metabolize n-Butyrate, glucose, and glutamine[END_REF]. This Markov kernel only depends on the current state (c, ν).

We have shown in Sect. 2.2 and Appendix B that the deterministic trajectory of the process between jumps is defined globally. The well-posedness of the PDMP will follow once we prove the stability of the PDMP ( [START_REF] Jacobsen | Point process theory and applications: marked point and piecewise deterministic processes[END_REF], theorem 3.2.1 and 7.3.1), that is to say that an almost surely finite number of jumps occurs during any finite time interval. This is done in the Appendix C and in [START_REF] Darrigade | Modélisation du dialogue hôte-microbiote au voisinage de l'épithélium de l'intestin distal[END_REF].

Deterministic limit of the PDMP

When the size of the cell population becomes large (1/C ≤ ⟨ν N , 1⟩/N ≤ C with N → ∞ for some C > 1), the corresponding sequence of models with a rescaled cell population (c N , ν N /N) and rescaled parameters

F N a = F a N , D N a = D a N and ψ N a = ψ a N converges in the space C([0, T ] , H 1 (]0, z max [ a )) 2 × D([0, T ] , M F + (X )),
where

D([0, T ] , M F + (X ))
is the space of càdlag functions from [0, T ] in the space M F + (X ) of finite positive measures on X , equipped with the Prokhorov metric. The convergence takes place in the sense of the weak topology, provided that

(c N0 , ν N0 /N) converges weakly to (c 0 , ξ 0 ) in H 1 (]0, z max [ a ) 2 × M F + (X ). The limit (c, ξ ) satisfies a self-aggregation equation ⟨ξ t , f t ⟩ = ⟨ξ 0 , f 0 ⟩ + t 0 ⟨ξ s , φ (•)ξ s * F a ∂ z f s + ∂ t f s ⟩ds + ∑ k∈E t 0 ∑ ℓ∈T Z ⟨µ k (z, ℓ), f s ⟩q ∞ k q k (z, ℓ, d(ξ s , z), c bs (z))ξ s (dz, ℓ)ds (19) for all f = ( f (1) , • • • , f (|T |) ) ∈ (C 1 (R + × [0, z max ], R)) |T | , where ⟨ξ t , f t ⟩ = ∑ ℓ∈T z max 0 f (ℓ) (t, z)ξ ℓ t (dz)
and ξ ℓ t is the restriction of ξ t to cells of type ℓ. Moreover, for i ∈ {butyrate, O 2 } we have

∂ t c i -σ i ∂ zz c i = -s i γ ∞ β c 4 o c b c 4 o c b + K 5 β (ξ ent t + ξ gc t ) * ψ a (z)
with the same boundary conditions as in [START_REF] Batlle | β -catenin and TCF mediate cell positioning in the intestinal epithelium by controlling the expression of EphB/EphrinB[END_REF][START_REF] Bourgin | Exploring the bacterial impact on cholesterol cycle: A numerical study[END_REF]. The proof of this result is technical and not provided here but can be found in Chapter 4 of [START_REF] Darrigade | Modélisation du dialogue hôte-microbiote au voisinage de l'épithélium de l'intestin distal[END_REF]. This deterministic limit is however quite natural. Indeed, the partial differential operator t 0 ⟨ξ s , φ (ξ s * F a )∂ z f s + ∂ t f s ⟩ds has the same form as the deterministic flow in [START_REF] Buchwald | A review of human salmonellosis: Ii. duration of excretion following infection with nontyphi salmonella[END_REF] where Ãt (c, ν) is replaced by its limit ξ s . For the stochastic part of the model, this limit is a standard result for punctual jump processes (see e.g. Chapter 7 in [START_REF] Bansaye | Stochastic models for structured populations: scaling limits and long time behavior[END_REF]). The fact that this stochastic part is coupled to a PDE-based deterministic flow in a PDMP leads to specific technical issues that are rigorously dealt with in [START_REF] Darrigade | Modélisation du dialogue hôte-microbiote au voisinage de l'épithélium de l'intestin distal[END_REF].

Let Div = {div,sc ; div,pc} ⊂ E denotes the set of division events, Ex = {ex,gc ; ex,ent} denotes the set of extrusion events and Dif = {sc, pc ; pc, gc ; pc, ent} denotes the set of differentiation events. If we assume that for all cellular type ℓ, the restricted measures ξ l t (dz) have a density denoted ρ ℓ (z,t) on [0, z max ], then ( 19) is a weak formulation of the system

∂ t ρ ℓ + ∂ z φ ρ ℓ (ρ * F a )) = ∑ k∈E \Div η k,ℓ q ∞ k q k (z, ℓ, d(ρ, z), c b (z)) (20)
+ ∑ k∈Div η k,ℓ (z -λ (z))q ∞ k q k (z -λ (z), ℓ, (d × ρ * D a )(z -λ (z)), c b (z -λ (z))
where ρ = ∑ ℓ ρ ℓ is the total local cell density and

η k,ℓ =                    ρ ℓ if k ∈ Div, -ρ ℓ if k ∈ Ex, ρ j if k ∈ Dif and the jump type k is a differentiation from cell type j into type ℓ, -ρ ℓ if k ∈ Dif and k is a differentiation from type ℓ into any other type. (21)
We can go one step further, as suggested in [START_REF] Darrigade | Modélisation du dialogue hôte-microbiote au voisinage de l'épithélium de l'intestin distal[END_REF], and derive a formal approximation of the model when the cell size a goes to zero.

We introduce the potential

W a (z) = - z -a
F a (y)dy.

From the definition of F a in (3), we see that

F a (z) = k a k 1 F 1 ( z a ) and therefore W a (z) = a k a k 1 W 1 ( z a ). From the definition of D a in (7), we see that D a (z) = d a d 1 D 1 ( z a ). We assume that k a = k 1 a 2 and d a = 1 a .
If the densities in ( 20) are smooth and using the change of variables ỹ = y/a in the convolutive terms we have after some computation

ρ * F a (z) = -ρ * W ′ a (z) = - a -a ρ ′ (z -y)W a (y)dy = -a 2 k a k 1 1 -1 W 1 ( ỹ)ρ ′ (z -a ỹ)d ỹ = - 1 -1 W 1 ( ỹ)ρ ′ (z -a ỹ)d ỹ and d(ρ, z) = d(z) (ρ * D a ) (z) = d(z) 1 -1 ρ(z -a ỹ)D 1 ( ỹ)d ỹ.
The function d converges pointwise towards 1 when a → 0 on the open interval ]0,

z max [. So ρ * F a (z) -→ a→0 -W ∂ z ρ(z), d(ρ, z) -→ a→0 Dρ(z) with W = 1 -1 y -1 F 1 (z)dzdy and D = 1 -1 D 1 (z)dz. Likewise, c b (z) -→ a→0 c b (z), ρ ℓ * ψ a (z) -→ a→0 ρ ℓ (z) and λ (z) -→ a→0 0.
It follows that when a goes to zero, the self-aggregation equations for the densities have the following formal approximation for all ℓ in T

∂ t ρ ℓ -W ∂ z (φ (z)ρ ℓ ∂ z ρ) = ∑ k η k,ℓ q ∞ k q k (z, ℓ, Dρ, c b ), (22)
These equations are complemented with the following mixed boundary conditions, chosen for their simplicity. They reflect the absence of cells at the bottom of the crypt but stem cells, and zero flux condition at the top.

ρ ℓ (0,t) = 0 if ℓ ̸ = sc, ρ sc (0,t) = ρ bot sc , ∂ z ρ ℓ (z max ,t) = 0.
The coupled reaction-diffusion equations are, for i ∈ {butyrate, O 2 },

∂ t c i -σ i ∂ zz c i = -s i γ ∞ β c 4 o c b c 4 o c b + K 5 β × (ρ ent + ρ gc )
with the same boundary conditions as in [START_REF] Batlle | β -catenin and TCF mediate cell positioning in the intestinal epithelium by controlling the expression of EphB/EphrinB[END_REF][START_REF] Bourgin | Exploring the bacterial impact on cholesterol cycle: A numerical study[END_REF]. Equation ( 22) is a porous media equation (PME) type. More precisely, the nonlinear diffusion term W ∂ z (φ (z)ρ ℓ ∂ z ρ) can be interpreted as a Darcy law for incompressible fluids in a porous media (according to which the advection speed is proportional to the local pressure gradient) where the pressure is proportional to the local total density of cells. Such equations have already been suggested to model densities of interacting cells [START_REF] Lushnikov | Macroscopic dynamics of biological cells interacting via chemotaxis and direct contact[END_REF][START_REF] Murray | From a discrete to a continuum model of cell dynamics in one dimension[END_REF][START_REF] Murray | Comparing a discrete and continuum model of the intestinal crypt[END_REF]. The equations present a non linear cross-diffusion term coming from the presence of ρ(z,t) = ∑ ℓ∈T ρ ℓ (z,t), the total cell density. The source terms are formally derived from the expressions given above for the self-aggregation equations. The spatial shift λ for introduction of daughter cells during cell division now vanishes in the deterministic limit, since the cell radius is assumed to be small. In other words, daughter cells are introduced at the same place as mother cells during division in the deterministic model.

Numerical implementation

The porous media type equations [START_REF] Curtain | An introduction to infinite dimensional linear systems theory[END_REF] are not trivial to solve, due to the non linear cross-diffusion terms, and requires adapted numerical tools. The kinetic diffusive explicit scheme described in [START_REF] Aregba-Driollet | Explicit diffusive kinetic schemes for nonlinear degenerate parabolic systems[END_REF] and used in [START_REF] Twarogowska | Numerical approximation and analysis of mathematical models arising in cells movement[END_REF] designed for nonlinear transport diffusion equations allowed us to obtain satisfactory simulations. Therefore, letting ρℓ = ∑ k̸ =ℓ ρ k and G(z) = 1 2 W φ (z), we modified ( 22) to obtain the system of transportdiffusion equations with source terms

∂ t ρ ℓ -∂ z W φ ∂ z ( ρℓ )ρ ℓ + G ′ ρ 2 ℓ -G∂ zz ρ 2 ℓ = -G ′′ ρ 2 ℓ + ∑ k η k,ℓ q ∞ k q k (z, ℓ, D × ρ, c b ) (23
) which is the model structure needed for the kinetic diffusive explicit scheme. The python code is available at https://forgemia.inra.fr/decrypt/deterministic-model.

The deterministic approximation is fastest to compute and allows us to perform more intensive numerical explorations than with the PDMP. On the other hand, the stochastic model is relevant to extract features directly comparable with biological observations, such as migration speed, cell fate tracking or renewal of the epithelium. In the sequel, we will first explore the PDMP model and show that it renders accurately the main features of an healthy crypt in term of population distribution and turn-over. Then, we will verify that the deterministic limit [START_REF] Curtain | An introduction to infinite dimensional linear systems theory[END_REF] provides an accurate approximation of the PDMP model for sufficiently small cells and we will carry out a systematic sensitivity analysis of the deterministic limit using the Elementary Effects method of Morris. This sensitivity analysis would not have been tractable with the original PDMP model due to an excessive computational load.

Numerical exploration of the PDMP model

We now show that simulations reproduce qualitatively the behaviour of a distal colonic mouse crypt. To understand the behaviour of our model, we did not start by simulating the full model presented in Sect. 2. We rather conducted a sequence of numerical experiments on submodels of increasing complexity, adding at each step new components from the full model. In 4.1, we explore the sole effect of cell-cell mechanical interactions. In 4.2, we reproduce the migration of epithelial cells by adding spatially regulated division and extrusion to the mechanical interactions. In 4.3, we explore the spatial distribution of cell populations by adding cell types and spatial regulation of differentiation. In 4.4, we add butyrate and oxygen diffusion and consumption by epithelial cells, and finally, in 4.5, we compute the whole model by adding butyrate regulation on cell fate.

Cell layers formation due to mechanical interactions

We begin our experiments with a population of identical cells subject only to mechanical interactions, with no stochastic events whatsoever. Parameters values are given in Table 7. Starting with a uniform cell distribution in the crypt, we observe that cells rapidly gather into layers. This phenomenon is allowed by the specific shape of the interaction function F a , as F a (0) = 0. As expected, the spacing between layers is influenced by the cell diameter a which is the unique parameter of the interaction kernel (Fig. 3).

Migration effect

Migration of epithelial cells is a hallmark behaviour of the intestinal crypt. As explained in Sect. 1.1, epithelial cells are produced by division of undifferentiated cells in the bottom third of the crypt and migrate to the upper part of the crypt where they are extruded. We can reproduce this phenomenon by simply adding to the previous sub-model stochastic division and extrusion that are regulated spatially and by density.

We consider a population of identical cells, interacting mechanically through Eq. ( 5), which divide and are extruded with respective rates q div (•) and q ex (•):

q ∞ div q div (z, ν * D a (z)) = q ∞ div regulation by position R(z, K div [z], κ div [z]) × regulation by density R(ν * D a (z), K div [dens], κ div [dens]), (24)
q ∞ ex q ex (z, ν * D a (z)) = q ∞ ex R(z, K ex [z], κ ex [z])
regulation by position

× (1 -R(ν * D a (z), K div [dens], κ div [dens])) regulation by density . ( 25
)
Figure 4a shows the trajectory of cells in the crypt, during one realization of the PDMP. The initial positions of 700 cells are drawn from a spatial uniform distribution, as in Sect. 4.1, and parameter values are given in Table 8. With the addition of spatially regulated division and extrusion events, we see that cells migrate from the bottom to the top of the crypt. In the lower zone of the crypt (z ≤ 75µm), cell division occurs. When a new cell is introduced at a distance a/2 of the mother cell, it rapidly joins a close cell layer as illustrated in the previous Sect. 4.1. This corresponds to the comb-like patterns in Fig. 4a. Figure 4b shows the number of division and extrusion events according to the position in the crypt, for the same realization as Fig. 4a. As in normal crypts, division occurs at the bottom of the crypt and extrusion at the top, fueling the migration phenomena. We finally display the time evolution of the number of cells, and the mean over 50 realizations of the PDMP in Fig. 4c. We can observe that the total population reaches a stationary state of around 740 cells after about 200 hours, with little stochastic variations around this value. This stationary state indicates an equilibrium between division and extrusion.

Interesting biomarkers evaluated during experiments to quantify the crypt activity are the total cell population (TCP), the average migration speed (MS) and the individual division rate (DR). In our model, possible drivers of these quantities are the interaction intensity k a , the extrusion rate q ∞ ex and the sensitivity to cell density of cell division K div [dens]. We conducted a local sensitivity analysis of the model to these parameters by changing their value by ± 30 and 50 % respectively to the reference value given in Table 8. For each condition, we ran 50 simulations starting from 700 cells drawn in a uniform law on [0, z max ] and followed the evolution of TCP, MS and DR at steady state. For one set of parameters, to compute MS, we computed the time needed by a cell to cover the distance between z = 40 and z = 175 µm to avoid boundary effects. Namely, let t f be the final time of simulation and z i (t) the position of cell i for t ∈ [b i , d i ] with b i and d i respectively the birth and death time of cell i (and b i = 0 for a cell present at t = 0, and d i = t f for cells still alive at t = t f ), we defined the start time

u i = inf{t ∈ [b i , d i ] : z i (t) ∈ [40, 175]} and the exit time e i = sup{t ∈ [b i , d i ] : z i (t) ∈ [40, 175]}.
Eventually, setting T i = e iu i and Z i = z i (e i )z i (u i ), the migration speed of cell i was defined as the ratio Z i T i . Note that this speed does not necessarily match the fastest or the lowest speed reached by cell i during the simulation. The division rate of cell i is n div,i d i -b i with n div,i the number of divisions of cell i. We then computed the average of these values among cells at steady state (after t = 300 h) for one realization of the PDMP, and then averaged over the 50 realizations to obtain respectively MS and DR. The TCP is taken at final time and averaged over the 50 realizations of the PDMP.

The effect of the three parameters tested is shown in Fig. 5. We observe that TCP is particularly impacted by the sensitivity of the division rate to local cell density: it shows a linear dependence to the K div [dens] parameter. In contrast, the extrusion rate and the repulsion force intensity have much lighter, and negative, effect.

On the contrary, the three parameters show a very similar influence on MS. Increasing the repulsion force, the extrusion rate or the division sensitivity all result in a comparable increase of cell velocity. A marked drop of the migration speed is however observed when the sensitivity of the division to local cell density is strongly reduced due to a specific behavior explained in Appendix D. The mechanisms by which the migration speed is increased seem however different. Increasing the repulsion force directly impacts mechanical interactions all along the crypt, so that cells move faster when they are not at mechanical equilibrium. By contrast, a higher K div [dens] results in a higher division rate leading to a higher cell density in the bottom of the crypt and thus pushing the whole cells upward. A higher q ∞ ex has a similar effect by depleting the cell density in the upper crypt. The greater number of cells for higher values of K div [dens] also increases the intensity of mechanical interactions. DR exhibits the same dependence to the three parameters as MS. This is expected because at steady-state, the total number of cells is almost constant (Fig. 4c) and therefore divisions must keep up with the migration to maintain the population. Therefore, q ∞ ex and k a seem to increase the activity of the epithelium, i.e. the migration speed and the number of divisions, while maintaining a constant population. On the contrary, K div [dens] strongly impacts the total cell population together with the epithelial activity.

Spatial distribution and cell turn-over in the crypt

We now add cell types to the previous model. The 5 different cell types and cell differentiation events described in Sect. 2 are introduced, with the objective to reproduce the population size of each epithelial cell type in the crypt as well as their spatial distribution, as explained in Sect. 1.1 and depicted in Fig. 1a. We assume that there is a uniform concentration of butyrate in the crypt, so that cell differentiation rates only depend on spatial and density terms. The precise jump rates we use are detailed in Appendix A.1 and summarized in Table 2. Parameter values are given in Table 9. We ran 50 independent simulations starting from an initial condition where cells are already spatially organised according to their types. For each simulation, positions of the 12 DCS cells are drawn randomly according to the following linear decreasing probability density

ρ dcs (z)dz = 2 12 (1 -z 12 ) if z ≤ 12, 0 otherwise. ( 26
)
We can observe that the population of each cell types rapidly reaches a steady state after about 4 days (Fig. 6a). These steady state population levels match particularly well available biological data for the distal colonic crypt of mice (see Table 3). Total cell number, number of stem cells, and cell type ratios are correctly recovered, showing that our model is able to capture the phenomena insuring the population level regulations. Moreover, the spatial distribution of each cell types in the crypt also reproduces the typical spatial organization of a crypt described in Sect. 1.1 (Fig. 6b). The stem cells remain in the distal part of the crypt (z ≤ 20µm), replaced upward by progenitor cells that are mainly located under z ≤ 80µm. Then progenitor cells are progressively replaced for z ≥ 50µm by goblet cells and enterocytes, with the approximate ratio 1 : 3. DCS cells are kept in the stem cell niche (z ≤ 20µm).

We also evaluate the renewal of the epithelium. To do so, we mimic the experimental procedure which relies on labeling mitotic cells and tracking their elimination over time [START_REF] Darwich | Meta-analysis of the turnover of intestinal epithelia in preclinical animal species and humans[END_REF]. We label at steady-state (t = 800h) the stem and progenitor cells, and we screen over time the ratio of labeled cells that have been eliminated respectively to the number of initially labeled cells. Upon division of a labeled cell, the label is randomly passed to the mother or the daughter cell with probability p = 0.5. We observe a fast renewal of epithelial cells (Fig. 7): after 2 days, on average 66% of labeled cells have been extruded, and it rises to 95% after 3.6 days. This is in agreement with experimental values for the renewal time of the mouse colonic epithelium which is 2.67 ± 2.04 days [START_REF] Darwich | Meta-analysis of the turnover of intestinal epithelia in preclinical animal species and humans[END_REF]. The migration speed is below reported experimental value. However, the reference values we use here for our observables come from different experiments which might explain some of the discrepancies.

We also note that the proportion of labeled cells that are eventually extruded seems to reach a plateau above 95% in our model, indicating that a small part of the cells retain their label and divide rarely. Indeed, a portion of the stem cells in the lowest part of the crypt have a very low division rate due to the contact inhibition exerted by neighbouring cells. Moreover, these cells survive for a very long time as they stay at the bottom of the crypt where they move very slowly due to the curvature (the φ (z) term in the force F i (t) in ( 2)) and the pressure exerted by fast dividing Regulation pathway ( j) Table 3: Simulations results and experimental reference for the submodel of Sect. 4.3, mean of 50 simulations, t f = 1000h. SI = small intestine.

Jump index (k) Wnt [z] density [dens] div, sc - - div, pc - - sc, pc + ∅ pc, gc + ∅ pc, ent + ∅ ex, gc + + ex, ent + +
progenitor cells located above them. We can see that about 17% of the stem cells are renewed (Fig. 7, blue curve) if we only label the stem cells, and that the average division rate of stem cells in our simulations is lower than for in vivo observations (Table 3). However, the latter lays within the simulated standard deviation interval. Progenitors average division rate is higher that experimental observation, but closer than for stem cells. As stem cells and progenitor cells influence each other, this is partially explained by the low division rate of stem cells. The variation of the stem cells renewal rate (Fig. 7) can be partly explained by the variation of the DCS cells positions. Indeed, at the beginning of each simulation these positions are drawn randomly in the distribution [START_REF] Degirmenci | GLI1-expressing mesenchymal cells form the essential Wnt-secreting niche for colon stem cells[END_REF]. DCS cells will then both impact contact inhibition of the stem cells and protect the stem cell niche from an invasion by the rapidly dividing progenitor cells accordingly.

4.4 Reproduction of the butyrate and O 2 gradient.

In the previous section, the spatial distribution of epithelial cells in the crypt was reproduced based on physical interactions between cells and spatial regulation, discarding cell metabolism. We now add butyrate and oxygen concentrations, whose dynamics follow equations (6,9) and (6,10) respectively, including the cell consump-tion in the source terms. The corresponding parameters are given in Table 10. The derivation of the chemical kinetics parameters is detailed in Appendix E.1. We do not yet take into account the feedback of butyrate on cell fate events (that is to say we use the exact same jump rates as in Sect.4.3, see Table 2 and Appendix A.1). At first, we aimed at reproducing the gradient of dioxygen and butyrate observed in the crypt at homeostasis. We ran 50 independent simulations starting from an initial condition where the number of cells from each cell type corresponds to the average number obtained at steady-state in Sect. 4.3, and position of cells were randomly drawn from the average steady-state distribution of their cell type obtained in Sect. 4.3 (and plotted in Fig. 6). Initial concentrations of butyrate and oxygen were respectively linearly increasing and decreasing with z. Concentrations of butyrate and O 2 obtained are plotted in Fig. 8. There is a strong spatial gradient of concentration for both molecules, due to the consumption by differentiated cells in the upper part of the crypt. Oxygen concentration decreases with z from 10 au at the bottom of the crypt to 1 au at the top while butyrate concentration increases from around 1 mM in the stem cell niche to 5 mM at the top. Few validation criteria are available for these results as concentrations cannot yet be measured precisely experimentally in the crypt, despite being considered to play a major role in crypt organization. We used the fact that concentration of O 2 at the top of the crypt is around 10 times lower than at the bottom, as reported in [START_REF] Rivera-Chávez | Oxygen as a driver of gut dysbiosis[END_REF]. This ratio of 10:1 is what we observe in our simulations, along with a concentration gradient of butyrate. These results confirm the hypothesis that the metabolism of differentiated epithelial cells in the upper part of the crypt induce a concentration gradient of butyrate at homeostasis and protect stem cells from deleterious effects of butyrate [START_REF] Kaiko | The Colonic Crypt Protects Stem Cells from Microbiota-Derived Metabolites[END_REF][START_REF] Snippert | Colonic Crypts: Safe Haven from Microbial Products[END_REF]. In vitro experiments have shown that a 1 mM concentration of butyrate is sufficient to strongly inhibit stem cells and progenitors proliferation [START_REF] Kaiko | The Colonic Crypt Protects Stem Cells from Microbiota-Derived Metabolites[END_REF]. In our model, the concentration of butyrate in the stem cell niche is 0.9 mM. To explain this value close to the inhibition value, one can argue that the experiments we mention were conducted in vitro. On the other hand, our model might overestimate the concentration of butyrate as we did not take into account the ketone production which is another metabolic pathway of butyrate consumption that does not involve oxygen [START_REF] Cherbuy | Effect of germfree state on the capacities of isolated rat colonocytes to metabolize n-Butyrate, glucose, and glutamine[END_REF].

Complete model

Finally, we simulated the complete IBM described in Sect. 2, which takes into account the effect of butyrate on the rates of stem cell division and progenitor differentiation as summarized in Table 1. We assessed that the capability to reproduce the homeostatic crypt is conserved, regardless of the additional mechanism that is added. We ran 50 independent simulations using the same initial condition as in Sect. 4.4. In figure 9, we plotted the trajectories of some individual cells. We observe the expected behaviour of birth and proliferation at the bottom of the crypt, followed by differentiation and then extrusion along their migration to the top of the crypt. Figure 10 represents cell repartitions, numbers, and butyrate and oxygen concentrations. Figure 10c represents the renewal rate of the epithelium. Simulation results are summarized in Table 4. We observe a good recovery of the total number of cells, the different cell population ratios and the concentration gradients (Fig. 10 and Table 4), except for the number of progenitors which is higher than expected. The division rate of progenitors is still in the experimental range, and interestingly, the division rate of stem cells is improved compared to Sect. 4.3 and within experimental range (Table 4). The renewal of the epithelium (Fig. 10c) is however slower than in Sect. 4.3 and than experimental values. This might be due to the excess of progenitor cells which increases cell density in the progenitor zone and therefore decreases the division rate. Now, we investigate the effect of a lower availability of butyrate in the lumen, which mimics a dysbiotic situation where butyrate-producing bacteria are reduced. We changed the butyrate concentration at the top of the crypt c b,lum , testing values from 0.1 to 5 mM (Fig. 11b). In this case, the O 2 concentration at the top of the crypt rises and doubles for a luminal butyrate concentration of 1 mM, meaning that the luminal part of the crypt becomes microaerobic. This kind of dysbiotic situation where butyrate concentration is low, oxygen is high and cell type composition of the epithelium is similar to homeostasis may occur after an antibiotic treatment [START_REF] Rivera-Chávez | Oxygen as a driver of gut dysbiosis[END_REF].

All together, the complete model includes the main regulation mechanisms of crypt turn-over and renders accurately most of the main biomarkers of an healthy crypt at homeostasis. However, at this level of complexity, the IBM becomes challenging to manipulate and the deterministic PDE limit model turns out useful.

Numerical exploration and analysis based on the deterministic approximation

Accuracy of the deterministic approximation

We aim to compare the outputs of the complete PDMP model with the deterministic approximation ( 22) derived in Sect. 3. The crypt homeostasis being described by the model steady states, we compare complete PDMP and its deterministic approximation when an equilibrium is reached. Starting from an empty crypt (i.e. ρ l,0 = 0 ∀l ∈ T ) and from a DCS cells distribution that approximates the distribution observed in the PDMP, the PDE approximation is simulated with a given boundary condition for the stem cells (see below) until steady state and compared to the PDMP outputs.

To approximate the DCS cells density found in the PDMP model, we chose the following shape :

ρ dcs (z) = N dcs z max 0 ρdcs (z)dz ρdcs (z) (27)
with

ρdcs (z) =          d(z -z d ) + 1 if z d -1 d < z < z d , 1 if z d ≤ z ≤ z u , 1 -u(z -z u ) if z u < z < z u + 1 u , 0 otherwise,
where 0 ≤ z d ≤ z u , d, u > 0 and N dcs is the predefined number of DCS cells. This function is a piecewise-linear approximation of a step-function summing to N dcs , controlled by 4 parameters (d, u, z d , z u ). As the PDMP is far from the asymptotic condition (cell size close to zero) for which the approximation was derived, we adapted the values of these parameters, and additionally the values of W , D and γ ∞ β as well as the boundary value for the stem cell density at the bottom, to obtain a reasonable fit of the mean cell distribution and concentration profiles predicted by the PDMP at steady state.

We observed that, although the DCS cells are "inactive" their shape and position influence cell distributions and proportions. When varying the values of (d, u, z d , z u ) which mimics the random draws of DCS cells in the PDMP model simulation, we observed that some parameter sets produced undesired behaviours such as no cell division at all. This is consistent with the high variability of stem cells renewal that we observed in the IBM and that we attributed to the randomness of the DCS cells position. In certain cases, DCS cells blocked stem cells division by contact inhibition (see end of Sect. 4.3).

The parameters that provide a suitable approximation of the PDMP model are gathered in Table 12. They will be used for numerical simulation in the sequel. To compare PDMP and PME model behaviours, we display steady-state densities of cells for both models in Fig. 12. We furthermore display the comparison of cell fate events by plotting the right-hand side of Eq. ( 22), i.e. the spatialized rate of division, differentiation and extrusion at steady-state. The spatialization of the cell fate events of the PDMP is correctly recovered, with extrusion events at the top of the crypt, division events near the stem cell niche in the lower crypt, and differentiation events distributed along the first half of the crypt. Spatial distribution of cells and cell fate events of PME are in good agreement with the PDMP.

Limit model outputs are gathered in Table 5 and can be compared with the complete PDMP outputs in Table 4. We evaluate the average rate of the different cell fate events according to the following formula

τ k = z max 0 ρ ℓ (z)q ∞ k q k (z, ℓ, D × ρ, c b )dz z max 0 ρ ℓ (z)dz (28)
for (k, ρ ℓ ) ∈ {((div, sc), ρ sc ), ((sc, pc), ρ sc ), ((div, pc), ρ pc ), ((pc, gc), ρ pc ), ((pc, ent), ρ pc ), ((ex, gc), ρ gc ), ((ex, ent), ρ ent )}. We observe that total cell and stem cells numbers are in the range of the PDMP results (709 vs 694 ± 22, and 20 vs 17 ± 8), as well as the ratio of globet cells/enterocytes ( 0.3 vs 0.32 ± 0.03) or stem cell division rate (0.041 vs 0.021 ± 0.013 h -1 ). Although progenitor cells division rate is in the right order (0.041 vs 0.05 ± 0.0 h -1 ) the ratio between progenitor and total cell numbers is lower (0.14 vs 0.26 ± 0.02). The parameters at our hand (DCS density shape, W , D, γ ∞ β and boundary condition for the stem cells) do not allow an appropriate value while keeping a good fit for the other observables in Table 5. Our sensitivity analysis in Sect. [START_REF] Bansaye | Stochastic models for structured populations: scaling limits and long time behavior[END_REF] 28) for the formula used to compute rates.

Sensitivity analysis

We use the deterministic approximation of the IBM to perform a Sensitivity Analysis (SA) using the global method of the Elementary Effects (EE) of Morris [START_REF] Morris | Factorial Sampling Plans for Preliminary Computational Experiments[END_REF][START_REF] Saltelli | Why so many published sensitivity analyses are false: A systematic review of sensitivity analysis practices[END_REF].

Indeed, this method uses variance-based techniques and exhibits reasonable computational time, which makes it suitable for SA on models with high number of inputs or with high computational load [START_REF] Campolongo | An effective screening design for sensitivity analysis of large models[END_REF]. We briefly recall that Morris method estimates sensitivity indices by generating random trajectories in the rescaled parameter space Ω , such that along each trajectory, parameter values are changed one by one following a predefined variation ∆ . More precisely, let X = (X 1 , • • • , X p ). For each trajectory k, we draw a random perturbation of the set of integers {1 ≤ i ≤ p} denoted by {σ k (i)}. Then, the i-th step of the trajectory k is, ∀1 ≤ i ≤ n,

X k i = (X k 1,0 + 1 {σ k (1)≤i} × ∆ , • • • , X k p,0 + 1 {σ k (p)≤i} × ∆ ) and X k 0 = (X k 1,0 , • • • , X k p,0
). According to [START_REF] Morris | Factorial Sampling Plans for Preliminary Computational Experiments[END_REF], we used ∆ = 2/3. Each trajectory k allows to calculate the elementary effect EE k i of the parameter X i on the output Y

EE k i = Y (X k σ (i) ) -Y (X k σ (i)-1 ) ∆ . (29)
A sample distribution of EE values for each parameter is generated by simulating r trajectories with random starting points, from which the following sensitivity indices are computed for each parameter X i

µ i = 1 r r ∑ k=1 EE k i ; σ i = 1 r r ∑ k=1 (EE k i -µ i ) 2 ; µ * i = 1 r r ∑ k=1 |EE k i |. (30)
Indices are the mean µ i , the standard deviation σ i and the index µ * i proposed in [START_REF] Campolongo | An effective screening design for sensitivity analysis of large models[END_REF] to avoid compensation of effects with opposite signs in the case of a non-monotonic influence. Following [START_REF] Awad | Convergence of sensitivity analysis methods for evaluating combined influences of model inputs[END_REF] (see Table 6) parameters can subsequently be classified by their effect (negligible, significant, non linear...) on the output of interest.

Effect Morris index

Decreasing µ < 0 Increasing µ > 0 Negligible µ * ∼ 0 and σ ∼ 0 Linear and additive µ * high and σ ∼ 0 Non linear or non monotone σ high In interactions of order ≥ 2 σ high Table 6: Morris index interpretation [START_REF] Awad | Convergence of sensitivity analysis methods for evaluating combined influences of model inputs[END_REF] The output of the model selected for the SA were the total density of cells in the crypt, the division rate of progenitor cells and the ratio of O 2 concentration at the bottom and at the top of the crypt, as for the PDMP. To carry out the analysis we vary the parameters over a range of 50% of their value described in Tables 11 and12 except for W which varies at 30%. Moreover, to select the number of trajectories r in order to obtain stable results, we applied two criteria given in [START_REF] Sarrazin | Global Sensitivity Analysis of environmental models: Convergence and validation[END_REF]. The first criteria (Ranking) is close to 1 when the ranking of influential parameters is stabilized and the second (Screening) reaches values below 0.05 when the screening (separation between influential and non-influential parameters) is stable. The authors of [START_REF] Sarrazin | Global Sensitivity Analysis of environmental models: Convergence and validation[END_REF] propose another criteria based on indices values, that we discarded because we are interested in large effects. We used two Python libraries, SALib and SAFE to conduct the SA [START_REF] Herman | SALib: An open-source Python library for Sensitivity Analysis[END_REF][START_REF] Pianosi | A Matlab toolbox for Global Sensitivity Analysis[END_REF].

Results of the sensitivity analysis are presented in Fig. 13. In the majority of cases the ranking and screening criteria thresholds are reached after running 200 trajectories. We plot the standard deviation of EEs σ as a function of µ * , the mean of the absolute values. We recall that high values of µ * indicate a high impact on the output, while high values of σ is representative of non-monotone or non-linear effects, or possible interaction with another parameter. Fig. 13 allows a better understanding of the effects of parameters on three selected outputs, described below. For all of them, we can see that the most influential parameters also exhibit high values of σ , hence complex effects.

The total density is mainly influenced by 2 groups of parameters. The first one gathers the most influential parameters. It is composed of regulation parameters with a positive effect (K div,pc [dens] as for the IBM of Sect. 4.2, K ex,gc [z] and K div,pc [z]). Parameters regulating progenitor cells division are particularly involved as progenitor cells are the main type of proliferating cells. Parameter W , which characterizes repulsion between cells, is also part of this group but it has a negative effect which is consistent with the influence of k a on the PDMP model in Fig. 5. In the second group of still important but less influential parameters, we find parameters related to stem cells and concentration boundary conditions.

The division rate of progenitor cells is mainly influenced by W and K div,pc [z], results in line with those obtained for total density. When K div,pc [z] increases, progenitor cells differentiate higher in the crypt and therefore epithelial cells spend more time in a progenitor state. However, they don't make more divisions as most likely the epithelium is saturated and density inhibition is fully active. From there comes the negative impact of K div,pc [z] on the division rate. We also note that an increase of the butyrate concentration at the top of the crypt c b,lum has a positive effect on the division rate. This could be explained by a faster differentiation (stimulated by butyrate, see Table 1) of progenitor cells into goblet cells or enterocytes. As expected, stem cell parameters have a negative effect on progenitor cell due to the competition for space between progenitors and stem cells.

The parameters with highest impact on the ratio of O 2 concentration at the bottom and at the top of the crypt are those related to the diffusion equations, as it could be expected. We find in first place boundary conditions c o,bot and c b,lum that both have positive effect. An augmentation of butyrate concentration at the top of the crypt increase the consumption of oxygen and lead to the increase of the ratio. Parameter K β has a negative effect. Indeed, this parameter models an activation threshold in the reaction term and therefore limits the oxygen consumption when its value is higher, which results in a decrease of the ratio. Interestingly, this output is also influenced by parameters linked to cell fate, namely those controlling the extrusion rate of differentiated cells in space and density (K ex,gc [z], and to a lesser extent K ex,gc [dens]) as well as those controlling the proliferation of progenitor cells in density (K div,pc [dens]) and the transition between proliferation and differentiation of progenitor cells in space (K div,pc [z]). The sign of these effects has to be considered with care, due to the high values of σ .

Additionally, Fig. 14 summarizes the sensitivity analysis for 20 different outputs : total density (ρ), density of each cell type (ρ ℓ ), ratio of progenitor cells to total density, medium position ( zℓ ) given by Eq. ( 31), rate of cell fate events, concentration of O 2 (resp. butyrate) at the top (resp. bottom) of the crypt, ratio of O 2 concentration (at the bottom / at the top). To visualize results for all the outputs on the same figure we normalized for each output the values µ * i of each parameter X i by max i (µ * i ). We also indicated with a color code if their influence on the output is positive or negative according the the sign of µ i . zℓ =

z max 0 ρ ℓ (z,t)zdz z max 0 ρ ℓ (z,t) dz (31)
One parameter in particular emerges from all these analyses, W , which is highly influent for a majority of the outputs presented. Other parameters are also influential such as

K div, [z], K div,pc [z], K div,pc [dens], K ex,gc [z].
6 Exploring crypt-microbiota-pathogen interactions for treatment design It has been evidenced that several intestinal pathogens hijack the Wnt signaling pathway [START_REF] Rogan | Bacterial Manipulation of Wnt Signaling: A Host-Pathogen Tug-of-Wnt[END_REF]. Salmonella enterica and Citrobacter rodentium -which is a mouse pathogen used as a model for human ulcerative colitis and infection by enteropathogenic and enterohemorrhagic Escherischia coli -both increase Wnt signaling in the crypt epithelium of the mouse, which promotes proliferation of undifferentiated cells. Additionally, these two pathogens increase respiratory electron acceptors concentrations in the lumen, including O 2 , to fuel their own growth [START_REF] Rivera-Chávez | Oxygen as a driver of gut dysbiosis[END_REF]. These infections also impact the microbiota structure. Infection with Salmonella has been shown to deplete butyrate producers of the Clostridia class. The presence of butyrate producers and high butyrate levels increase resistance of mice to Citrobacter infection [START_REF] Osbelt | Variations in microbiota composition of laboratory mice influence Citrobacter rodentium infection via variable short-chain fatty acid production[END_REF], and human patients suffering from ulcerative colitis show decrease expression of genes related to butyrate metabolism [START_REF] Zhou | Evidence for the butyrate metabolism as key pathway improving ulcerative colitis in both pediatric and adult patients[END_REF]. We highlight in this section the insights brought by our model in this context, by checking the ability of our model to reproduce these pathogenic situations in section 6.1 and by addressing treatment design issues in section 6.2.

Pathogenic steady states

Our model can be readily used to mimic a disruption of butyrate or Wnt equilibrium. All other parameters being unchanged, we set a lower concentration of butyrate c b,lum = 0.5 mM (instead of 5 mM for a healthy crypt) or a higher level of Wnt signaling K div,pc [z] = 160 µm (instead of 40 µm). The two steady states we obtain are plotted respectively in Fig. 15a and Fig. 15b. The two situations are modelled with both the PDMP and the PME, with again a very good agreement even in these pathological cases. In both situations, we observe an increased number of undifferentiated cells and luminal O 2 concentration compared to the healthy state. The phenotype associated to low butyrate, where luminal O 2 concentration is the highest and there is a moderate increase of undifferentiated cells is similar to what is observed in a Salmonnella infection. In the case of high Wnt signaling, the increase of luminal O 2 concentration is less pronounced but the crypt is almost totally filled with undifferentiated cells. This is consistent with the phenomenon of epithelial hyperplasia observed in the case of a Citrobacter infection. In this latter case, we also observe a strong increase of butyrate concentration at the bottom of the crypt, which could be deleterious for stem cells. In vivo, this increase could be transient (or even unobserved) as changes triggered by the infection (including the increase of luminal O 2 ) could kill butyrate producers in the microbiota, and this cannot be taken into account here as there is no dynamical model for the microbiota. This suggests the need to investigate in more details and exploit the interactions between the pathogen, the crypt cells and the microbiota, which is done in the next section.

Investigating a treatment by butyrate supplementation

We have seen in the previous section that infections by Salmonella or Citrobacter are characterised by a higher level of undifferentiated cells in the epithelium and a higher level of luminal O 2 than in the healthy case. A vicious circle can set in with the microbiota: luminal O 2 is used by the pathogen for its growth, and in the meantime kills the butyrate producers, which in turn reduces butyrate production and therefore O 2 consumption and differentiation of epithelial cells. However, we may wonder if this dynamic can be mitigated by a supplementation of butyrate in the lumen. We therefore address the following questions: 1) Does butyrate supplementation allow to accelerate recovery after enteric infection ? 2) Is recovery impacted by dosage and intake time schemes ?

To investigate this treatment strategy, we coupled the crypt model with a very simple representation of the microbiota and pathogen in the following way. Since the presence of O 2 in the lumen is toxic to butyrate producers and is generally a marker of dysbiosis, we consider that the fraction of active butyrate producer b(t) in the microbiota depends on the level of O 2 at the top of the crypt (c o,lum )

db dt (t) = µ b 1 {c o,lum (t)≤c * o } -1 {c o,lum (t)>c * o } (b max -b(t))(b(t) -b min ). ( 32
)
The pathogen p(t) fills the ecological niche left by butyrate producers:

p(t) = p max 1 - b(t) b max ( 33
)
where p max is the maximal quantity of pathogen. The luminal butyrate c b,lum (t) can either be produced by the butyrate producers (c b,endo (t)) or come from a treatment (c b,treat (t)):

c b,lum (t) = c b,endo (t) + c b,treat (t) with c b,endo (t) = k b b(t). (34)
Eventually, we consider that the pathogen, above a certain level p * , triggers an increase of Wnt signaling from the host:

K div,pc [z](t) = K div,pc [z] min if p(t) ≤ p * , K div,pc [z] max otherwise. ( 35
)
Equation ( 34) and ( 33) can be used to eliminate the butyrate producers and the pathogen from the model, so that, considering only the dynamic of butyrate and Wnt signaling, we obtain:

c b,lum (t) = c b,endo (t) + c b,treat (t), ( 36
) dc b,endo dt (t) = µ 1 {c o,lum (t)≤c * o } -1 {c o,lum (t)>c * o } (c b,max -c b,endo (t)) ×(c b,endo (t) -c b,min ), (37)
K div,pc [z](t) = K div,pc [z] min if c b,endo (t) ≥ c * b , K div,pc [z] max otherwise. ( 38
)
We used both butyrate and Wnt signaling dysbiotic levels introduced in previous section to set the model parameters values to c b,min = 0.5 mM and K div,pc [z] max = 160 µm so that this situation corresponds to a fully infected crypt with low endogenous butyrate concentration and high Wnt signaling. The corresponding high steadystate value of oxygen at the top of the crypt is c o,max = 3.5 au. We set c b,max = 5 mM and K div,pc [z] was arbitrarily set to 0.01 h -1 so that the evolution of the endogenous butyrate shows a realistic response time of approximately 10 days. This is consistent with the time scale of an infection by the enteric pathogens since the infection cycle of Citrobacter lasts two to three weeks [START_REF] Collins | Citrobacter rodentium: infection, inflammation and the microbiota[END_REF], and Salmonella excretion after infection lasts three weeks in average [START_REF] Buchwald | A review of human salmonellosis: Ii. duration of excretion following infection with nontyphi salmonella[END_REF].

Starting from the fully infected crypt, we tested different simple treatment strategies for c b,treat (t) which might move the system towards a healthy state:

butyrate supplementation for one hour every ten hours at concentrations of 2, 6.5, 11, 15.5, 20 mM, butyrate supplementation at 2 mM for 3.25, 5.5, 7.75 hours every ten hours, or a constant supplementation.

The simulation spanned 400 hours and the treatment was applied during the first 300 hours. The strategy consisting in drastically increasing butyrate concentrations for short time periods (1 h) was never successful within the range of concentrations that we tested (see Fig. 16). Interestingly, the only treatments that overcame the infection were the periodic butyrate supplementation for 7.75 hours every 10 hours or the constant supplementation at 2 mM (see Fig. 17). When butyrate is supplemented, O 2 concentration drops due to the metabolic activity of differentiated cells. As O 2 drops below c * o during supplementation, its allows for the growth of butyrate producers and the rise of endogeneous butyrate. It also increases the amount of differentiated cells. This increase seems almost negligible but, as differentiated cells consume a lot of butyrate and oxygen, it can still impact the concentrations. When endogeneous butyrate reaches a threshold of c * b (after 200 h), the epithelium switches back to a healthy state: Wnt signaling instantaneously decreases which induces a massive production of differentiated cells. This in turn decreases O 2 concentrations which drops below c * o and butyrate producers of the microbiota grow back to their healthy level.

Discussion and perspective

A model of cell fate in the intestinal crypt, including interactions with microbiotaderived metabolites has been introduced and studied. We now discuss the new features included in the model. We finally raise and discuss model limitations.

A PDMP model of crypt dynamics, including interaction with the microbiota

The model hereby introduced adds up to a long list of models in the literature studying the cell fate in the intestinal crypt. It however improves existing models in two ways. First, it provides new biological insights by focusing on interactions with the gut microbiota. Indeed, existing models mainly address issues related to morphogenesis or oncology: they then put a strong effort on modeling mechanical and biological mechanisms involved in cell turn-over and crypt formation. In this study, we rather focus on regulations of the crypt homeostasis that involve important interactions with the microbiota. These interactions, secondary for morphogenesis or abnormal crypt development during cancer, are crucial at homeostasis: in an healthy crypt, microbederived metabolites such as butyrate are used by enterocytes as an energetic source, the metabolism of which consumes oxygen and maintains anaerobia in the intestinal lumen. This microaerobic environment is beneficial for butyrate-producing microorganisms. In our model, the butyrate concentration presents a strong gradient in the crypt and we estimate the concentration of butyrate at the bottom of the crypt to be around 1 mM. This concentration gradient is commonly used as a working hypothesis, but it cannot be properly measured in vivo yet. Our model seems however to favour this hypothesis. Moreover, we could estimate the relationship between the butyrate and the O 2 concentrations at the top of the crypt. A microaerobic situation (doubling of the O 2 concentration from homeostasis) was reached at a luminal butyrate concentration of 1 mM.

Second, along with the model, a rigorous mathematical framework has been proposed based on PDMPs, which is in itself a novelty for a crypt model. This frame-work allows for a theoretical analysis of the system. In particular, it ensures the wellposedness of the model, guarantying locally bounded positive solutions. Furthermore, asymptotic analysis arguments lead to the formal derivation of a limit deterministic model of the PDMP. This limit model facilitates and fasten the numerical exploration of the system, by replacing multiple replicates of the same stochastic IBM by one unique deterministic simulation of an approximate average solution.

An extensive numerical exploration of both models (the PDMP model and its deterministic approximation) has been conducted and allowed identifying the main parameters driving different outputs, among them total cell population, progenitor divisions and O 2 distribution. This numerical exploration provides a complete overview of the main model drivers. Interestingly, it suggest that the top/bottom ratio of oxygen in the crypt is sensitive, among others, to cell related parameters, with a complex, nonlinear dependence. This highlights the role of crypt cells in the maintenance of anaerobic conditions in the gut.

Finally the crypt model is coupled with a simplified model of pathogen and microbiota, which demonstrates its potential for crypt-microbiota-pathogen interactions analysis and treatment design.

Improving microbiota integration and mathematical analysis

In this model, the impact of microbe-derived metabolites on the crypt physiology has been added, with in particular the oxygen consumption that drives luminal anaerobia beneficial for the main part of commensals. If the most important mechanisms mediating the interactions between the host and the microbiota have been taken into account, only extremely simplified models of the endogenous microbiota and pathogens were considered, and the microbial populations were not modeled in themselves in this study. In the future, a microbial compartment should be added to accurately model butyrate as well as other metabolites and short chain fatty acids production in an anaerobic environment, as well as inflammatory response and oxygen harmful impact on bacterial populations. This crypt model could also be embedded in macroscopic model of gut physiology, such as for example in [START_REF] Bourgin | Exploring the bacterial impact on cholesterol cycle: A numerical study[END_REF][START_REF] Labarthe | A mathematical model to investigate the key drivers of the biogeography of the colon microbiota[END_REF].

To include the impact of the microbe-derived metabolites we chose to decrease the complexity of other mechanisms such as the mechanical interactions. Some crypt models, for instance [START_REF] Buske | On the biomechanics of stem cell niche formation in the gut -Modelling growing organoids[END_REF], involve a very detailed biophysical models of cell-cell interactions that are crucial for studying crypt morphogenesis or hyperplasia. However, we chose to simplify the description of the mechanics in order to reduce the computational load that was placed instead on the description of butyrate regulations. It also allowed a simpler derivation of the self-aggregation limit model. This simplified version of the mechanical interactions however rendered accurately the main biomarkers of an healthy crypt, in particular migration speed, renewal rate or spatial structure.

Another potential improvement of the present study would be a rigorous derivation of the PME limit model. In this paper, we limited ourselves to a formal derivation, focusing on the main simplifications leading to the limit. To our knowledge, a rigorous derivation of the PME from the self-aggregation model is still an open problem.

All together, this model provides an important building block for modeling hostmicrobiota symbiosis and crypt homeostasis. It could be used to study specifically the interplay between microbial metabolism and host cell turnover, which is crucial for epithelial maturation and microbiota.

A Jump rates definitions

We detail the jump rates used in the model as summed up in Table 1.

A.1 Jump rates for Sect. 4.3 and 4.4

For stem cells:

q ∞ div,sc q div,sc (z, ℓ, d(ν, z)) =1 {ℓ=sc} q ∞ div,sc 1 -R(z, K div,sc [z], κ div,sc [z]) × 1 -R(d(ν, z), K div,sc [dens], κ div,sc [dens]) , (39)
q ∞ sc,pc q sc,pc (z, ℓ) = 1 {ℓ=sc} q ∞ sc,pc R(z, K sc,pc [z], κ sc,pc [z]). (40)
For progenitor cells:

q ∞ div,pc q div,pc (z, ℓ, d(ν, z)) =1 {ℓ=pc} q ∞ div,pc 1 -R(z, K div,pc [z], κ div,pc [z]) × 1 -R(d(ν, z), K div,pc [dens], κ div,pc [dens]) , (41)
q ∞ pc,ent q pc,ent (z, ℓ) = 1 {ℓ=pc} q ∞ pc,ent R(z, K pc,ent [z], κ pc,ent [z]), (42)
q ∞ pc,gc q pc,gc (z, ℓ) = 1 {ℓ=pc} q ∞ pc,gc R(z, K pc,gc [z], κ pc,gc [z]). (43)
For enterocytes:

q ∞ ex,ent q ex,ent (z, ℓ, d(ν, z)) = 1 {ℓ=ent} q ∞ ex,ent R(z, K ex,ent [z], κ ex,ent [z]) × R(d(ν, z), K ex,ent [dens], κ ex,ent [dens]). (44)
For goblet cells:

q ∞ ex,gc q ex,gc (z, ℓ, d(ν, z)) = 1 {ℓ=gc} q ∞ ex,gc R(z, K ex,gc [z], κ ex,gc [z]) × R(d(ν, z), K ex,gc [dens], κ ex,gc [dens]). (45)
A.2 Jump rate definition of the complete model

The jump rate dependency to butyrate, space and cell density is the following. For stem cells:

q ∞ div,sc q div,sc (z, ℓ, d(ν, z)) = 1 {ℓ=sc} q ∞ div,sc 1 -R(z, K div,sc [z], κ div,sc [z]) (46)
× 1 -R(d(ν, z), K div,sc [dens], κ div,sc [dens]) × 1 -R(c b (z), K div,sc [but], κ div,sc [but]) , q ∞ sc,pc q sc,pc (z, ℓ) = 1 {ℓ=sc} q ∞ sc,pc R(z, K sc,pc [z], κ sc,pc [z]). (47)
For progenitor cells:

q ∞ div,pc q div,pc (z, ℓ, d(ν, z)) =1 {ℓ=pc} q ∞ div,pc 1 -R(z, K div,pc [z], κ div,pc [z]) × 1 -R(d(ν, z), K div,pc [dens], κ div,pc [dens]) , (48)
q ∞ pc,ent q pc,ent (z, ℓ) = 1 {ℓ=pc} q ∞ pc,ent R(z, K pc,ent [z], κ pc,ent [z]) × R(c b (z), K pc,ent [but], κ pc,ent [but]), (49)
q ∞ pc,gc q pc,gc (z, ℓ) = 1 {ℓ=pc} q ∞ pc,gc R(z, K pc,gc [z], κ gc,pc [z]) × R(c b (z), K pc,gc [but], κ pc,gc [but]). (50)
For enterocytes:

q ∞ ex,ent q ex,ent (z, ℓ, d(ν, z)) = 1 {ℓ=ent} q ∞ ex,ent R(z, K ex,ent [z], κ ex,ent [z]) × R(ν * D a (z), K ex,ent [dens], κ ex,ent [dens]). (51)
For goblet cells:

q ∞ ex,gc q ex,gc (z, ℓ, d(ν, z)) = 1 {ℓ=gc} q ∞ ex,gc R(z, K ex,gc [z], κ ex,gc [z]) × R(d(ν, z), K ex,gc [dens], κ ex,gc [dens]). (52)

B Existence and regularity of deterministic trajectories

The global existence and C 1 regularity of the ODE part describing cells trajectories comes from the Cauchy-Lipschitz theorem. We briefly sketch existence and regularity results for the PDE part describing concentrations. For the butyrate concentration, the operator

A 1 = σ ∂ zz with domain D(A 1 ) = {c ∈ H 2 (]0, z max [ a ), ∂ z c(-a/2) = 0, c(z max + a/2) = 0}
is a Riesz spectral operator on ]0, z max [ a , that generates a strongly continuous, analytic semigroup S 1 on L 2 (]0, z max [ a ) (see e.g. [START_REF] Curtain | An introduction to infinite dimensional linear systems theory[END_REF]). For the oxygen concentration, the same holds for

A 2 = σ ∂ zz with domain D(A 2 ) = {c ∈ H 2 (]0, z max [ a ), c(-a/2) = 0, ∂ z c(z max + a/2) = 0}, with associated semigroup S 2 .
It fits the framework proposed in [START_REF] Lasiecka | Unified theory for abstract parabolic boundary problems-a semigroup approach[END_REF], so that setting

W 2,1 L 2 (]0,zmax[ a ) = f ∈ L 2 ([0, T ], H 2 (]0, z max [ a )); ∂ f ∂t ∈ L 2 ([0, T ], L 2 (]0, z max [ a )
where H α (]0, z max [ a ) are the standard Sobolev spaces, and denoting S = (S 1 , S 2 ) the semigroup induced on the product space L 2 (]0, z max [ a ) 2 , the following property holds:

Proposition 1 Let c 0 = (c 01 , c 02 ) ∈ H 1 (]0, z max [ a ) 2 , (c lum , c bot ) ∈ H 1 [0, T ] 2 and u = (u 1 , u 2 ) ∈ L 2 ([0, T ], L 2 (]0, z max [ a ) 2 ,
such that c 01 (z max + a/2) = c lum (0) and c 02 (-a/2) = c bot (0), then the evolution problem

∂ t c i -σ ∆ c i = u i , for i=1,2 ∂ z c 1 (-a/2,t) = 0, c 2 (-a/2,t) = c bot , c 1 (z max + a/2,t) = c lum , ∂ z c 2 (z max + a/2,t) = 0, c(•, 0) = c 0 (53)
admits a unique weak solution in W 2,1

L 2 (]0,zmax[ a ) ∩C([0, T ], H 1 (]0, z max [ a )) 2 given by c(z,t) = (1 ]0,zmax[ a (z)c lum (t),1 ]0,zmax[ a (z)c bot (t)) + S(t)(c 0 -(1 ]0,zmax[ a (z)c lum (0), 1 ]0,zmax[ a (z)c bot (0))) - t 0 S(t -s)(1 ]0,zmax[ a (z) dc lum dt (s), 1 ]0,zmax[ a (z) dc bot dt (s))ds + t 0 S(t -s)u(z, s)ds.
Now the semilinear PDEs ( 6) can be viewed as [START_REF] Pearce | Intestinal enteroids recapitulate the effects of short-chain fatty acids on the intestinal epithelium[END_REF] with u replaced by a nonlinear perturbation taking the form of an operator defined on

L 2 (]0, z max [ a ) 2 × [0, T ] by F(c,t) = (-4γ ∞ c 4 o c b c 4 o c b + K β × (ν gc t + ν ent t ) * ψ a (z), -γ ∞ c 4 o c b c 4 o c b + K β × (ν gc t + ν ent t ) * ψ a (z))
This operator is locally Lipschitz in c on L 2 (]0, z max [ a ) 2 , uniformly in t on [0, T ], when ν t is a fixed, measurable function from [0, T ] to M F + (X ) endowed with the weak convergence topology and the associated Borel σ -algebra, such that sup t≤T ⟨ν t , 1⟩ < ∞. Using proposition 1, together with a standard fixed point argument we obtain the local existence (in time) of the solutions of PDEs [START_REF] Barker | Adult intestinal stem cells: Critical drivers of epithelial homeostasis and regeneration[END_REF], that can be extended to a global existence on [0, T ]. More precisely, Proposition 2 Let ν t as above, and c 0 = (c 01 , c 02 )

∈ H 1 (]0, z max [ a ) 2 , (c lum , c bot ) ∈ H 1 [0, T ]
2 such that c 01 (z max + a/2) = c lum (0) and c 02 (-a/2) = c bot (0). The non-linear evolution problem (6) admits a unique

solution c in W 2,1 L 2 (]0,zmax[ a ) ∩C([0, T ], H 1 (]0, z max [ a )) 2 .

C PDMP stability

We aim at showing that an almost surely finite number of jump occurs during any finite time interval. We first introduce the independent Poisson probability measures N k (dθ , i, dt) associated to each jump type

k ∈ E , defined on [0, 1] × N * × R + with intensity n k (dθ , i, dt) = q ∞ k dθ dt ∑ k∈N δ k (i). (54)
The process introduced in Sect. 2 is solution of the following SDE, which is called the pathwise representation of the process

(c t , ν t ) = A t (c 0 , ν 0 ) + ∑ k∈E t 0 N * 1 0 A t-s (c s , ν s -+ µ k (x i s -))) -A t-s (c s , ν s -) 1 {i≤⟨ν s -,1⟩} (i)1 {θ ≤q k (x i s -,νs- * D k (x i s -),cs * ψa(z i s -))} (θ )N k (dθ , i, ds). ( 55
)
Fig. 18 sketches the rationale behind this formula, built upon the remark that if T n and T n-1 are successive jump times, then for t > T n ,

A t-T n-1 (c T n-1 , ν T n-1 ) = A t-Tn (c Tn , ν T - n ).
Note that, when possible without ambiguity, we simplify the notation 1 {θ ≤q k (x i s -,νs- * D k (x i s -),cs * ψa(z i s -))} (θ ) by 1 {θ ≤q k (•)} (θ ).

For all n ∈ N, we define the stopping time τ n = inf{t ≥ 0, ⟨ν t , 1⟩ = n}. Let T n (t) := min(t, τ n ).

(56)

The stopped process at n cells is (c Tn(t) , ν Tn(t) ) t≥0 . C p,q,• (R + × X , R) denotes the space of functions f (t, z, ℓ) that are p times derivable with respect to t and q times derivable with respect to z for all ℓ ∈ T .

Lemma C1 Let f ∈ C 1,1,• (R + × X , R), g ∈ H 1 (]0, z max [ a ) 2 , ψ ∈ C 1 (R 2 ), t ∈ R.
We denote by ⟨c, g⟩ L 2 the scalar product in L 2 . Then, for t > 0 and n ≥ ⟨ν 0 , 1⟩, ψ(⟨c Tn(t) , g⟩ L 2 , ⟨ν Tn(t) , f Tn(t) ⟩) = ψ(⟨c 0 , g⟩, ⟨ν 0 , f 0 ⟩)

+ t 0 ∂ 1 ψ(⟨c s , g⟩ L 2 , ⟨ν s , f s ⟩)⟨∂ s c s , g⟩ L 2 1 {s≤Tn(t)} (s) ds + t 0 ∂ 2 ψ(⟨c s , g⟩ L 2 , ⟨ν s , f s ⟩)⟨ν s , φ ν s * F a ∇ z f s + ∂ s f s ⟩1 {s≤Tn(t)} (s) ds + ∑ k∈E t 0 N * 1 0 ψ(⟨c s , g⟩ L 2 , ⟨ν s -+ µ k (x i s -), f s ⟩) -ψ(⟨c s , g⟩ L 2 , ⟨ν s -, f s ⟩) × 1 {s≤Tn(t)} 1 {i≤⟨ν s -,1⟩} 1 {θ ≤q k (•)} N k (dθ , di, ds). ( 57
)
Proof Using the pathwise representation of the PDMP (55), we get

⟨ν Tn(t) , f Tn(t) ⟩ = ⟨ ÃTn(t) (c 0 , ν 0 ), f Tn(t) ⟩ + ∑ k∈E Tn(t) 0 N * 1 0 ⟨ ÃTn(t)-s (c s , ν s -+ µ k (x i s -)), f Tn(t) ⟩ -⟨ ÃTn(t)-s (c s , ν s -), f Tn(t) ⟩ × 1 {i≤⟨ν s -,1⟩} 1 {θ ≤q k (•)} N k (dθ , di, ds)
where the stochastic integral is correctly defined since the population size is majored by n, and f is continuous on the compact X . Then, applying formula [START_REF] Buchwald | A review of human salmonellosis: Ii. duration of excretion following infection with nontyphi salmonella[END_REF] , the first term of the RHS can be written

⟨ ÃTn(t) (c 0 , ν 0 ), f Tn(t) ⟩ = ⟨ν 0 , f 0 ⟩ + Tn(t) 0 ⟨ Ãu (c 0 , ν 0 ), φ Ãu (c 0 , ν 0 ) * F a ∇ z f u + ∂ u f u ⟩du,
and the second one

∑ k∈E Tn(t) 0 N * 1 0 ⟨ν s -+ µ k (x i s -), f s ⟩ -⟨ν s -, f s ⟩ × 1 {i≤⟨ν s -,1⟩} 1 {θ ≤q k (•)} N k (dθ , di, ds) + ∑ k∈E Tn(t) 0 N * 1 0 Tn(t) s ⟨ Ãu-s (c s , ν s -+ µ k (x i s -)), φ Ãu-s (c s , ν s -+ µ k (x i s -)) * F a ∇ z f u + ∂ u f u ⟩ -⟨ Ãu-s (c s , ν s -), φ Ãu-s (c s , ν s -) * F a ∇ z f u + ∂ u f u ⟩ × 1 {i≤⟨ν s -,1⟩}
1 {θ ≤q k (•)} (θ ) N k (dθ , di, ds) (59)
≤ ⟨ν 0 , 1⟩ p + ∑ k∈E Tn(t)

0 N * 1 0 [(⟨ν s -, 1⟩ + 1) p -⟨ν s -, 1⟩ p ] 1 {i≤⟨ν s -,1⟩} (i)
1 {θ ≤q k (•)} (θ ) N k (dθ , di, ds). (60)
The last integral grows with time since the integrand is non-negative, so that the supremum can be applied in the left hand side. Furthermore, for any non-negative x, (1 + x) px p ≤ C(p)(1 + x p-1 ). Hence, sup

s≤Tn(t) ⟨ν s , 1⟩ p ≤ ⟨ν 0 , 1⟩ p + ∑ k∈E Tn(t) 0 N * 1 0 C(p) ⟨ν s -, 1⟩ p-1 + 1 1 {i≤⟨ν s -,1⟩} (i)
1 {θ ≤q k (•)} (θ )N k (dθ , di, ds). (61)
Switching to the expectations:

E sup s≤Tn(t) ⟨ν s , 1⟩ p ≤ E[⟨ν 0 , 1⟩ p ] + ∑ k∈E E   Tn(t) 0 N * 1 0 C(p) ⟨ν s -, 1⟩ p-1 + 1 1 {i≤⟨ν s -,1⟩} (i)1 {θ ≤q k (•)} (θ )q ∞ k dθ ∑ k≥1 δ k (di) ds   ≤ E[⟨ν 0 , 1⟩ p ] + ∑ k∈E E   Tn(t) 0 ⟨ν s -,1⟩ ∑ i=1 C(p) ⟨ν s -, 1⟩ p-1 + 1 q ∞ k q k (•) ds   ≤ E[⟨ν 0 , 1⟩ p ] + ∑ k∈E E   t 0 q ∞ k C(p) (⟨ν s∧τn , 1⟩ p + ⟨ν s∧τn , 1⟩) ds   .
Recalling that n + n p ≤ 2n p for all p ≥ 1, Fubini's theorem leads to

E sup s≤Tn(t) ⟨ν s , 1⟩ p ≤ E[⟨ν 0 , 1⟩ p ] + 2C(p) ∑ k∈E q ∞ k t 0 E sup u≤s∧τn ⟨ν u , 1⟩ p ds.
and Gronwall lemma provides an estimate that does not depend on n:

E sup s≤Tn(t) ⟨ν s , 1⟩ p ≤ E[⟨ν 0 , 1⟩ p ] exp 2C(p) ∑ k∈E q ∞ k t = C(p,t). (62)
Hence τ n → ∞ a.s. Indeed, suppose that there exists T > 0 such that, for all n ∈ N, P(τ n ≤ T ) > ε > 0.

We remark that

τ n ≤ T ⇒ sup s≤Tn(T )
⟨ν s , 1⟩ p ≥ n p so that P( sup The stability of the PDMP is a direct consequence of this lemma.

s≤Tn(T ) ⟨ν s , 1⟩ p ≥ n p ) ≥ P(τ n ≤ T ) Hence,
Corollary C3 If there exists an integer p ≥ 1 such that E[⟨ν 0 , 1⟩ p ] < ∞, then:

1. The PDMP is stable. 2. T n (t) can be replaced by t in formula [START_REF] Rogan | Bacterial Manipulation of Wnt Signaling: A Host-Pathogen Tug-of-Wnt[END_REF].

Proof (1): Let (J n ) n be the sequence of jump times. We aim at showing that J n → n→∞ ∞ a.s. Assume that there exists T < +∞ and a set Ω T of non-null probability such that, ∀ ω ∈ Ω T , sup n∈N J n (ω) < T . We recall that, as proven in the proof of Lemma C2, that E sup s≤t ⟨ν s , 1⟩ < ∞ implies that τ n → ∞ a.s. Hence, it exists N ∈ N and a set Ω T,N ⊂ Ω T of non-null probability such that, ∀ ω ∈ Ω T,N , T < τ N (ω). Let C(N) be an upper bound of the jump rate on M P + (X ) N = {ν ∈ M P + (X ), ⟨ν, 1⟩ ≤ N} ⊂ M P + (X ). Such an upper bound exists since the individual jump rates are bounded by q ∞ k by definition. Hence, for all ω ∈ Ω T,N , the jump time sequence (J n (ω)) n∈N can be built as an extraction of the sequence (J ′ n (ω)) n∈N of the jump time of a Poisson process with jump rate C(N). But the sequence (J ′ n (ω)) n∈N is a.s. divergent, which comes in contradiction with sup n∈N J n (ω) < T for all ω ∈ Ω T .

(2): As the stochastic integrals in the proof of Lemma C1 are correctly defined (see e.g. theorem A.3 [START_REF] Bansaye | Stochastic models for structured populations: scaling limits and long time behavior[END_REF]) due to the hypothesis on ν 0 and Lemma C2.

D A case of bistability

Here we explain the unexpected simulations values we obtained for K div [dens] = -50% in Sect. 4.2. Starting from an initial number of 700 cells, we see in Fig. 19a that either the number of cells drops below 400 cells and activity in the crypt goes on, or the number of cells stabilises around 500 and the crypt freezes, meaning that no additional extrusion or division events are observed, so that the total cell population does not evolve any longer. We evaluated the fraction of crypt that freezes upon variations of K div [dens], always starting with 700 cells (Fig. 19b). This fraction increases when K div [dens] decreases. Indeed, for low values of K div [dens], cell division is totally inhibited by contact inhibition and the division rate is 0. At first, cells are extruded at the top of the crypt as density of cells (starting from 700 cells) is way above what K div [dens] prescribes. Once exceeding cells have been extruded, the population can reach a mechanical equilibrium where no extrusion occurs at the top of the crypt (as density equilibrium has been reached). Meanwhile, no division occurs in the lower part of the crypt as density is still too high there, and the division rate is therefore zero.

E Parameters values for the different models

PDMP model parameters are given for the different sub-models of increasing complexity studied during the PDMP exploration. Namely, the parameters used for the model with the sole deterministic part, studied in Sect. 4.1, are given in Table 7. When the spatialized division and extrusion events are added in Sect. 10. Finally, the complete model including the feedback of butyrate on stem cell division as studied in Sect. 4.5 is implemented with the parameters presented in Table 11.

4.2, parameters of

Parameters used for the deterministic PDE approximation of the PDMP are given in Table 12.

E.1 Parameters of chemical kinetic

We could not find value for the molar concentration of oxygen. Therefore, we switched to normalized concentrations. We set co = co C with C a normalizing constant. The evolution equation of oxygen concentration 6 becomes

∂ t co -σ i ∂ zz co = - s o C γ ∞ β co 4 c b co 4 c b + K 5 β C 4 (ν ent t + ν gc t ) * ψ a (z). (63)
We suppose that experimental measures of the reaction speed in [START_REF] Clausen | Kinetic studies on the metabolism of short-chain fatty acids and glucose by isolated rat colonocytes[END_REF] are done at a standard (normoxic) oxygen concentration c * o . The half of the maximal reaction speed is reached experimentally at a butyrate concentration of 0.184 mole/L. Therefore,

K 5 β = 0.184 × c * 4 o .
We suppose as well that oxygen concentration at the bottom of the crypt is c

  a) Model features and expected outcome. (b) Crypt geometry.(c) Effect of the curvature on mechanical interactions between cells.

Fig. 1 :

 1 Fig. 1: Summary of the model. (a) First, the differentiation path between stem cells, progenitor cells and differentiated cells (enterocytes and goblet cells) is displayed on the left. Next, the spatial distribution of crypt cell populations at homeostasis (expected model outcome) is represented, from the stem cell niche at the bottom of the crypt, to the differentiated cells at the top. Finally, regulation mechanisms (mechanical, metabolic and population dynamics mechanisms) driving the spatial structure in the model are listed, including stochastic events (cell division, differentiation and extrusion), cell migration, cell density, enterocyte metabolism and metabolite gradients (butyrate, O 2 and W nt). (b) We also indicate the geometrical assumptions underlying the definition of f (z) including the profiles of the crypt domain and of the function φ defined in Eq. (4) modeling curvature effects in the mechanical model. (c) Note that φ (z 3 )F a (z * ) < φ (z 1 )F a (z * ) due to the curvature of the basal membrane, modeled by the phenomenological term φ (•).

Fig. 2 :

 2 Fig. 2: Model function profiles. Profiles of interaction kernels and regulation functions taken into account in the model

Fig. 3 :

 3 Fig. 3: Trajectory of cells due to repulsive mechanical interactions only. Initial positions of 700 cells are drawn in a law U(0, z max ), with same initial condition for both simulations.

Fig. 4 :

 4 Fig. 4: Cell migration. For an initial condition of 700 cells uniformly distributed along the crypt, we indicate the cells trajectories (a) and the spatial distribution of stochastic events (b) for one simulation, and the number of cells for 16 randomly selected simulations, as well as the mean and standard deviation for 50 repetitions of the PDMP (c). Division events are visible in (a) with comb-like patterns in the lower zone of the crypt corresponding to daughter cell that joins a close layer. As expected, division and extrusion are spatially segregated (b). A stationary number of cells is reached after about 200 hours, indicating an equilibrium between division and extrusion (c).

Fig. 5 :

 5 Fig. 5: Local sensitivity analysis of hallmark of crypt physiology to the driving parameters. Three parameters that potentialy impact the total number of cells at final time, the migration speed and the number of divisions in the crypt at steady state are tested. Parameters are shifted one by one, and model output variations are plotted. Mean (plain) and standard deviation (dashed lines) from 50 repetitions are plotted.

2 .

 2 min = 40 µm, which corresponds to a healthy crypt. The corresponding low steady-state value of oxygen at the top of the crypt is c o,min = 1 au. The threshold c * b and c * o for butyrate and oxygen respectively were then defined as the average of the infected and healthy state: c * b = c b,min + c b,max 2 and c * o = c o,min + c o,max The value for µ = µ b k b

10 . 5 β× 10 4 .

 1054 * o . The exact value of c * o is unknown to us, but we set c * o = 10. Then, C = Therefore, we can find a value for the parameter K 5 β /C 4 of (63): K As C is unknown, the parameter s o /C is still unknown, but we decided to use C = 1 in all simulations.

( a )Fig. 6 :Fig. 7 :Fig. 8 :Fig. 9 :Fig. 10 :

 a678910 Fig. 6: Crypt structure. The structure of the crypt by cell type is studied, by displaying the time evolution of the cell populations, and their spatial distribution at final time.

Fig. 12 :

 12 Fig.12: Comparison of densities and cell fate events between PDMP and PME model at steady-state. Results obtained using parameters from tables 12 for the PME and 11 for the PDMP (in particular, the cell size for the PDMP is a = 10 µm). For the PDMP, spatial distribution of cells (resp. cell fate events) are histograms of local cell (resp. jumps) counts with mean and standard deviation from the mean from 50 simulations.

Fig. 13 :Fig. 14 :

 1314 Fig. 13: Sensitivity analysis results and ranking indices Figures on the left represent the standard deviation of EEs against the mean of absolute values of EEs. Parameters are listed in order of influence (circle: increasing effect ; square: decreasing effects).On the right we represent the criteria convergence with a threshold represented by the dashed lines. We recall that convergences are obtained if the ranking index is below 1 and screening below 0.05 based on the work of ([START_REF] Sarrazin | Global Sensitivity Analysis of environmental models: Convergence and validation[END_REF]). In order to make the convergences visibly comparable, the ranking index is multiplied by 0.05 to share the same threshold as the screening index. We recall that K ext,gc [z] = K ext,ent [z] and K ext,gc [dens] = K ext,ent[dens]

( a )

 a Low luminal butyrate. (b) Increased Wnt signaling.

Fig. 15 :

 15 Fig. 15: Two pathogenic steady states: (a) low luminal butyrate c b,lum = 0.5 mM (b) increased Wnt signaling K div,pc = 160 µm. On the left, cell densities (PDMP as histogram, and PME model as blue line). Densities of the PME model have been rescaled by the bin size of the histogram for comparison with the PDMP. On the right, butyrate and O 2 concentrations along the crypt.

Fig. 16 :

 16 Fig. 16: Comparison of the impact of different treatment on the crypt. x-axis labels refer to treatment: "2 * 1" means supplementation of 2 mM of butyrate for 1 h every 10 h. Values are taken at final time.

Fig. 17 : 2 ( 1 ) 1 + 2 ) 2 +Fig. 18 : 1 +Fig. 19 :

 172112218119 Fig. 17: Simulated treatment of an infected crypt with butyrate supplementation. Starting from a pathogenic situation (c b,endo (0) = c b,min = 0.5 mM, K div,pc [z](0) = K div,pc [z] max = 160 µm), the butyrate concentration c b,treat takes the form of a constant supplementation of 2 mM for 7.75 hours every ten hours.

Table 1 :

 1 •, K k [ j], κ k [ j])for Enumeration The values of K k [ j] and κ k [ j] can be found in Table11. Detailed formulations can be found in Appendix A.2.

	Regulation pathway ( j)

Table 2 :

 2 Regulation of cell fate events used in Sect. 4.3. Note that this table is a subset of Table 1 as the model used in Sect. 4.3 is a submodel of the complete model presented in Sect. 2.

	Observable	Simulated value	Reference value Reference
	total number of cells	694 ± 21 (at t f )	700	[51]
	number of stem cells	23 ± 7 (at t f )	20	[61]
	ratio of progenitors/total	0.22 ± 0.01 (at t f )	0.22 (human)	[11]
	ratio goblets/enterocytes	0.33 ± 0.03 (at t f )	0.33	[70]
	95 % renewal time (day)	3.6 ± 0.02	2.67 ± 2.04	[25]
	division rate of stem cells (h -1 )	1/83 ± 1/83	1/24(SI)-1/62	[7, 70]
	division rates of progenitors (h -1 )	1/18 ± 0	1/12(SI)-1/62	[8, 70]
	migration speed (µm/h)	2.6 ± 0.03	4	[51]

Table 4 :

 4 Simulations results and experimental reference for the complete IBM (Sect. 4.5), mean of 50 simulations, t f = 500h, SI = small intestine. Note that in this table, two rows are added with respect to Table3: one related to butyrate and the other to oxygen.

Table 5 :

 5 .2 confirmed that this was hardly achievable without modifying other parameter values in the model. Deterministic approximation model outputs. See Eq. (

	Observable	Simulated value
	Total number of cells	709
	Number of stem cells	20
	Ratio progenitors/total	0.14
	Ratio goblets/enterocytes	0.3
	Division rate of stem cells (h -1 )	0.041
	Differentiation rate of stem cells (h -1 )	0.096
	Division rate of progenitor cells (h -1 )	0.041
	Differentiation rate of progenitor cells in goblet cells (h -1 )	0.015
	Differentiation rate of progenitor cells in enterocytes (h -1 )	0.046
	Extrusion rate of goblet cells (h -1 )	0.007
	Extrusion rate of enterocytes (h -1 )	0.007
	Ratio of O 2 concentration at bottom/at top	10.02
	Butyrate concentration at bottom (mM)	1.79

  1 {θ ≤q k (•)} duN k (dθ , di, ds).applied to ⟨ν Tn(t) , 1⟩ p gives ⟨ν Tn(t) , 1⟩ p =⟨ν 0 , 1⟩ p -+ µ k (x i s -), 1⟩ p -⟨ν s -, 1⟩ p 1 {i≤⟨ν s -,1⟩} (i)

		Tn(t)	1
	k∈E + ∑	0 N *	0	⟨ν s

  by Markov inequality, for all n ∈ N 1⟩ p ≥ n p n p ≥ P(τ n ≤ T ) n p > εn p

	E sup Hence, lim n→∞ sup s≤Tn(t) ⟨ν s , 1⟩ p = sup ⟨ν s , 1⟩ p ≥ P sup s≤t ⟨ν s , 1⟩ p a.s. Fatou's lemma finally gives ⟨ν s , which is in contradiction with (62). E sup s≤t ⟨ν s , 1⟩ p = E lim inf n→∞ sup s≤Tn(t) ⟨ν s , 1⟩ p ≤ lim inf s≤Tn(t) n→∞ E sup ⟨ν s , 1⟩ p ≤ C(p,t).
	s≤Tn(T )	s≤Tn(T )

  Table 8 are used. In Sect. 4.3 where cell types are introduced, parameters of Table 9 are used. Adding O 2 and butyrate in Sect. 4.4, we use parameters presented in Table

Table 7 :

 7 Parameters for Sect. 4.1.

		Value	Unit	Description	Ref
	z max	200	µm	height of the crypt	[51]
	r 0	25	µm	shape parameter of the crypt	-
	ε	= 0.1 × r 0	µm	shape parameter of the crypt	-
	a	10	µm	cell diameter	[14]
	k a	15	µm × h -1	maximal intensity of F a	-
	d a	1	/ 0	maximal value of D a	-

Table 8 :

 8 Parameters for Sect. 4.2.

		Value	Unit	Description	Ref
	z max	200	µm	height of the crypt	[51]
	r 0	25	µm	shape parameter of the crypt	-
	ε	= 0.1 × r 0	µm	shape parameter of the crypt	-
	a	10	µm	cell diameter	[14]
	k a	15	µm × h -1	maximal intensity of F a	-
	d a	1	/ 0	maximal value of D a	-
	q ∞ div	0.22	h -1	maximal rate of division	-
	K div [z]	40	µm	space regulation parameter for division	-
	κ div [z]	15	µm	space regulation parameter for division	-
	K div [dens]	40	cell	cell density regulation parameter for division	-
	κ div [dens]	6	cell	cell density regulation parameter for division	-
	q ∞ ex	0.34	h -1	maximal rate of extrusion	-
	K ex [z]	190	µm	space regulation parameter for extrusion	-
	κ ex [z]	15	µm	space regulation parameter for extrusion	-

Table 9 :

 9 Parameters for Sect. 4.3.

		Value	Unit	Description	Ref
	z max	200	µm	height of the crypt	[51]
	r 0	25	µm	shape parameter of the crypt	-
	ε	= 0.1 × r 0	µm	shape parameter of the crypt	-
	a	10	µm	cell diameter	[14]
	k a	15	µm × h -1	maximal intensity of F a	-
	d a	1	/ 0	maximal value of D a	-
	N dcs	12	cell	number of DCS cells	[61]
	ρ dcs (z)dz	see (26)		spatial distribution of DCS cells	-
	q ∞ div,sc	0.09	h -1	maximal rate of stem cell division	-
	K div,sc [z]	12	µm	space regulation parameter for stem cell div.	-
	κ div,sc [z]	5	µm	space regulation parameter for stem cell div.	-
	K div,sc [dens]	53	cell	cell density regulation parameter for stem cell div.	-
	κ div,sc [dens]	6	cell	cell density regulation parameter for stem cell div.	-
	q ∞ sc,pc	0.2	h -1	max. rate of stem cell to progenitor diff.	-
	K sc,pc [z]	= K div,sc [z]	µm	space regulation parameter for stem cell to progenitor diff.	-
	κ sc,pc [z]	= κ div,sc [z]	µm	space regulation parameter for stem cell to progenitor differentiation	-
	q ∞ div,pc	0.22	h -1	max. rate of progenitor division	-
	K div,pc [z]	40	µm	space regulation parameter for progenitor division	-
	κ div,pc [z]	40	µm	space regulation parameter for progenitor division	-
	K div,pc [dens]	41	cell	cell density regulation parameter for progenitor division	-
	κ div,pc [dens]	= κ div,sc [dens]	cell	cell density regulation parameter for progenitor division	-
	q ∞ pc,ent	0.15	h -1	maximal rate of progenitor to enterocyte differentiation	-
	K pc,ent [z]	= K div,pc [z]	µm	space regulation parameter for progenitor to enterocyte diff.	-
	κ pc,ent [z]	15	µm	space regulation parameter for progenitor to enterocyte diff.	-
	q ∞ pc,gc	= 0.33 × q ∞ pc,ent	h -1	maximal rate of progenitor to goblet differentiation	-
	K pc,gc [z]	= K pc,ent [z]	µm	space regulation parameter for progenitor to goblet diff.	-
	κ pc,gc [z]	= κ pc,ent [z]	µm	space regulation parameter for progenitor to goblet diff.	-
	q ∞ ex,ent	0.34	h -1	maximal rate of enterocyte extrusion	-
	K ex,ent [z]	190	µm	space regulation parameter for enterocyte extrusion	-
	κ ex,ent [z]	15	µm	space regulation parameter for enterocyte extrusion	-
	K ex,ent [dens]	20	cell	cell density regulation parameter for enterocyte extrusion	-
	κ ex,ent [dens]	= κ div,sc [dens]	cell	cell density regulation parameter for enterocyte extrusion	-
	q ∞ ex,gc	= q ∞ ex,ent	h -1	maximal rate of goblet cell extrusion	-
	K ex,gc [z]	= K ex,ent [z]	µm	space regulation parameter for goblet extrusion	-
	κ ex,gc [z]	= κ ex,ent [z]	µm	space regulation parameter for goblet extrusion	-
	K ex,gc [dens]	= K ex,ent [dens]	cell	cell density regulation parameter for goblet extrusion	-
	κ ex,gc [dens]	= κ ex,ent [dens]	cell	cell density regulation parameter for goblet extrusion	-

Table 10 :

 10 Parameters for Sect. 4.4.

		Value	Unit	Description	Ref
	z max	200	µm	height of the crypt	[51]
	r 0	25	µm	shape parameter of the crypt	-
	ε	= 0.1 × r 0	µm	shape parameter of the crypt	-
	a	10	µm	cell diameter	[14]
	k a	15	µm × h -1	maximal intensity of F a	-
	d a	1	/ 0	maximal value of D a	-
	N dcs	12	cell	number of DCS cells	[61]
	ρ dcs (z)dz	see (26)		spatial distribution of DCS cells	-
	q ∞ div,sc	0.09	h -1	maximal rate of stem cell division	-
	K div,sc [z]	12	µm	space regulation parameter for stem cell div.	-
	κ div,sc [z]	5	µm	space regulation parameter for stem cell div.	-
	K div,sc [dens]	53	cell	cell density regulation parameter for stem cell div.	-
	κ div,sc [dens]	6	cell	cell density regulation parameter for stem cell div.	-
	q ∞ sc,pc	0.2	h -1	max. rate of stem cell to progenitor diff.	-
	K sc,pc [z]	= K div,sc [z]	µm	space regulation parameter for stem cell to progenitor diff.	-
	κ sc,pc [z]	= κ div,sc [z]	µm	space regulation parameter for stem cell to progenitor differentiation	-
	q ∞ div,pc	0.22	h -1	max. rate of progenitor division	-
	K div,pc [z]	40	µm	space regulation parameter for progenitor division	-
	κ div,pc [z]	40	µm	space regulation parameter for progenitor division	-
	K div,pc [dens]	41	cell	cell density regulation parameter for progenitor division	-
	κ div,pc [dens]	= κ div,sc [dens]	cell	cell density regulation parameter for progenitor division	-
	q ∞ pc,ent	0.15	h -1	maximal rate of progenitor to enterocyte differentiation	-
	K pc,ent [z]	= K div,pc [z]	µm	space regulation parameter for progenitor to enterocyte diff.	-
	κ pc,ent [z]	15	µm	space regulation parameter for progenitor to enterocyte diff.	-
	q ∞ pc,gc	= 0.33 × q ∞ pc,ent	h -1	maximal rate of progenitor to goblet differentiation	-
	K pc,gc [z]	= K pc,ent [z]	µm	space regulation parameter for progenitor to goblet diff.	-
	κ pc,gc [z]	= κ pc,ent [z]	µm	space regulation parameter for progenitor to goblet diff.	-
	q ∞ ex,ent	0.34	h -1	maximal rate of enterocyte extrusion	-
	K ex,ent [z]	190	µm	space regulation parameter for enterocyte extrusion	-
	κ ex,ent [z]	15	µm	space regulation parameter for enterocyte extrusion	-
	K ex,ent [dens]	20	cell	cell density regulation parameter for enterocyte extrusion	-
	κ ex,ent [dens]	= κ div,sc [dens]	cell	cell density regulation parameter for enterocyte extrusion	-
	q ∞ ex,gc	= q ∞ ex,ent	h -1	maximal rate of goblet cell extrusion	-
	K ex,gc [z]	= K ex,ent [z]	µm	space regulation parameter for goblet extrusion	-
	κ ex,gc [z]	= κ ex,ent [z]	µm	space regulation parameter for goblet extrusion	-
	K ex,gc [dens]	= K ex,ent [dens]	cell	cell density regulation parameter for goblet extrusion	-
	κ ex,gc [dens]	= κ ex,ent [dens]	cell	cell density regulation parameter for goblet extrusion	-
	s o	-4	/ 0	stoechiometric coefficient for O 2	-
	s b	-1	/ 0	stoechiometric coefficient for butyrate	-
	σ o	3.6 × 10 6	µm 2 /h	diffusion coefficient for O 2	[44]
	σ but	= σ o	µm 2 /h	diffusion coefficient for butyrate	-
	γ ∞ β	1.6 × 10 4	10 -15 × mmol h×µm 2	max. reaction speed of β -oxydation	[18]
	K β	(0.184 × 10 4 ) 1/5	mM	affinity of β -oxydation (supposed experimental normoxic condition) [19]
	c b,lum	5	mM	luminal concentration of butyrate at homeostasis	[33]
	c o,bot	10	a.u.	Homestatic O 2 concentration at the bottom of the crypt	[56]

Table 11 :

 11 Parameters for Sect. 4.5.

		Value	Unit	Description	Ref
	z max	200	µm	height of the crypt	[51]
	r 0	25	µm	shape parameter of the crypt	-
	ε	= 0.1 × r 0	µm	shape parameter of the crypt	-
	a	10	µm	cell diameter	[14]
	k a	15	µm × h -1	maximal intensity of F a	-
	d a	1	/ 0	maximal value of D a	-
	N dcs	12	cell	number of DCS cells	[61]
	ρ dcs (z)dz	see (26)		spatial distribution of DCS cells	-
	q ∞ div,sc	0.15	h -1	maximal rate of stem cell division	-
	K div,sc [z]	12	µm	space regulation parameter for stem cell div.	-
	κ div,sc [z]	5	µm	space regulation parameter for stem cell div.	-
	K div,sc [dens]	53	cell	cell density regulation parameter for stem cell div.	-
	κ div,sc [dens]	6	cell	space regulation parameter for stem cell div.	-
	K div,sc [but]	2		butyrate regulation parameter for stem cell div.	-
	κ div,sc [but]	5		butyrate regulation parameter for stem cell div.	-
	q ∞ sc,pc	0.2	h -1	max. rate of stem cell to progenitor diff.	-
	K sc,pc [z]	= K div,sc [z]	µm	space regulation parameter of stem cell to progenitor diff.	-
	κ sc,pc [z]	= κ div,sc [z]	µm	space regulation parameter of stem cell to progenitor diff.	-
	q ∞ div,pc	0.22	h -1	maximal rate of progenitor division	-
	K div,pc [z]	40	µm	space regulation parameter for progenitor division	-
	κ div,pc [z]	40	µm	space regulation parameter for progenitor division	-
	K div,pc [dens]	41	cell	cell density regulation parameter for progenitor division	-
	κ div,pc [dens]	= κ div,sc [dens]	cell	cell density regulation parameter for progenitor division	-
	q ∞ pc,ent	0.25	h -1	maximal rate of progenitor to enterocyte diff.	-
	K pc,ent [z]	= K div,pc [z]	µm	space regulation parameter for progenitor to enterocyte diff.	-
	κ pc,ent [z]	15	µm	space regulation parameter for progenitor to enterocyte diff.	-
	K pc,ent [but]	1.5		butyrate regulation parameter for progenitor to enterocyte diff.	-
	κ pc,ent [but]	5		butyrate regulation parameter for progenitor to enterocyte diff.	-
	q ∞ pc,gc	= 0.33 × q ∞ pc,ent	h -1	maximal rate of progenitor to goblet differentiation	-
	K pc,gc [z]	= K pc,ent [z]	µm	space regulation parameter for progenitor to goblet diff.	-
	κ pc,gc [z]	= κ pc,ent [z]	µm	space regulation parameter for progenitor to goblet diff.	-
	K pc,gc [but]	K pc,ent [but]		butyrate regulation parameter for progenitor to goblet diff.	-
	κ pc,gc [but]	κ pc,ent [but]		butyrate regulation parameter for progenitor to goblet diff.	-
	q ∞ ex,gc	0.34	h -1	maximal rate of goblet cell extrusion	-
	K ex,gc [z]	190	µm	space regulation parameter for goblet extrusion	-
	κ ex,gc [z]	15	µm	space regulation parameter for goblet extrusion	-
	K ex,gc [dens]	20	cell	cell density regulation parameter for goblet extrusion	-
	κ ex,gc [dens]	κ div,sc [dens]	cell	cell density regulation parameter for goblet extrusion	-
	q ∞ ex,ent	= q ∞ ex,gc	h -1	maximal rate of enterocyte extrusion	-
	K ex,ent [z]	K ex,gc [z]	µm	space regulation parameter for enterocyte extrusion	-
	κ ex,ent [z]	κ ex,gc [z]	µm	space regulation parameter for enterocyte extrusion	-
	K ex,ent [dens]	= K ex,gc [dens]	cell	cell density regulation parameter for enterocyte extrusion	-
	κ ex,ent [dens]	= κ ex,gc [dens]	cell	space reuglation parameter for enterocyte extrusion	-
	s o	-4	/ 0	stoechiometric coefficient for O 2	-
	s b	-1	/ 0	stoechiometric coefficient for butyrate	-
	σ o	3.6 × 10 6	µm 2 /h	diffusion coefficient for O 2	[44]
	σ but	= σ o	µm 2 /h	diffusion coefficient for butyrate	-
	γ ∞ β	1.6 × 10 4	10 -15 × mmol h×µm 2	max. reaction speed of β -oxydation	[18]
	K β	(0.184 × 10 4 ) 1/5	mM	affinity of β -oxydation	[19]
	c b,lum	5	mM	luminal concentration of butyrate at homeostasis	[33]
	c o,bot	10	a.u.	homestatic O 2 concentration at the bottom of the crypt	[56]

Marie Haghebaert's contribution was permitted in part by funding from the European Research Council (ERC) under the European Union's Horizon 2020 research and innovation programme (Grant agreement ERC-2017-AdG No. 788191 -Homo.symbiosus)

Due to the smoothness of f , ∂ t f and ∇ z f , Fubini's theorem and formula [START_REF] Purcell | Life at low Reynolds number[END_REF] are applied to the last integral term to get:

We then have:

We also have ⟨c Tn(t) , g⟩ L 2 = ⟨c 0 , g⟩ + t 0 ⟨∂ s c s , g⟩ L 2 1 {s≤Tn(t)} ds.

Then, multidimensional Itô's formula with jumps is applied to get the desired expression.

We now have a tool to study the stability of the PDMP. We must show that the population does not explode in finite time.

where C(p,t) < ∞ is a constant that only depends on p and t.

Proof For all n ∈ N, let T n (t) be as defined in [START_REF] Rivera-Chávez | Oxygen as a driver of gut dysbiosis[END_REF]. As the deterministic evolution does not impact the population size and the stochastic events modify it with at most one individual by jump, formula (57)