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

Abstract

Human health and physiology is strongly influenced by interactions between human cells and intestinal microbiota in the gut. In mammals, the host-microbiota crosstalk is mainly mediated by regulations at the intestinal crypt level: the epithelial cell turnover in crypts is directly influenced by metabolites produced by the microbiota. Conversely, enterocytes maintain hypoxia in the gut, favorable to anaerobic bacteria which dominate the gut microbiota. We constructed an individual-based model of epithelial cells interacting with the microbiota-derived chemicals diffusing in the crypt lumen. This model is formalized as a piecewise deterministic Markov process (PDMP). It accounts for local interactions due to cell contact (among which are mechanical interactions), for cell proliferation, differentiation and extrusion which are regulated spatially or by chemicals concentrations. It also includes chemicals diffusing and reacting with cells. A deterministic approximated model is also introduced for a large population of small cells, expressed as a system of porous media type equations. Both models are extensively studied through numerical exploration. Their biological relevance is thoroughly assessed by recovering bio-markers of an healthy crypt, such as cell population distribution along the crypt or population turn-over rates. Simulation results from the deterministic model are compared to the PMDP model and we take advantage of its lower computational cost to perform a sensitivity analysis by Morris method. We finally use the crypt model to explore butyrate supplementation to enhance recovery after infections by enteric pathogens.

Appendices

Jump rates definitions

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

1.1 Jump rates for Sects. 4.3 and 4.4

For stem cells:

$$\begin{aligned} q^{\infty }_{\text {div,sc}}q_{\text {div,sc}}(z,\ell ,{\mathbf {d}}(\nu ,z))&= \mathbbm {1}_{\{ \ell =\text {sc} \}}q^{\infty }_{\text {div,sc}} \Big ( 1 - R(z,K_{\text {div,sc}}[z],\kappa _{\text {div,sc}}[z]) \Big ) \nonumber \\&\quad \times \Big ( 1 - R({\mathbf {d}}(\nu ,z),K_{\text {div,sc}}[dens],\kappa _{\text {div,sc}}[dens]) \Big ), \end{aligned}$$

(39)

$$\begin{aligned} q^{\infty }_{\text {sc},\text {pc}} q_{\text {sc},\text {pc}}(z,\ell )&= \mathbbm {1}_{\{ \ell =\text {sc} \}}q^{\infty }_{\text {sc},\text {pc}} R(z,K_{\text {sc},\text {pc}}[z],\kappa _{\text {sc},\text {pc}}[z]). \end{aligned}$$

(40)

For progenitor cells:

$$\begin{aligned} q^{\infty }_{\text {div},\text {pc}}q_{\text {div},\text {pc}}(z,\ell ,{\mathbf {d}}(\nu ,z))&= \mathbbm {1}_{\{ \ell =\text {pc} \}}q^{\infty }_{\text {div},\text {pc}} \Big ( 1 - R(z,K_{\text {div},\text {pc}}[z],\kappa _{\text {div},\text {pc}}[z]) \Big )\nonumber \\&\quad \times \Big ( 1 - R({\mathbf {d}}(\nu ,z),K_{\text {div},\text {pc}}[dens],\kappa _{\text {div},\text {pc}}[dens]) \Big ), \end{aligned}$$

(41)

$$\begin{aligned} q^{\infty }_{\text {pc},\text {ent}}q_{\text {pc},\text {ent}}(z,\ell )&= \mathbbm {1}_{\{ \ell =\text {pc} \}}q^{\infty }_{\text {pc},\text {ent}} R(z,K_{\text {pc},\text {ent}}[z],\kappa _{\text {pc},\text {ent}}[z]), \end{aligned}$$

(42)

$$\begin{aligned} q^{\infty }_{\text {pc},\text {gc}} q_{\text {pc},\text {gc}}(z,\ell )&= \mathbbm {1}_{\{ \ell =\text {pc} \}}q^{\infty }_{\text {pc},\text {gc}} R(z,K_{\text {pc},\text {gc}}[z],\kappa _{\text {pc},\text {gc}}[z]). \end{aligned}$$

(43)

For enterocytes:

$$\begin{aligned} q^{\infty }_{\text {ex},\text {ent}} q_{\text {ex},\text {ent}}(z,\ell ,{\mathbf {d}}(\nu ,z))= & {} \mathbbm {1}_{\{ \ell =\text {ent} \}}q^{\infty }_{\text {ex},\text {ent}} R(z,K_{\text {ex},\text {ent}}[z],\kappa _{\text {ex},\text {ent}}[z])\nonumber \\&\times R({\mathbf {d}}(\nu ,z),K_{\text {ex},\text {ent}}[dens],\kappa _{\text {ex},\text {ent}}[dens]). \end{aligned}$$

(44)

For goblet cells:

$$\begin{aligned} q^{\infty }_{\text {ex},\text {gc}}q_{\text {ex},\text {gc}}(z,\ell ,{\mathbf {d}}(\nu ,z))= & {} \mathbbm {1}_{\{ \ell =\text {gc} \}}q^{\infty }_{\text {ex},\text {gc}} R(z,K_{\text {ex},\text {gc}}[z],\kappa _{\text {ex},\text {gc}}[z])\nonumber \\&\times R({\mathbf {d}}(\nu ,z),K_{\text {ex},\text {gc}}[dens],\kappa _{\text {ex},\text {gc}}[dens]). \end{aligned}$$

(45)

1.2 Jump rate definition of the complete model

The jump rate dependency to butyrate, space and cell density is the following.

For stem cells:

$$\begin{aligned} q^{\infty }_{\text {div},\text {sc}} q_{\text {div,sc}}(z,\ell ,{\mathbf {d}}(\nu ,z))&= \mathbbm {1}_{\{ \ell =\text {sc} \}}q^{\infty }_{\text {div,sc}} \Big ( 1 - R(z,K_{\text {div,sc}}[z],\kappa _{\text {div,sc}}[z]) \Big ) \nonumber \\&\quad \times \Big ( 1 - R({\mathbf {d}}(\nu ,z),K_{\text {div,sc}}[dens],\kappa _{\text {div,sc}}[dens]) \Big )\nonumber \\&\quad \times \Big (1 - R({\mathbf {c}}_{{{\mathbf {b}}}}(z),K_{\text {div},\text {sc}}[but],\kappa _{\text {div},\text {sc}}[but])\Big ), \end{aligned}$$

(46)

$$\begin{aligned} q^{\infty }_{\text {sc},\text {pc}} q_{\text {sc},\text {pc}}(z,\ell )&= \mathbbm {1}_{\{ \ell =\text {sc} \}}q^{\infty }_{\text {sc},\text {pc}} R(z,K_{\text {sc},\text {pc}}[z],\kappa _{\text {sc},\text {pc}}[z]). \end{aligned}$$

(47)

For progenitor cells:

$$\begin{aligned} q^{\infty }_{\text {div},\text {pc}} q_{\text {div},\text {pc}}(z,\ell ,{\mathbf {d}}(\nu ,z))&= \mathbbm {1}_{\{ \ell =\text {pc} \}}q^{\infty }_{\text {div},\text {pc}} \Big ( 1 - R(z,K_{\text {div},\text {pc}}[z],\kappa _{\text {div},\text {pc}}[z]) \Big ) \nonumber \\&\quad \times \Big ( 1 - R({\mathbf {d}}(\nu ,z),K_{\text {div},\text {pc}}[dens],\kappa _{\text {div},\text {pc}}[dens]) \Big ), \end{aligned}$$

(48)

$$\begin{aligned} q^{\infty }_{\text {pc},\text {ent}} q_{\text {pc},\text {ent}}(z,\ell )&= \mathbbm {1}_{\{ \ell =\text {pc} \}}q^{\infty }_{\text {pc},\text {ent}} R(z,K_{\text {pc},\text {ent}}[z],\kappa _{\text {pc},\text {ent}}[z]) \nonumber \\&\quad \times R({\mathbf {c}}_{{{\mathbf {b}}}}(z),K_{\text {pc},\text {ent}}[but],\kappa _{\text {pc},\text {ent}}[but]), \end{aligned}$$

(49)

$$\begin{aligned} q^{\infty }_{\text {pc},\text {gc}} q_{\text {pc},\text {gc}}(z,\ell )&= \mathbbm {1}_{\{ \ell =\text {pc} \}}q^{\infty }_{\text {pc},\text {gc}} R(z,K_{\text {pc},\text {gc}}[z],\kappa _{\text {gc},\text {pc}}[z]) \nonumber \\&\quad \times R({\mathbf {c}}_{{{\mathbf {b}}}}(z),K_{\text {pc},\text {gc}}[but],\kappa _{\text {pc},\text {gc}}[but]). \end{aligned}$$

(50)

For enterocytes:

$$\begin{aligned} q^{\infty }_{\text {ex},\text {ent}} q_{\text {ex},\text {ent}}(z,\ell ,{\mathbf {d}}(\nu ,z))= & {} \mathbbm {1}_{\{ \ell =\text {ent} \}}q^{\infty }_{\text {ex},\text {ent}} R(z,K_{\text {ex},\text {ent}}[z],\kappa _{\text {ex},\text {ent}}[z]) \nonumber \\&\times R(\nu * D_a(z),K_{\text {ex},\text {ent}}[dens],\kappa _{\text {ex},\text {ent}}[dens]). \end{aligned}$$

(51)

For goblet cells:

$$\begin{aligned} q^{\infty }_{\text {ex},\text {gc}} q_{\text {ex},\text {gc}}(z,\ell ,{\mathbf {d}}(\nu ,z))= & {} \mathbbm {1}_{\{ \ell =\text {gc} \}}q^{\infty }_{\text {ex},\text {gc}} R(z,K_{\text {ex},\text {gc}}[z],\kappa _{\text {ex},\text {gc}}[z])\nonumber \\&\times R({\mathbf {d}}(\nu ,z),K_{\text {ex},\text {gc}}[dens],\kappa _{\text {ex},\text {gc}}[dens]). \end{aligned}$$

(52)

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=\sigma \partial _{zz}$ with domain ${\mathscr {D}}(A_1)=\{c \in H^2(\left]0,z_{\text {max}}\right[_a),\, \partial _z c(-a/2)=0,\,c(z_{\text {max}}+a/2)=0\}$ is a Riesz spectral operator on $\left]0,z_{\text {max}}\right[_a$, that generates a strongly continuous, analytic semigroup $S_1$ on $L^2(\left]0,z_{\text {max}}\right[_a)$ (see e.g. Curtain and Zwart 1995). For the oxygen concentration, the same holds for $A_2=\sigma \partial _{zz}$ with domain ${\mathscr {D}}(A_2)=\{c \in H^2(\left]0,z_{\text {max}}\right[_a),\,c(-a/2)=0,\, \partial _z c(z_{\text {max}}+a/2)=0 \}$, with associated semigroup $S_2$. It fits the framework proposed in Lasiecka (1980), so that setting

$$\begin{aligned} W^{2,1}_{L^2(\left]0,z_{\text {max}}\right[_a)} = \left\{ f \in L^2([0,T],H^{2}(\left]0,z_{\text {max}}\right[_a)); \frac{\partial f}{\partial t} \in L^2([0,T],L^2(\left]0,z_{\text {max}}\right[_a) \right\} \end{aligned}$$

where $H^{\alpha }(\left]0,z_{\text {max}}\right[_a)$ are the standard Sobolev spaces, and denoting $S=(S_1,S_2)$ the semigroup induced on the product space $L^2(\left]0,z_{\text {max}}\right[_a)^2$, the following property holds:

Proposition 1

Let $c_0=(c_{01},c_{02}) \in H^1(\left]0,z_{\text {max}}\right[_a)^2$, $(c_{\text {lum}},c_{\text {bot}}) \in \left( H^1[0,T]\right) ^2$ and $u=(u_1,u_2) \in \left( L^2([0,T],L^2(\left]0,z_{\text {max}}\right[_a)\right) ^2$, such that $c_{01}(z_{\text {max}}+a/2)=c_{\text {lum}}(0)$ and $c_{02}(-a/2)=c_{\text {bot}}(0)$, then the evolution problem

$$\begin{aligned}&\partial _t c_i - \sigma \varDelta c_i = u_i,\text { for i=1,2}\nonumber \\&\partial _z c_1(-a/2,t) = 0,\quad c_2(-a/2,t) = c_{bot},\nonumber \\&c_1(z_{\text {max}}+a/2,t)=c_{lum},\quad \partial _z c_2(z_{\text {max}}+a/2,t)=0,\nonumber \\&c(\cdot ,0) = c_0 \end{aligned}$$

(53)

admits a unique weak solution in $\left( W^{2,1}_{L^2(\left]0,z_{\text {max}}\right[_a)} \cap C([0,T],H^1(\left]0,z_{\text {max}}\right[_a))\right) ^2$ given by

$$\begin{aligned} c(z,t)&= (\mathbbm {1}_{\left]0,z_{\text {max}}\right[_a}(z) c_{\text {lum}}(t), \mathbbm {1}_{\left]0,z_{\text {max}}\right[_a}(z)c_{\text {bot}}(t)) \\&\quad +S(t)(c_0 - (\mathbbm {1}_{\left]0,z_{\text {max}}\right[_a}(z)c_{\text {lum}}(0),\mathbbm {1}_{\left]0,z_{\text {max}}\right[_a}(z)c_{\text {bot}}(0))) \\&\quad - \int _0^t S(t-s)(\mathbbm {1}_{\left]0,z_{\text {max}}\right[_a}(z)\frac{dc_{\text {lum}}}{dt}(s),\mathbbm {1}_{\left]0,z_{\text {max}}\right[_a}(z)\frac{dc_{\text {bot}}}{dt}(s)) ds \\&\quad + \int _0^t S(t-s) u(z,s)ds. \end{aligned}$$

Now the semilinear PDEs (6) can be viewed as (53) with u replaced by a nonlinear perturbation taking the form of an operator defined on $L^2(\left]0,z_{\text {max}}\right[_a)^2 \times [0,T]$ by

$$\begin{aligned} F(c,t)= & {} (-4\gamma ^\infty \frac{c_o^4 c_b}{c_o^4 c_b + K_{\beta }} \times (\nu _t^{\text {gc}} + \nu _t^{\text {ent}})*\psi _{a}(z),\\&-\gamma ^\infty \frac{c_o^4 c_b}{c_o^4 c_b + K_{\beta }}\times (\nu _t^{\text {gc}} + \nu _t^{\text {ent}})*\psi _{a}(z)) \end{aligned}$$

This operator is locally Lipschitz in c on $L^2(\left]0,z_{\text {max}}\right[_a)^2$, uniformly in t on [0, T], when $\nu _t$ is a fixed, measurable function from [0, T] to ${\mathscr {M}}_{F^+}({\mathscr {X}})$ endowed with the weak convergence topology and the associated Borel $\sigma$-algebra, such that $\sup _{t \le T} \langle \nu _t,1\rangle < \infty$. Using proposition 1, together with a standard fixed point argument we obtain the local existence (in time) of the solutions of PDEs (6), that can be extended to a global existence on [0, T]. More precisely,

Proposition 2

Let $\nu _t$ as above, and $c_0=(c_{01},c_{02}) \in H^1(\left]0,z_{\text {max}}\right[_a)^2$, $(c_{\text {lum}},c_{\text {bot}}) \in \left( H^1[0,T]\right) ^2$ such that $c_{01}(z_{\text {max}}+a/2)=c_{\text {lum}}(0)$ and $c_{02}(-a/2)=c_{\text {bot}}(0)$. The non-linear evolution problem (6) admits a unique solution c in $\Big (W^{2,1}_{L^2(\left]0,z_{\text {max}}\right[_a)} \cap C([0,T],H^1(\left]0,z_{\text {max}}\right[_a)) \Big )^2$.

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\theta ,i,dt)$ associated to each jump type $k \in {\mathscr {E}}$, defined on $[0,1]\times {\mathbb {N}}^* \times {\mathbb {R}}_+$ with intensity

$$\begin{aligned} n_k(d\theta ,i,dt) = q^{\infty }_k d\theta dt \sum \limits _{k \in {\mathbb {N}}} \delta _k(i). \end{aligned}$$

(54)

The process introduced in Sect. 2 is solution of the following SDE, which is called the pathwise representation of the process

$$\begin{aligned} (c_{t},\nu _{t})= & {} A_{t}(c_0,\nu _0) + \sum _{k \in {\mathscr {E}}} \int \limits _0^t \int \limits _{{\mathbb {N}}^*} \int \limits _0^1 \left[ A_{t-s}(c_s,\nu _{s^-} + \mu _k(x^i_{s^-}))) - A_{t-s}(c_s,\nu _{s^-}) \right] \nonumber \\&\mathbbm {1}_{\{ i \le \langle \nu _{s^-},1 \rangle \}}(i) \mathbbm {1}_{\{ \theta \le q_k(x^i_{s^-},\nu _{s-}*D_k (x^i_{s^-}),c_s*\psi _{a}(z^i_{s^-})) \}}(\theta ) N_k(d\theta ,i,ds). \end{aligned}$$

(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$,

$$\begin{aligned} A_{t - T_{n-1}}(c_{T_{n-1}},\nu _{T_{n-1}}) = A_{t - T_n}(c_{T_n},\nu _{T_n^-}). \end{aligned}$$

Note that, when possible without ambiguity, we simplify the notation $\mathbbm {1}_{\{ \theta \le q_k(x^i_{s^-},\nu _{s-}*D_k (x^i_{s^-}),c_s*\psi _{a}(z^i_{s^-})) \}}(\theta )$ by $\mathbbm {1}_{\{ \theta \le q_k(\cdot ) \}}(\theta )$.

Fig. 18

Sketch of the process trajectory. The trajectory $A_t(c_0,\nu _0)$ is first plotted. After the first jump, the trajectory $A_t(c_0,\nu _0)$ for $t \ge T_1$, which equals $A_{t - T_1}(c_{T_1},\nu _{T_1^{-}})$, is replaced by $A_{t - T_1}(c_{T_1},\nu _{T_1^{-}}+\mu _k(T_1^{-}))$. This process is iterated at each jump

For all $n \in {\mathbb {N}}$, we define the stopping time $\tau _n = \inf \{ t \ge 0, \langle \nu _t ,1 \rangle = n \}$. Let

$$\begin{aligned} T_n(t):= \min (t,\tau _n). \end{aligned}$$

(56)

The stopped process at n cells is $(c_{T_n(t)},\nu _{T_n(t)})_{t \ge 0}$. ${\mathscr {C}}^{p,q,\bullet }({\mathbb {R}}^+ \times {\mathscr {X}},{\mathbb {R}})$ denotes the space of functions $f(t,z,\ell )$ that are p times derivable with respect to t and q times derivable with respect to z for all $\ell \in {\mathscr {T}}$.

Lemma C1

Let $f \in {\mathscr {C}}^{1,1,\bullet }({\mathbb {R}}^+ \times {\mathscr {X}},{\mathbb {R}})$, $g \in H^1(\left]0,z_{\text {max}}\right[_a)^{2}$ , $\psi \in {\mathscr {C}}^1({\mathbb {R}}^2)$, $t \in {\mathbb {R}}$. We denote by $\langle c,g \rangle _{L^2}$ the scalar product in $L^2$. Then, for $t>0$ and $n \ge \langle \nu _0,1 \rangle$,

$$\begin{aligned}&\psi (\langle c_{T_n(t)},g \rangle _{L^2},\langle \nu _{T_n(t)},f_{T_n(t)} \rangle ) \nonumber \\&\quad = \psi (\langle c_0,g \rangle ,\langle \nu _0,f_0 \rangle ) \nonumber \\&\qquad + \int \limits _0^{t} \partial _1 \psi (\langle c_s,g \rangle _{L^2},\langle \nu _s,f_s \rangle ) \langle \partial _s c_s, g \rangle _{L^2} \mathbbm {1}_{\{ s \le T_n(t) \}}(s)\ ds \nonumber \\&\quad \quad + \int \limits _0^{t} \partial _2 \psi (\langle c_s,g \rangle _{L^2},\langle \nu _s,f_s \rangle ) \langle \nu _s, \phi \nu _s * F_a\nabla _z f_s + \partial _s f_s \rangle \mathbbm {1}_{\{ s \le T_n(t) \}}(s)\ ds \nonumber \\&\quad \quad + \sum _{k \in {\mathscr {E}}} \int \limits _0^{t} \int \limits _{{\mathbb {N}}^*} \int \limits _0^1 \left[ \psi (\langle c_s,g \rangle _{L^2},\langle \nu _{s^-} + \mu _k(x^i_{s^-}), f_s \rangle ) - \psi (\langle c_s,g \rangle _{L^2},\langle \nu _{s^-}, f_s \rangle ) \right] \nonumber \\&\quad \quad \times \,\mathbbm {1}_{\{ s \le T_n(t) \}}\mathbbm {1}_{\{ i \le \langle \nu _{s^-},1 \rangle \}} \mathbbm {1}_{\{ \theta \le q_k(\cdot ) \}} N_k(d\theta ,di,ds). \end{aligned}$$

(57)

Proof

Using the pathwise representation of the PDMP (55), we get

$$\begin{aligned}&\langle \nu _{T_n(t)},f_{T_n(t)} \rangle \\&\quad = \langle {\tilde{A}}_{T_n(t)}(c_0,\nu _0),f_{T_n(t)} \rangle \\&\quad \quad + \sum _{k \in {\mathscr {E}}} \int \limits _0^{T_n(t)} \int \limits _{{\mathbb {N}}^*} \int \limits _0^1 \left[ \langle {\tilde{A}}_{T_n(t)-s}(c_s,\nu _{s^-} + \mu _k(x^i_{s^-})),f_{T_n(t)} \rangle - \langle {\tilde{A}}_{T_n(t)-s}(c_s,\nu _{s^-}),f_{T_n(t)} \rangle \right] \\&\quad \quad \times \mathbbm {1}_{\{ i \le \langle \nu _{s^-},1 \rangle \}} \mathbbm {1}_{\{ \theta \le q_k(\cdot ) \}} N_k(d\theta ,di,ds) \end{aligned}$$

where the stochastic integral is correctly defined since the population size is majored by n, and f is continuous on the compact ${\mathscr {X}}$.

Then, applying formula (12), the first term of the RHS can be written

$$\begin{aligned}&\langle {\tilde{A}}_{T_n(t)}(c_0,\nu _0),f_{T_n(t)} \rangle = \langle \nu _0,f_0 \rangle \\&\qquad + \int \limits _0^{T_n(t)} \langle {\tilde{A}}_u(c_0,\nu _0), \phi {\tilde{A}}_u(c_0,\nu _0) * F_a\nabla _z f_u + \partial _u f_u \rangle du, \end{aligned}$$

and the second one

$$\begin{aligned}&\sum \limits _{k \in {\mathscr {E}}}\int \limits _0^{T_n(t)} \int \limits _{{\mathbb {N}}^*} \int \limits _0^1 \left( \langle \nu _{s^-} + \mu _k(x^i_{s^-}) ,f_s \rangle - \langle \nu _{s^-} , f_s \rangle \right) \\&\qquad \times \mathbbm {1}_{\{ i \le \langle \nu _{s^-}, 1 \rangle \}} \mathbbm {1}_{\{ \theta \le q_k(\cdot ) \}} N_k(d\theta ,di,ds) \\&\qquad + \sum \limits _{k \in {\mathscr {E}}}\int \limits _0^{T_n(t)} \int \limits _{{\mathbb {N}}^*} \int \limits _0^1 \int \limits _s^{T_n(t)} \left[ \langle {\tilde{A}}_{u-s}(c_s,\nu _{s^-} + \mu _k(x^i_{s^-})), \phi {\tilde{A}}_{u-s}(c_s,\nu _{s^-} + \mu _k(x^i_{s^-})) \right. \\&\qquad * F_a\nabla _z f_u + \partial _u f_u \rangle - \langle {\tilde{A}}_{u-s}(c_s,\nu _{s^-}), \phi {\tilde{A}}_{u-s}(c_s,\nu _{s^-}) * F_a\nabla _z f_u + \partial _u f_u \rangle \big ] \\&\qquad \times \,\mathbbm {1}_{\{ i \le \langle \nu _{s^-}, 1 \rangle \}} \mathbbm {1}_{\{ \theta \le q_k(\cdot ) \}}\ du N_k(d\theta ,di,ds). \end{aligned}$$

Due to the smoothness of f, $\partial _t f$ and $\nabla _z f$, Fubini’s theorem and formula (55) are applied to the last integral term to get:

$$\begin{aligned}&\int \limits _0^{T_n(t)} \sum \limits _{k \in {\mathscr {E}}} \int \limits _0^u \int \limits _{{\mathbb {N}}^*} \int \limits _0^1 \left[ \langle {\tilde{A}}_{u-s}(c_s,\nu _{s^-} + \mu _k(x^i_{s^-})), \phi {\tilde{A}}_{u-s}(c_s,\nu _{s^-} + \mu _k(x^i_{s^-})) * F_a\nabla _z f_u + \partial _u f_u \rangle \right. \\&\qquad - \langle {\tilde{A}}_{u-s}(c_s,\nu _{s^-}), \phi {\tilde{A}}_{u-s}(c_s,\nu _{s^-}) * F_a\nabla _z f_u + \partial _u f_u \rangle \Big ] \\&\qquad \times \mathbbm {1}_{\{ i \le \langle \nu _{s^-}, 1 \rangle \}} \mathbbm {1}_{\{ \theta \le q_k(\cdot ) / q^{\infty }_k \}}\ N_k(d\theta ,di,ds)du \\&\quad = \int \limits _0^{T_n(t)} \langle \nu _u,\phi \nu _u * F_a\nabla _z f_u + \partial _u f_u \rangle - \langle {\tilde{A}}_u(c_0,\nu _0), \phi {\tilde{A}}_u(c_0,\nu _0) * F_a\nabla _z f_u + \partial _u f_u \rangle \ du. \end{aligned}$$

We then have:

$$\begin{aligned}&\langle \nu _{T_n(t)},f_{T_n(t)} \rangle \\&\quad = \langle \nu _0,f_0 \rangle + \int \limits _0^{T_n(t)} \langle {\tilde{A}}_u(c_0,\nu _0), \phi {\tilde{A}}_u(c_0,\nu _0) * F_a\nabla _z f_u + \partial _u f_u \rangle du \\&\qquad + \int \limits _0^{T_n(t)} \langle \nu _u,\phi \nu _u * F_a\nabla _z f_u + \partial _u f_u \rangle \\&\qquad - \int \limits _0^{T_n(t)} \langle {\tilde{A}}_u(c_0,\nu _0), \phi {\tilde{A}}_u(c_0,\nu _0) * F_a\nabla _z f_u + \partial _u f_u \rangle \ du \\&\qquad +\sum \limits _{k \in {\mathscr {E}}}\int \limits _0^{T_n(t)} \int \limits _{{\mathbb {N}}^*} \int \limits _0^1\left( \langle \nu _{s^-} + \mu _k(x^i_{s^-}) ,f_s \rangle - \langle \nu _{s^-} , f_s \rangle \right) \\&\qquad \mathbbm {1}_{\{ i \le \langle \nu _{s^-}, 1 \rangle \}} \mathbbm {1}_{\{ \theta \le q_k(\cdot ) \}}\ N_k(d\theta ,di,ds) \\&\quad = \langle \nu _0,f_0 \rangle + \int \limits _0^t \langle \nu _u,\phi \nu _u * F_a\nabla _z f_u + \partial _u f_u \rangle \mathbbm {1}_{\{ u \le T_n(t) \}} \ du \\&\qquad + \sum \limits _{k \in {\mathscr {E}}}\int \limits _0^t \int \limits _{{\mathbb {N}}^*} \int \limits _0^1 \langle \mu _k(x^i_{s^-}) ,f_s \rangle \mathbbm {1}_{\{ s \le T_n(t) \}} \mathbbm {1}_{\{ i \le \langle \nu _{s^-}, 1 \rangle \}} \mathbbm {1}_{\{ \theta \le q_k(\cdot ) \}}\ N_k(d\theta ,di,ds). \end{aligned}$$

We also have

$$\begin{aligned} \langle c_{T_n(t)},g \rangle _{L^2} = \langle c_0,g \rangle + \int _{0}^{t} \langle \partial _s c_s,g \rangle _{L^2}\mathbbm {1}_{\{ s \le T_n(t) \}}\ ds. \end{aligned}$$

Then, multidimensional Itô’s formula with jumps is applied to get the desired expression. $\square$

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

Lemma C2

For any $t \ge 0$, if an integer $p \ge 1$ such that ${\mathbb {E}} \left[ \langle \nu _0,1 \rangle ^p \right] < \infty$ exists then

$$\begin{aligned} {\mathbb {E}} \left[ \underset{s \le t}{\sup } \langle \nu _s,1 \rangle ^p \right] \le C(p,t) \end{aligned}$$

(58)

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

Proof

For all $n \in {\mathbb {N}}$, let $T_n(t)$ be as defined in (56). 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) applied to $\langle \nu _{T_n(t)},1 \rangle ^p$ gives

$$\begin{aligned} \langle \nu _{T_n(t)},1 \rangle ^p&= \langle \nu _0,1 \rangle ^p \nonumber \\&\quad + \sum \limits _{k \in {\mathscr {E}}} \int \limits _0^{T_n(t)} \int \limits _{{\mathbb {N}}^*} \int \limits _0^1 \left( \langle \nu _{s^-} + \mu _k(x^i_{s^-}) ,1 \rangle ^p - \langle \nu _{s^-} ,1 \rangle ^p \right) \mathbbm {1}_{\{ i \le \langle \nu _{s^-},1\rangle \}}(i) \nonumber \\&\quad \quad \quad \quad \quad \mathbbm {1}_{\{ \theta \le q_k(\cdot ) \}}(\theta )\ N_k(d\theta ,di,ds) \end{aligned}$$

(59)

$$\begin{aligned}&\le \langle \nu _0,1 \rangle ^p + \sum \limits _{k \in {\mathscr {E}}} \int \limits _0^{T_n(t)} \int \limits _{{\mathbb {N}}^*} \int \limits _0^1 \left[ \left( \langle \nu _{s^-} ,1 \rangle + 1 \right) ^p - \langle \nu _{s^-} ,1 \rangle ^p \right] \mathbbm {1}_{\{ i \le \langle \nu _{s^-},1\rangle \}}(i) \nonumber \\&\quad \quad \quad \quad \quad \mathbbm {1}_{\{ \theta \le q_k(\cdot ) \}}(\theta )\ N_k(d\theta ,di,ds). \end{aligned}$$

(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)^p -x^p \le C(p) (1 + x^{p-1})$. Hence,

$$\begin{aligned} \underset{s \le T_n(t)}{\sup }\langle \nu _s,1 \rangle ^p \le \langle \nu _0,1 \rangle ^p +&\sum \limits _{k \in {\mathscr {E}}} \int \limits _0^{T_n(t)} \int \limits _{{\mathbb {N}}^*} \int \limits _0^1 C(p) \left( \langle \nu _{s^-} ,1 \rangle ^{p-1} + 1 \right) \mathbbm {1}_{\{ i \le \langle \nu _{s^-},1\rangle \}}(i) \nonumber \\&\quad \quad \quad \quad \quad \mathbbm {1}_{\{ \theta \le q_k(\cdot ) \}}(\theta ) N_k(d\theta ,di,ds). \end{aligned}$$

(61)

Switching to the expectations:

$$\begin{aligned}&{\mathbb {E}} \left[ \underset{s \le T_n(t)}{\sup }\langle \nu _s,1 \rangle ^p \right] \le {\mathbb {E}} \left[ \langle \nu _0,1 \rangle ^p \right] \\&\qquad + \sum \limits _{k \in {\mathscr {E}}} {\mathbb {E}} \left[ \int \limits _0^{T_n(t)} \int \limits _{{\mathbb {N}}^*} \int \limits _0^1 C(p) \left( \langle \nu _{s^-} ,1 \rangle ^{p-1} + 1 \right) \mathbbm {1}_{\{ i \le \langle \nu _{s^-},1\rangle \}}(i) \mathbbm {1}_{\{ \theta \le q_k(\cdot ) \}}(\theta ) q^{\infty }_k d\theta \left( \sum \limits _{k \ge 1}\delta _k(di) \right) ds \right] \\&\quad \le {\mathbb {E}} \left[ \langle \nu _0,1 \rangle ^p \right] + \sum \limits _{k \in {\mathscr {E}}} {\mathbb {E}} \left[ \int \limits _0^{T_n(t)} \sum \limits _{i=1}^{\langle \nu _{s^-},1 \rangle }C(p) \left( \langle \nu _{s^-} ,1 \rangle ^{p-1} + 1 \right) q^{\infty }_k q_k(\cdot )\ ds \right] \\&\quad \le {\mathbb {E}} \left[ \langle \nu _0,1 \rangle ^p \right] + \sum \limits _{k \in {\mathscr {E}}} {\mathbb {E}} \left[ \int \limits _0^t q^{\infty }_k C(p) \left( \langle \nu _{s \wedge \tau _n} ,1 \rangle ^p + \langle \nu _{s \wedge \tau _n} ,1 \rangle \right) ds \right] . \end{aligned}$$

Recalling that $n + n^p \le 2n^p$ for all $p\ge 1$, Fubini’s theorem leads to

$$\begin{aligned}&{\mathbb {E}} \left[ \underset{s \le T_n(t)}{\sup }\langle \nu _s,1 \rangle ^p \right] \le {\mathbb {E}} \left[ \langle \nu _0,1 \rangle ^p \right] + 2 C(p) \left( \sum \limits _{k \in {\mathscr {E}}} q^{\infty }_k \right) \int \limits _0^t {\mathbb {E}} \left[ \underset{u \le s \wedge \tau _n }{\sup } \langle \nu _u ,1 \rangle ^p \right] ds. \end{aligned}$$

and Gronwall lemma provides an estimate that does not depend on n:

$$\begin{aligned} {\mathbb {E}} \left[ \underset{s \le T_n(t)}{\sup }\langle \nu _s,1 \rangle ^p \right] \le {\mathbb {E}} \left[ \langle \nu _0,1 \rangle ^p \right] \exp \left( 2 C(p) \left( \sum \limits _{k \in {\mathscr {E}}} q^{\infty }_k \right) t\right) = C(p,t). \end{aligned}$$

(62)

Hence $\tau _n \rightarrow \infty$ a.s. Indeed, suppose that there exists $T>0$ such that, for all $n\in {\mathbb {N}}$, ${\mathbb {P}}\left( \tau _n \le T\right)> \varepsilon > 0$. We remark that

$$\begin{aligned} \tau _n \le T \Rightarrow \underset{s \le T_n(T)}{\sup }\langle \nu _s,1 \rangle ^p \ge n^p \quad \text { so that } \quad {\mathbb {P}}(\underset{s \le T_n(T)}{\sup }\langle \nu _s,1 \rangle ^p \ge n^p ) \ge {\mathbb {P}}(\tau _n \le T) \end{aligned}$$

Hence, by Markov inequality, for all $n\in {\mathbb {N}}$

$$\begin{aligned} {\mathbb {E}} \left[ \underset{s \le T_n(T)}{\sup }\langle \nu _s,1 \rangle ^p \right] \ge {\mathbb {P}}\left( \underset{s \le T_n(T)}{\sup }\langle \nu _s,1 \rangle ^p \ge n^p\right) n^p \ge {\mathbb {P}}\left( \tau _n \le T \right) n^p > \varepsilon n^p \end{aligned}$$

which is in contradiction with (62).

Hence, $\underset{n \rightarrow \infty }{\lim }\ \underset{s \le T_n(t)}{\sup }\langle \nu _s,1 \rangle ^p = \underset{s \le t }{\sup }\langle \nu _s,1 \rangle ^p$ a.s. Fatou’s lemma finally gives

$$\begin{aligned} {\mathbb {E}} \left[ \underset{s \le t }{\sup }\langle \nu _s,1 \rangle ^p \right]= & {} {\mathbb {E}} \left[ \underset{n \rightarrow \infty }{\liminf } \underset{s \le T_n(t)}{\sup }\langle \nu _s,1 \rangle ^p \right] \\&\le \underset{n \rightarrow \infty }{\liminf }\ {\mathbb {E}} \left[ \underset{s \le T_n(t)}{\sup }\langle \nu _s,1 \rangle ^p \right] \le C(p,t). \square \end{aligned}$$

The stability of the PDMP is a direct consequence of this lemma.

Corollary C3

If there exists an integer $p \ge 1$ such that ${\mathbb {E}} \left[ \langle \nu _0,1 \rangle ^p \right] < \infty$, then:

• 1. The PDMP is stable.

• 2. $T_n(t)$ can be replaced by t in formula (57).

Proof

(1): Let $(J_n)_n$ be the sequence of jump times. We aim at showing that $J_n \underset{n \rightarrow \infty }{\rightarrow } \infty$ a.s. Assume that there exists $T < +\infty$ and a set $\varOmega _T$ of non-null probability such that, $\forall \ \omega \in \varOmega _T$, $\sup \limits _{n \in {\mathbb {N}}} J_n(\omega ) < T$. We recall that, as proven in the proof of Lemma C2, that ${\mathbb {E}} \left[ \sup \limits _{s \le t} \langle \nu _s,1 \rangle \right] < \infty$ implies that $\tau _n \rightarrow \infty$ a.s. Hence, it exists $N \in {\mathbb {N}}$ and a set $\varOmega _{T,N}\subset \varOmega _T$ of non-null probability such that, $\forall \ \omega \in \varOmega _{T,N}$, $T < \tau _N(\omega )$. Let $C(N)$ be an upper bound of the jump rate on ${\mathscr {M}}_{P^+}({\mathscr {X}})_N = \left\{ \nu \in {\mathscr {M}}_{P^+}({\mathscr {X}}), \langle \nu ,1 \rangle \le N \right\} \subset {\mathscr {M}}_{P^+}({\mathscr {X}})$. Such an upper bound exists since the individual jump rates are bounded by $q_k^\infty$ by definition. Hence, for all $\omega \in \varOmega _{T,N}$, the jump time sequence $(J_n(\omega ))_{n \in {\mathbb {N}}}$ can be built as an extraction of the sequence $(J'_n(\omega ))_{n \in {\mathbb {N}}}$ of the jump time of a Poisson process with jump rate $C(N)$. But the sequence $(J'_n(\omega ))_{n \in {\mathbb {N}}}$ is a.s. divergent, which comes in contradiction with $\sup \limits _{n \in {\mathbb {N}}} J_n(\omega ) < T$ for all $\omega \in \varOmega _T$.

(2): As the stochastic integrals in the proof of Lemma C1 are correctly defined (see e.g. theorem A.3 Bansaye and Méléard 2015) due to the hypothesis on $\nu _0$ and Lemma C2. $\square$

A case of bistability

Here we explain the unexpected simulations values we obtained for $K_{\text {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_{\text {div}}[dens]$, always starting with 700 cells (Fig. 19b). This fraction increases when $K_{\text {div}}[dens]$ decreases. Indeed, for low values of $K_{\text {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_{\text {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.

Fig. 19

Illustration of bistability phenomena for high $K_{\text {div}}[dens]$. We plot the evolution of the total number of cells in 50 repetitions of the model in Sect. 4.2 for $K_{\text {div}}[dens]=20.5$ (left panel). We observe that some simulations show a constant total cell population (around 500 cells) for $t>250$ without additional stochastic variations, whereas other simulations show stochastic evolution without interruption, but a significantly lower total cell population (around 300 cells). Then, we compute 50 repetitions of the model for different values of $K_{\text {div}}[dens]$ and count the fraction of simulations keeping showing stochastic events all along the simulation

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. 4.2, parameters of Table 8 are used. In Sect. 4.3 where cell types are introduced, parameters of Table 9 are used. Adding $\hbox {O}_2$ and butyrate in Sect. 4.4, we use parameters presented in Table 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.

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

Table 7 Parameters for Sect. 4.1

Table 8 Parameters for Sect. 4.2

Table 9 Parameters for Sect. 4.3

Table 10 Parameters for Sect. 4.4

Table 11 Parameters for Sect. 4.5

Table 12 Parameters for the deterministic PDE model

1.1 Parameters of chemical kinetic

We could not find value for the molar concentration of oxygen. Therefore, we switched to normalized concentrations. We set $\bar{c_o} = \frac{c_o}{C}$ with C a normalizing constant. The evolution equation of oxygen concentration 6 becomes

$$\begin{aligned} \partial _t {\bar{c}}_o - \sigma _i \partial _{zz} {\bar{c}}_o = - \frac{s_o}{C} \gamma ^\infty _{\beta } \frac{{\bar{c_o}}^4c_b}{{\bar{c_o}}^4c_b+ \frac{K_{\beta }^5}{C^4}} (\nu _t^{\text {ent}} + \nu _t^{\text {gc}}) * \psi _{a}(z). \end{aligned}$$

(63)

We suppose that experimental measures of the reaction speed in Clausen and Mortensen (1994) 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,

$$\begin{aligned} K_{\beta }^5 = 0.184 \times c_o^{*4}. \end{aligned}$$

We suppose as well that oxygen concentration at the bottom of the crypt is $c_o^*$. The exact value of $c_o^*$ is unknown to us, but we set ${\bar{c}}_o^* = 10$. Then, $C = \frac{c_o^*}{{\bar{c}}_o^*} = \frac{c_o^*}{10}$. Therefore, we can find a value for the parameter $K_{\beta }^5/C^4$ of (63):

$$\begin{aligned} \frac{K_{\beta }^5}{C^4} = \frac{0.184 \times c_o^{*4} \times 10^4}{c_o^{4*}} = 0.184 \times 10^4. \end{aligned}$$

As C is unknown, the quantity $s_0/C$ in equation (63) is still unknown, but we decided to use $C=1$ in all simulations.