Probabilistic and Piecewise Deterministic models in Biology

TL;DR

This paper reviews recent advances in Piecewise Deterministic Markov Processes (PDMPs) and their applications in biological modeling, highlighting their theoretical challenges and diverse use cases in population dynamics, neuroscience, and gene expression.

Contribution

It provides new examples and insights into the long-term behavior, sampling methods, time scale separation, and moment calculus of PDMPs in various biological contexts.

Findings

Analysis of long time behavior in population models

Application of PDMPs to structured population dynamics

Use of moment calculus in gene expression models

Abstract

We present recent results on Piecewise Deterministic Markov Processes (PDMPs), involved in biological modeling. PDMPs, first introduced in the probabilistic literature by Davis (1984), are a very general class of Markov processes and are being increasingly popular in biological applications. They also give new interesting challenges from the theoretical point of view. We give here different examples on the long time behavior of switching Markov models applied to population dynamics, on uniform sampling in general branching models applied to structured population dynamic, on time scale separation in integrate-and-fire models used in neuroscience, and, finally, on moment calculus in stochastic models of gene expression.

Figures

Equations

$$\forall t\ge 0,{\partial }_t{X}_t=F^{(i)}(X_t),$$

$$F^{(i)}:x\mapsto a_{i}x,a_i\in \mathbb{R}.$$

$$F^{(i)}(x,y)=\left(\begin{array}{c}{\alpha }_{i}x(1-a_{i}x-b_{i}y)\\ {\beta }_{i}x(1-c_{i}x-d_{i}y)\end{array}\right),a_i,b_i,c_i,d_i,{\alpha }_i,{\beta }_i\in {\mathbb{R}}_{+}.$$

$$\forall t\ge 0,{\partial }_t{X}_t=a_{I_t}{X}_t.$$

$$\forall t\ge 0,X_t=X_0{e}^{{\int }_0^t{a}_{I_s}ds}.$$

$$\underset{t\to +\infty }{\lim }\frac{1}{t}{\int }_0^t{a}_{I_s}ds=\sum _{i\in F}{a}_i\nu (\{i\}).=\nu (a).$$

$$\forall q>0,\underset{t\to +\infty }{\lim }\mathbb{E}[X_t^q]=+\infty .$$

$$\forall q\le p,\underset{t\to +\infty }{\lim }\mathbb{E}[X_t^q]=0.$$

$$\mathbb{E}[X_t^p]=\mathbb{E}[X_0^p]\mathbb{E}\left[e^{{\int }_0^{t}p{a}_{I_s}ds}\right]=\mathbb{E}[X_0^p]{\mu }_0{e}^{t(A+pD)}1,$$

$$c_p{e}^{{\lambda }_{p}t}\le {\mu }_0{e}^{t(A+pD)}1\le C_p{e}^{{\lambda }_{p}t},$$

$$v_p(A+pD)={\lambda }_p{v}_p.$$

$$v_{p}D+{\partial }_p{v}_p(A+pD)=v_p{\partial }_p{\lambda }_p+{\lambda }_p{\partial }_p{v}_p.$$

$$\forall n\ge 0,Y_{n+1}={\Theta }_n{Y}_n,$$

$$\forall n\ge 0,Y_n=Y_0\prod _{k=0}^{n-1}{\Theta }_k=X_0{e}^{{\sum }_{k=0}^{n-1}\ln ({\Theta }_k)}.$$

$$F^{(0)}:v\mapsto \left(\begin{array}{cc}-1& 4\\ -1/4& -1\end{array}\right)\cdot v,F^{(1)}:v\mapsto \left(\begin{array}{cc}-1& -1/4\\ 4& -1\end{array}\right)\cdot v.$$

$$\forall t\geq 0,\quad\left\{\begin{array}[]{c @{=} c}x(t)&e^{-t}(\cos(t)x(0)-\sin(t)y(0)/4)\\ y(t)&e^{-t}(4\sin(t)x(0)+\cos(t)y(0)).\\ \end{array}\right.$$

$$\Vert (x(t),y(t))\Vert \le C{e}^{-t}\Vert (x(0),y(0))\Vert .$$

$$\underset{t\to +\infty }{\lim }\Vert X(t)\Vert =0,$$

$$⟨x-y,F^{(i)}(x)-F^{(i)}(y)⟩\le -\rho (i)\Vert x-y{\Vert }^2.$$

$$\forall t\ge 0,\Vert u^1(t)-u^2(t){\Vert }^2\le e^{-\rho (i)t}\Vert u^1(0)-u^2(0){\Vert }^2.$$

$$⟨u,{\mathbf{J}}_{x}G\cdot u⟩\le -\rho \Vert u{\Vert }^2,$$

$$⟨x-y,G(x)-G(y)⟩\le -\rho \Vert x-y{\Vert }^2.$$

$${\text{Hess}}_{x}V\ge \rho I,$$

$$\phi :t\mapsto \frac{⟨x-y,G(x)-G(ty+(1-t)x)⟩}{\Vert x-y{\Vert }^2}.$$

$$\sum _{i\in F}\rho (i)\nu (i)>0,$$

$$\forall t\ge 0,\mathbb{E}[V(X_t,I_t)]\le C{e}^{-\gamma t}V(X_0,I_0)+K.$$

$$\forall t\ge 0,{\partial }_t{X}_t=\sum _{i\in F}\nu (i)F^{(i)}(X).$$

$$Z_t=\sum _{u\in V_t}{\delta }_{X_t^u},t\ge 0,$$

$$\displaystyle m(x,s,t)=\mathbb{E}\left(N_{t}\big{|}Z_{s}=\delta_{x}\right),$$

$$B(x)\le b_1{\left|x\right|}^{\gamma }+b_2.$$

Full text

Introduction

The piecewise deterministic Markov processes (denoted PDMPs) were first introduced in the literature by Davis [30, 31]. The key point of his work was to endow the PDMP with rather general tools, similar as such that already exist for diffusion processes. Indeed, PDMPs form a family of non-diffusive càdlàg Markov processes, involving a deterministic motion punctuated by random jumps. The motion of the PDMP $\{X(t)\}_{t\geq 0}$ depends on three local characteristics, namely the jump rate $\lambda$, the deterministic flow $\varphi$ and the transition measure $Q$ according to which the location of the process at the jump time is chosen. The process starts from $x$ and follows the flow $\varphi(x,t)$ until the first jump time $T_{1}$ which occurs either spontaneously in a Poisson-like fashion with rate $\lambda(\varphi(x,t))$ or when the flow $\varphi(x,t)$ hits the boundary of the state-space. In both cases, the location of the process at the jump time $T_{1}$, denoted by $Z_{1}=X(T_{1})$, is selected by the transition measure $Q(\varphi(x,T_{1}),\cdot)$ and the motion restarts from this new point as before. This fully describes a piecewise continuous trajectory for $\{X(t)\}_{t\geq 0}$ with jump times $\{T_{k}\}$ and post jump locations $\{Z_{k}\}$, and which evolves according to the flow $\varphi$ between two jumps.

Since the seminal work of Davis, PDMPs have been heavily studied from the theoretical perspective, we may refer for instance the readers to [27, 2, 1, 57, 7, 9] among many others.

From the applied point of view, and more precisely in biological applications, these processes are sometimes referred as hybrid, and are especially appealing for their ability to capture both continuous (deterministic) dynamics and discrete (probabilistic or deterministic) events. First applications date back at least to the studies of the cell cycle model [53, 52], and more recent works, to name just a few, in neurobiology [63, 35, 35], cell population and branching models [4, 22, 60, 36, 51, 20], gene expression models [72, 56, 33], food contaminants [14, 12] or multiscale chemical reaction network models [28, 3, 50, 46]. See also [47, 16, 66] for selected reviews on hybrid models in biology.

The paper is organized as follows. In Section 1, we present long time behavior results for switched flows in population dynamics. In Section 2, we present many-to-one formulas and description of the trait of an individual uniformly sampled in a general branching model applied to structured population dynamic. In Section 3, we present limit theorems and time scale separation in integrate-and-fire models used in neuroscience. Finally, in Section 4, we derive asymptotic moments formulas in a stochastic model of gene expression.

1. Long time behavior of some PDMP for population dynamics

In this section, we present recent results on long time behavior and exponential convergence of some PDMP used in population dynamics, in particular, for modeling population growth in a varying environment.

We consider processes $(X_{t},I_{t})_{t\geq 0}$, where the first component $X_{t}$ has continuous paths in a continuous space $E$, namely a subset of $\mathbb{R}^{d}$, for some integer $d$, and the second component $I_{t}$ is a pure jump process on a finite state space $F$ or in $\mathbb{N}$. The continuous variable represents the number of individuals (or more precisely the population density or the relative abundance) of a population and the discrete variable models the population environment. In a fixed environment $I_{t}=i\in F$, we assume that $(X_{t})_{t\geq 0}$ evolves as the solution of an ordinary differential equation, that is

[TABLE]

where $F^{(i)}$ is some smooth function. For instance, the oldest and maybe the most famous model to represent the evolution of a population is the choice $E=\mathbb{R}_{+}$ and

[TABLE]

When there is no variability in the environment, solutions are given by exponential functions and there is either explosion or extinction according to the sign of $a_{i}$. We develop in Subsection 1.1 this toy model when the environment varies. Almost all properties (almost-sure convergence, moments…) can be derived and it gives a first hint to understand more complicated models such as multi-dimensional linear models or non-linear models. These generalizations are described in Subsection 1.2, and 1.3, respectively. For non-linear population models, maybe one of the most famous models is given by the Lotka-Volterra equations (also known as the predator-prey equations), with state-space $E=\mathbb{R}_{+}^{2}$ and

[TABLE]

To end the introduction of this section, let us mention that all of these models are particular instance of a class of switching Markov model, see for instance [9, 71] for general references. Some other examples of application are developed in [57, 38, 61].

1.1. Malthus model

Let us consider, in this subsection, that $(I_{t})_{t\geq 0}$ is an irreducible continuous time Markov chain on a finite state space $F$. We denote by $\nu$ its unique invariant probability measure. The process $(X_{t})_{t\geq 0}$ will be the solution of

[TABLE]

This means that $F^{(i)}$ is given by (1.2). Thus, we have

[TABLE]

The almost-sure behavior of this process is easy to understand. Indeed, by the ergodic theorem, we have that

[TABLE]

Then, we have the following dichotomy: either $\nu(a)>0$ and $(X_{t})_{t\geq 0}$ tends almost-surely to infinity at an exponential rate, or $\nu(a)<0$ and $(X_{t})_{t\geq 0}$ tends to [math]. The second statement can be understood as the extinction of the population, even if, in contrast with classical birth and death processes, $(X_{t})_{t\geq 0}$ never hits [math]; that is the extinction time (or the hitting time of [math]) is almost-surely infinite.

The convergence of moments of $(X_{t})_{t\geq 0}$ is more tricky than the almost-sure convergence. Indeed, one cannot use a dominated convergence theorem or a Jensen inequality (which is in the wrong way). Nevertheless, we have

Lemma 1.1** (Convergence of moments).**

If $\nu(a)>0$ then

[TABLE]

If $\nu(a)<0$ then there exists $p>0$ such that

[TABLE]

The proof which follows comes from [6]; see also [32].

Proof.

If we denote by $A$ the generator of the process $(I_{t})_{t\geq 0}$ and $\mu_{0}$ the vector corresponding to its initial law, then the Feynman-Kac formula states that, for every $p>0$,

[TABLE]

where it is the classical exponential of matrices, $D$ is the diagonal matrix such that $D_{i,i}=a_{i}$ and $\mathbf{1}$ is the vector made of ones. Now, classical results on matrices (including Perron-Frobenius Theorem) give the existence of two constants $c_{p},C_{p}>0$ such that

[TABLE]

where $\lambda_{p}$ is the (unique) eigenvalue of largest real part of the matrix $(A+pD)$. It then remains to understand the sign of $\lambda_{p}$ to understand the moment behaviour. By definition of $\lambda_{p}$, there exists $v_{p}$ such that

[TABLE]

Moreover, for $p=0$, we have $\lambda_{0}=0$ and $v_{0}=\nu$. Again classical results on matrices ensure that $p\mapsto\lambda_{p}$ and $p\mapsto v_{p}$ are smooth and then differentiating in $p$, we obtain

[TABLE]

In particular, for $p=0$ and multiplying by $\mathbf{1}$ by the right, we have $\nu(a)=(\partial_{p}\lambda_{p})_{\scriptscriptstyle{|p=0}}$. As a consequence, if $\nu(a)>0$, then $p\mapsto\lambda_{p}$ is increasing in a neighborhood of [math] and for all $p>0$ small enough we have $\lambda_{p}>0$. Thus, $\mathbb{E}[X_{t}^{p}]$ tends to infinity for $p$ small enough, and then for all $p>0$ by the Jensen inequality. In contrary, if $\nu(a)<0$ then there exists $p>0$ such that $\lim_{t\to+\infty}\mathbb{E}[X_{t}^{p}]=0$ and also for all $q<p$, again by Jensen inequality. ∎

Understanding the long time behavior of moments is primordial when we are interested in the convergence to some invariant law more general than $\delta_{0}$. Indeed, moments provide Lyapunov-type functionals ensuring that the process returns frequently in a compact-set.

1.2. Others linear models

From a modeling point of view, the results of the previous section are easily understandable. Indeed, it is enough that the environment has in mean a positive effect to ensure an exponential growth, while if the environment is, in mean, not favorable, then the population goes to extinct.

However, note that this interpretation only holds because the growth parameter is well-chosen. Indeed, let us first consider a discrete time analogous of the Malthus model: the sequence $(Y_{n})_{n\geq 0}$ satisfying

[TABLE]

where $(\Theta_{n})_{n}$ is a sequence of i.i.d. non-negative random variables. Here $\Theta_{n}$ represents the mean number of children per individuals at generation $n$ (all individuals have the same number of children which depends on the environment). If $(\Theta_{n})_{n\geq 0}$ is a constant sequence then $Y_{n}=\Theta_{1}^{n}Y_{0}$ and a dichotomy occurs with respect to a critical value $\Theta_{1}=1$. When $(\Theta_{n})_{n\geq 0}$ is a sequence of i.i.d random variable, then there is also a dichotomy but the critical value (to be smaller or larger than $1$) is no longer $\mathbb{E}[\Theta_{1}]$ but $\exp\left(\mathbb{E}[\ln(\Theta_{1})]\right)$. Indeed, we have

[TABLE]

This is exactly the same phenomenon for Galton-Watson chain in random environment, see for instance [5] and [42, Section 2.9.2].

There is a similar (but more complex, maybe) problem in upper dimension. Let us consider that $E=\mathbb{R}^{2}$, $F=\{0,1\}$, and $(I_{t})_{t\geq 0}$ jumps from [math] to $1$ and from $1$ to [math] with the same rate $\lambda>0$, and

[TABLE]

Then the solutions of the equation $\partial_{t}\begin{pmatrix}x\\ y\end{pmatrix}=F^{(1)}\begin{pmatrix}x\\ y\end{pmatrix}$ are given by

[TABLE]

It is easy to see that there exists some constants $C>0$ such that for all $t\geq 0$,

[TABLE]

This property is also satisfied by the solutions of $\partial_{t}\begin{pmatrix}x\\ y\end{pmatrix}=F^{(0)}\begin{pmatrix}x\\ y\end{pmatrix}$. However, we have the following theorem.

Theorem 1.2** ([8]).**

Let $(X_{t})_{t\geq 0}$ be the solution of Equation (1.1) with flows defined in (1.3). There exists $\beta_{c}>0$ such that if $\lambda<\beta_{c}$ then

[TABLE]

and if $\lambda>\beta_{c}$ then $\lim_{t\to+\infty}\|X(t)\|=+\infty$. Moreover, the speed of convergence is exponential.

For this type of results for more general vector fields, see for instance [8, 57, 54]. Let us briefly explain this phenomenon. One can see that, whatever the initial conditions of the solutions $(x,y)$ of (1.4), the map $t\mapsto\|(x(t),y(t))\|$ is not decreasing; see for instance Figure (2).

Then, heuristically, we see that if $I$ jumps sufficiently fast (before the decay of the distance) then the distance may grow. However, if $I$ is too slow, then $(X_{t})_{t\geq 0}$ essentially follows one of the two flows and goes rapidly to [math].

Switching linear models can have even more unpredictable behavior than the previous one. Indeed, in [54], the author exhibits two different linear functions $F^{(0)},F^{(1)}$ on $\mathbb{R}^{2}$ and two critical values $\beta_{2}>\beta_{1}$ such that the switched process $(X_{t})_{t\geq 0}$ tends to infinity for at least one $\lambda\in(\beta_{1},\beta_{2})$ and to [math] if $\lambda\in(0,\beta_{1})\cup(\beta_{2},+\infty)$.

1.3. A general theorem for convergence

In this subsection, we expose a condition similar to Subsection 1.1 for ensuring convergence to an equilibrium (not necessary $\delta_{0}$). As pointed out in the previous section, we have to measure precisely the growth of each underlying flow. So we assume that for all $i\in F$, there exists $\rho(i)\in\mathbb{R}$ such that

[TABLE]

The parameter $\rho(i)$ measures the contraction of a solution $\dot{u}=F^{(i)}(u)$ with respect to the norm $\|\cdot\|$. Indeed, for two solutions, $u^{1}$ and $u^{2}$, we have

[TABLE]

Note that the next result does not rely on a specific norm, but it has to be the same for all flows. A sufficient condition to prove (1.5) is given by the next lemma.

Lemma 1.3**.**

Let $G:\mathbb{R}^{d}\to\mathbb{R}^{d}$ be a $C^{1}$ map and let $\mathbf{J}_{x}G$ be its Jacobian matrix at the point $x\in\mathbb{R}^{d}$. If for all $u,x\in\mathbb{R}^{d}$,

[TABLE]

then, for all $x,y\in\mathbb{R}^{d}$,

[TABLE]

In particular, if $G=\nabla V$, for some $C^{2}$ function $V$, then (1.6) holds under the classical condition that $V$ is strongly convex with constant $\rho$; namely for all $x\in\mathbb{R}^{2}$,

[TABLE]

in the sense of quadratic form.

Proof.

It is classic to derive the previous bounds using the mean value theorem for the following map

[TABLE]

∎

We can now state the main result of this subsection.

Theorem 1.4** ([24]).**

Assume that $(I_{t})_{t\geq 0}$ is an irreducible Markov process on a finite state space $F$ with invariant measure $\nu$, $(X_{t})_{t\geq 0}$ satisfies (1.1), and that Eq. (1.5) holds. If

[TABLE]

then $(X_{t},I_{t})$ converges in law to a unique invariant law.

The proof of this theorem is given in [24, 7]. The proof of [24] is based on a generalization (developed in [44]) of the classical Harris Theorem [43], which does not hold in general for PDMP (see for instance the introduction of [7]). One of the main point of the proof is the construction of a Lyapunov function $V$, in the sense that there exist, $C,\gamma,K>0$ such that

[TABLE]

The construction of $V$ and the control of its moments are based on Lemma 1.1.

Note that another proof of the previous theorem is given in [7]. Nevertheless, the proof of [24] takes the advantage to be generalizable in many contexts (dependent environment, non-deterministic underlying dynamics, asymptotic assumptions…). See [24] for details.

Closely related articles [10, 58, 59] deal with a fluctuating Lotka-Volterra model. They show that random switching between two environments both favorable to the same specie can lead to the extinction of this specie or coexistence of both species.

As Theorem 1.4, their proofs are based on the construction of a Lyapunov function. Nevertheless instead of controlling the exit from a compact set of $\mathbb{R}$ to infinity, this Lyapunov function is used to control the exit from a compact set included in $(0,\infty)^{2}$ to the axes $\{(0,y)|y\geq 0\}$ or $\{(x,0,)|x\geq 0\}$. A complete description of this model is given in [58]. Note however that their counter-intuitive result is similar to the one of Theorem 1.2. It is based on the fact that, when $I$ jumps sufficiently fast (namely $\lambda$ being large in Theorem 1.2), then $X$ is close to follow the following deterministic dynamics:

[TABLE]

Finally, we end this section mentioning that such hybrid framework can also be applied to models where the population is represented by a discrete variable, and/or the environment fluctuates continuously. One can cite for instance [26], where the author studies a model of interaction between a population of insects and a population of trees. Another example is the chemostat model (a type of bioreactor), introduced in [29] and studied in [25, 19, 18, 23]. Lyapunov functions are again a key theoretical tool to study long time behavior for such models. However, we point out that extinction events may occur in discrete population models, and that the interplay between the discrete and continuous components makes the problem more difficult than standard absorbing time studies in purely discrete population models. Nevertheless, it was proved in [23, Theorem 3.4] that extinction event is almost-surely finite and has furthermore exponential tails. This generalizes part of a result of [25].

2. Uniform sampling in a branching structured population

In this section, which is a short version of [60], we give a characterization of the trait of a uniformly sampled individual among the population at time $t$ in a structured branching population.

2.1. Introduction

We consider a branching Markov process where each individual $u$ is characterized by a trait $\left(X_{t}^{u},\ t\geq 0\right)$ which dynamic follows an $\mathcal{X}$-valued Markov process, where $\mathcal{X}:=\widetilde{\mathcal{X}}\times\mathbb{R}_{+}$ and $\widetilde{\mathcal{X}}\subset\left(\mathbb{R}_{+}\right)^{d}$ for some $d\geq 1$. For any $x\in\mathcal{X}$, the last coordinate corresponds to a time coordinate and we denote by $\widetilde{x}\in\widetilde{\mathcal{X}}$ the vector of the first $d$ coordinates. We assume that this trait influences the lifecycle of each individual (in terms of lifetime, number of descendants and inheritance of the trait). Then, an interesting problem consists in the characterization of the trait of a ”typical” individual in the population. This is the question we address here.

Let us describe the process in more details. We consider a strongly continuous contraction semi-group with associated infinitesimal generator $\mathcal{G}:\mathcal{D}(\mathcal{G})\subset\mathcal{C}_{b}(\mathcal{X})\rightarrow\mathcal{C}_{b}(\mathcal{X})$, where $\mathcal{C}_{b}(\mathcal{X})$ denotes the space of continuous bounded functions from $\mathcal{X}$ to $\mathbb{R}$. Then, each individual $u$ has a trait $(X_{t}^{u},\ t\geq 0)$ which evolves as a Markov process defined as the unique $\mathcal{X}$-valued càdlàg solution of the martingale problem associated with $(\mathcal{G},\mathcal{D}(\mathcal{G}))$. An individual with trait $x$ dies at an instantaneous rate $B(x)$, where $B$ is a continuous function from $\mathcal{X}$ to $\mathbb{R}_{+}$. It is replaced by $A_{u}(x)$ children, where $A_{u}(x)$ is a $\mathbb{N}$-valued random variable with distribution $\left(p_{k}\left(x\right),k\geq 0\right)$. For convenience, we assume that $p_{1}(x)\equiv 0$ for all $x\in\mathcal{X}$. The trait of the descendants at birth depends on the trait of the mother at death. For all $k\in\mathbb{N}$, let $P^{(k)}(x,\cdot)$ be the probability measure on $\mathcal{X}^{k}$ corresponding to the trait distribution at birth of the $k$ descendants of an individual with trait $x$. We denote by $P_{j}^{(k)}\left(x,\cdot\right)$ the $j$th marginal distribution of $P^{(k)}$ for all $k\in\mathbb{N}$ and $j\leq k$.

Finally, we denote by $\mathcal{M}_{P}(\mathcal{X})$ the set of point measures on $\mathcal{X}$. Following Fournier and Méléard [39], we work in $\mathbb{D}\left(\mathbb{R}_{+},\mathcal{M}_{P}\left(\mathcal{X}\right)\right)$, the state of càdlàg measure-valued Markov processes. For any $Z\in\mathbb{D}\left(\mathbb{R}_{+},\mathcal{M}_{P}\left(\mathcal{X}\right)\right)$, we write:

[TABLE]

the measure-valued process describing the dynamic of the population where $V_{t}$ denotes the set of all individuals in the population at time $t$. We will denote by $N_{t}$ the cardinality of $V_{t}$ and by:

[TABLE]

its expected value, for $0\leq s\leq t$.

An example of such a process is the size-structured population where each individual grows exponentially fast. For more details we refer the reader to [36] or [13]. Then, at rate $B(X_{t}^{u})$, the individual $u$ splits at time $t$ in two daughter cells of size $X_{t}^{u}/2$.

Figure 3 is a realization of such a process. The diameter of each circle corresponds to the size of the individual at division and the length of the branch represents the lifetime of each individual. In particular, we notice the link between the lifetime and the size of each individual: the bigger the cell is, the shorter its lifetime is.

In Section 2.2, we give the chosen assumptions on the model to ensure that our process is well-defined. Then, in Section 2.3, we describe the process corresponding to the trait of a ”typical” individual in the population. Finally, in Section 2.4, we explain why this process corresponds to the trait of a uniformly sampled individual in a large population approximation.

2.2. Definition and existence of the structured branching process

To define rigorously the structured branching process on $\mathbb{R}_{+}$, we need to ensure that almost surely the number of individuals in the population does not blow up in finite time.

As explained in the previous description of the process, the lifetime of each individual depends on the dynamic of the trait through the division rate $B$. One way of dealing with this dependency is to assume that the division rate is upper bounded by some constant $\overline{B}$. Then, in the case of a binary division, the expected number of individuals in the population at time $t$ is upper bounded by the expected value of a Yule process with birth rate $\overline{B}$ at time $t$ which is equal to $\exp(\overline{B}t)$ and which is bounded on compact sets. In the case of an unbounded division rate, the same reasoning applies if we can control the excursions of the dynamic of the trait. Thus, we consider two sets of hypotheses: the first one controls what happens regarding divisions (in term of rate of division and of mass creation) and the second one controls the dynamic of the trait between divisions.

Assumption 1**.**

We consider the following assumptions:

1. (1)

There exist $b_{1},b_{2}\geq 0$ and $\gamma\geq 0$ such that for all $x\in\mathcal{X}$,

[TABLE] 2. (2)

There exists $\underline{x}\in\widetilde{\mathcal{X}}$ such that for all $x\in\mathcal{X}$, $k\in\mathbb{N}$:

[TABLE] 3. (3)

There exists $\overline{m}\geq 0$ such that for all $x\in\mathcal{X}$,

[TABLE]

The first point controls the division rate. The second point means that we consider a fragmentation process with a possibility of mass creation at division when the mass is small enough. In particular, clones are allowed in the case of bounded traits and bounded number of descendants and any finite type branching structured process can be considered.

As explained before, we make a second assumption to control the behavior of traits between divisions.

Assumption 2**.**

There exist $c_{1},c_{2}\geq 0$ such that for all $x\in\mathcal{X}$:

[TABLE]

where $\gamma$ is defined in Assumption 1 and for $x\in\left(\mathbb{R}_{+}\right)^{d}$, $h_{\gamma}(x)=|x|^{\gamma}=\left(\sum_{i=1}^{d}x_{i}\right)^{\gamma}$.

This assumption ensures in particular that the trait does not blow up in finite time in the case of a deterministic dynamic for the trait. Assumptions 1(1) and 2 are linked via the parameter $\gamma$ which controls the balance between the growth of the population and the dynamic of the trait. In particular, if Assumption 1 is satisfied for $\gamma=0$, the division rate is bounded and Assumption 2 is always satisfied.

Under Assumptions 1 and 2, the previously described measure-valued process is well-defined as the strongly unique solution of a stochastic differential equation.

The proof of existence relies on a recursive construction of a solution via the sequence of jumps time of the population. Then, we prove that this sequence is unique conditionally to the Poisson point measure determining the jumps in the population and to the stochastic flows corresponding to the dynamic of the trait. Finally, using Assumptions 1 and 2, we prove that the number of individuals in the population does not blow up in finite time. This is detailed in [60] (Theorem 2.3.).

2.3. The trait of sampled individuals at a fixed time : Many-to-One formula

In order to characterize the trait of a uniformly sampled individual, the spinal approach ([21],[55]) consists in following a ”typical” individual in the population whose behavior summarizes the behavior of the entire population. In this section, we give the dynamic of the trait of a typical individual: it is a Markov process called the auxiliary process.

From now on, we assume that for all $x\in\mathcal{X}$, $t\geq 0$ and $s\leq t$, $m(x,s,t)\neq 0$. With a slight abuse of notation, we denote by $X_{s}^{u}$ the trait of the unique ancestor living at time $s$ of $u$.

Let $f$ be a non-negative measurable function and $x\in\mathcal{X}$. The idea is to notice that the following operator:

[TABLE]

is a conservative (non-homogeneous) semi-group. Then, the auxiliary process is defined as the time-inhomogeneous Markov process with associated family of semi-groups $\left(P_{r,s}^{(t)},r\leq s\leq t\right)$. It describes the dynamic of the trait of a ”typical” individual in the population in the sense that it satisfies a so-called ”Many-to-One” formula. This formula is true under Assumptions 1, 2 and under the following technical assumptions:

Assumption 3**.**

There exists a function $C$ such that for all $j\leq k,\ j,k\in\mathbb{N}$ and $0\leq s\leq t$, we have:

[TABLE]

This assumption tells us that we control uniformly in $x$ the benefit or the penalty of a division.

Assumption 4**.**

For all $t\geq 0$ and $x\in\mathcal{X}$ we have:

$s\mapsto m(x,s,t)$* is differentiable on $[0,t]$ and its derivative is continuous on $[0,t]$,*

• -

for all $f\in\mathcal{D}(\mathcal{A}):=\left\{f\in\mathcal{D}(\mathcal{G})\text{ s.t. }m(\cdot,s,t)f\in\mathcal{D}(\mathcal{G})\ \forall t\geq 0,\ \forall s\leq t\right\}$, $s\mapsto\mathcal{G}(m(\cdot,s,t)f)(x)$ is continuous,

We can now state the Many-to-One formula.

Theorem 2.1**.**

Under Assumptions 1, 2, 3 and 4, for all $t\geq 0$, for all $x_{0}\in\mathcal{X}$ and for all non-negative measurable functions $f:\mathcal{X}\rightarrow\mathbb{R}_{+}$, we have:

[TABLE]

where $\left(Y_{s}^{(t)},s\leq t\right)$ is an inhomogeneous-Markov process with infinitesimal generator $\left(\mathcal{A}_{s}^{(t)},\mathcal{D}(\mathcal{A})\right)$ given for $f\in\mathcal{D}(\mathcal{A})$ and $x\in\mathcal{X}$ by:

[TABLE]

where:

[TABLE]

[TABLE]

[TABLE]

Comments on the proof.

To prove the Many-to-One formula (2.2), we first show that the family of semi-groups $\left(P_{r,s}^{(t)},r\leq s\leq t\right)$ is uniquely defined as the unique solution of an integro-differential equation (see Lemma 3.2 in [60]). In particular, we need Assumption 3 to prove the uniqueness. Then, the infinitesimal generator $\left(\left(\mathcal{A}_{s}^{(t)}\right)_{s\leq t},\mathcal{D}\left(\mathcal{A}\right)\right)$ of the auxiliary process is obtained by differentiation of the semi-group $\left(P_{r,s}^{(t)},r\leq s\leq t\right)$ using Assumption 4.

Remark 2.2**.**

We can also prove a pathwise version of (2.2) using a monotone class argument (Theorem 3.1. in [60]).

Unlike previous works on this subject ([41], [45], [22]), the existence of our auxiliary process does not rely on the existence of spectral elements for the mean operator of the branching process. In particular, we can apply this result to models with a varying environment. The dynamic of this auxiliary process heavily depends on the comparison between $m(x,s,t)$ and $m(y,s,t)$, for $x,y\in\mathcal{X}$. It emphasizes several bias due to growth of the population. First, the auxiliary process jumps more than the original process, if jumping is beneficial in terms of number of descendants. This phenomenon of time-acceleration also appears for examples in [21], [55] or [45]. Moreover, the reproduction law favors the creation of a large number of descendant as in [4] and the non-neutrality favors individuals with an ”efficient” trait at birth in terms of number of descendants. Finally, a new bias appears on the dynamic of the trait because of the combination of the random evolution of the trait and non-neutrality. Indeed, if the dynamic of the trait is deterministic, we have $\widehat{\mathcal{G}}_{s}^{(t)}f(x)=\mathcal{G}f(x)$.

2.4. Ancestral lineage of a uniform sampling at a fixed time in a large population

The Many-to-One formula (2.2) gives the law of the trait of a typical individual in the population. But so far, we did not prove that it corresponds to the trait of a uniformly sampled individual. In particular, we have to take into account that the number of individuals in the population is stochastic and depends on the dynamic of the trait. Using the law of large numbers, we can approximate the number of individuals in a population starting for $n$ individuals, divided by $n$, by the mean number of individual in the population. That is why we now look at the ancestral lineage of a uniform sampling in a large population.

It only makes sense to speak of a uniformly sampled individual at time $t$ if the population does not become extinct before time $t$. For all $t\geq 0$, let $\Omega_{t}=\left\{N_{t}>0\right\}$ denote the event of survival of the population. Let $\nu\in\mathcal{M}_{P}(\mathcal{X})$ be such that:

[TABLE]

We set

[TABLE]

where $X_{i}$ are i.i.d. random variables with distribution $\nu$. For $t\geq 0$, we denote by $U(t)$ the random variable with uniform distribution on $V_{t}$ conditionally on $\Omega_{t}$ and by $\left(X^{U(t)}_{s},\ s\leq t\right)$ the process describing the trait of a sampling along its ancestral lineage. If $X$ is a stochastic process, we denote by $X^{\nu}$ the process with initial distribution $\nu\in\mathcal{M}_{P}(\mathcal{X})$.

Theorem 2.3**.**

Under Assumptions 1,2, 3 and 4, for any $t\geq 0$, the following convergence in law in $\mathbb{D}\left([0,t],\mathcal{X}\right)$ holds:

[TABLE]

Remark 2.4**.**

If we start with $n$ individuals with the same trait $x$, we obtain:

[TABLE]

Therefore, the auxiliary process describes exactly the dynamic of the trait of a uniformly sampled individual in the large population limit, if all the starting individuals have the same trait. If the initial individuals have different traits at the beginning, the large population approximation of a uniformly sampled individual is a linear combination of the auxiliary process.

3. Averaging for some integrate-and-fire models

3.1. The model

We examine the following slow-fast hybrid version of integrate-and-fire models [17], used in mathematical neuroscience. In such a setting, $X$ represents the membrane potential of a neural cell which is increasing until it reaches some threshold $c$, corresponding to the time where a nerve impulse is triggered, and then the potential is reset to some lower value. More precisely, the process $(X(t),t\in[0,T])$, with $T$ some time horizon, obeys the following dynamic:

1. (1)

Initial state: At time $T^{\ast}_{0}=0$, the process starts at $X(T^{\ast}_{0})=\xi_{0}$, a random variable with support included in $(m,c)$ where $\{c\}$ is considered as a boundary and $m<c$ is some real. 2. (2)

First jumping time: Let $Y$ be a continuous time Markov chain valued in a countable space ${\mathcal{Y}}$. This chain starts at $Y(0)=\zeta$, a $\mathcal{Y}$-valued random variable. The first hitting time of the boundary occurs at time $T^{\ast}_{1}$ defined as

[TABLE]

where $\alpha$ is a positive measurable function such that $\alpha(\mathcal{Y})$ is bounded from above and $F$ is a positive continuous function. 3. (3)

Piecewise deterministic motion: For $t\in[T^{\ast}_{0},T^{\ast}_{1})$, we set

[TABLE]

The dynamic of $X$ is thus continuous here, and given by the differential equation:

[TABLE] 4. (4)

Jumping measure: Then, at time $T^{\ast-}_{1}$, the process is constrained to stay inside $(m,c)$ by jumping according to the $Y$-dependent measure $\mu_{Y_{T^{\ast-}_{1}}}$ whose support is included in $(m,c)$:

[TABLE] 5. (5)

And so on: Go back to step 1 in replacing $T^{\ast}_{0}$ by $T^{\ast}_{1}$ and $\xi_{0}$ by $\xi_{1}=X(T^{\ast}_{1})$.

This is not hard to see that the process $(X,Y)$ is a piecewise deterministic Markov process. Moreover, due to the presence of a boundary $c$, this is also a constrained Markov process. Our aim is, in this quite simple situation, to describe the behavior of the process $X$ at the boundary $c$, when the dynamic of the underlying process $Y$ is infinitely accelerated.

3.2. Acceleration and averaged process

From now on, we assume that the underlying process of celerities, that is the continuous time Markov chain $Y$, has a fast dynamic, by introducing a (small) timescale parameter $\varepsilon$ such that

[TABLE]

In the same time, to insure a limiting behavior, we assume that $Y$ is positive recurrent with intensity matrix $Q=(q_{zy})_{z,y\in\mathcal{Y}}$ and invariant probability measure $\pi$, considered as a row vector. For convenience, let us also define by $\alpha(V)$ the diagonal matrix such that $\text{diag}(\alpha(V))=\{\alpha(y)~{};~{}y\in{\mathcal{Y}}\}$.

As $\varepsilon$ goes to zero, the process $Y_{\varepsilon}$ converges towards the stationary state associated to $Y$ in the sense that, by the ergodic theorem,

[TABLE]

Therefore, as $\varepsilon$ goes to zero, the process $X_{\varepsilon}$, defined as $X$ by replacing $Y$ by $Y_{\varepsilon}$, should have its dynamic averaged with respect to the measure $\pi$. The behavior of the limiting process away from the boundary is indeed not hard to describe, see for example [49] and references therein. Here and later on, $\mathbb{D}\left(\left[0,\tau\right],\mathbb{R}\right)$ is the space of càdlàg functions from $[0,\tau]$ to $\mathbb{R}$, with $\tau$ a time horizon. For technical reasons, we assume that

• •

$\frac{1}{F}$ is integrable over $(m,c)$,

• •

$\alpha({\mathcal{Y}})$ is bounded from above,

• •

$\mathbb{E}(\sup_{\varepsilon\in(0,1]}p^{\ast}_{\varepsilon}(T))$ is finite, where $p^{\ast}_{\varepsilon}(T)$ is the counting measure at the boundary for the process $X_{\varepsilon}$.

Proposition 3.1**.**

Assume that $\xi_{0}$ is deterministic. Then, for any $\eta>0$, the process $X_{\varepsilon}$ converges in law towards a process $\bar{X}$ on $\mathbb{D}\left(\left[0,\frac{c-\xi_{0}}{\max\alpha(\mathcal{Y})}-\eta\right],\mathbb{R}\right)$ defined as:

[TABLE]

Proposition 3.1 gives the behavior of the limiting process away from the boundary: celerities are averaged against the measure $\pi$. But what happens at boundary? This is what we want to characterize. The following result holds.

Theorem 3.2**.**

The process $X_{\varepsilon}$ converges in law in $\mathbb{D}([0,T],\mathbb{R})$ towards a process $\bar{X}$ such that for any sufficiently regular function $f:(m,c)\to\mathbb{R}$ the process

[TABLE]

for $t\in[0,T]$, defined a martingale, with $\bar{p}^{\ast}$ the counting measure at the boundary for $\bar{X}$. The averaging measure at the boundary $\bar{\mu}$ is defined by

[TABLE]

where, for $y\in\mathcal{Y}$, $\pi^{*}$ is given by

[TABLE]

We can read in (3.1) the behavior of the averaged process $\bar{X}$. Away from the boundary, the dynamic of $\bar{X}$ is the one of the averaged flow, as indeed described in Proposition 3.1. At the boundary, the averaged jump measure is not given by the averaged version of the jumping measure, that is

[TABLE]

but the values of the celerity at the boundary are taken into account through the presence of $\alpha$. There is a balance between the celerity intensities and the invariant probabilities. For example, consider a piecewise linear motion whose speed switches at rate $1/\varepsilon$ between $1$ and $2$. When, $\varepsilon$ goes to zero, the process will spend as much time increasing with celerities $1$ and $2$, but when increasing with celerity $2$, the process obviously increases twice more than with celerity one in a same time window, and then approaches the boundary twice more rapidly, and therefore will have in fact twice more chance to hit the boundary with this celerity.

Note that because of the separation of variables in the form of the flow, the value of the process $\bar{X}$ at the boundary does not appear in the expression of $\pi^{\ast}$, as it could be the case in more general situations. Theorem such as 3.2 can be derived using the Prohorov program : prove the tightness of the family $\{X_{\varepsilon},\varepsilon\in(0,1)\}$ and then identify the limit. Tightness for constrained Markov processes, as the process described in the present section, has an interesting and rich literature, see [48] and references therein. We refer to [40] for more details about the proof of Theorem 3.2.

4. Stochastic models of protein production with cell division and gene replication

The production of proteins is the main mechanism of the prokaryotic cells (like bacteria). These functional molecules represent about $3.6\times 10^{6}$ molecules of $2000$ different types and compose about half of the dry mass of the cell [15] and this pool needs to be duplicated in the cell cycle time. In total, it is estimated that $67\%$ of the energy of the cell is dedicated to the protein production [67]. The protein production is a two steps process by which the genetic information is first transcribed into an intermediate molecule, the messenger-RNA (mRNA), and then translated into a protein. Experimental studies have shown that this production is subject to a large variability [37, 62, 68]. This heterogeneity can either be explained by the stochastic nature of the production mechanism (both transcription and translation are due to random encounter between molecules) or the effect of cell scale phenomenon like the DNA replication, cell division, or fluctuation in the resources.

In particular, the extensive experimental study of Taniguchi et al. (2010) [68] measured the variability of protein concentration for more than $1000$ different proteins in E. coli. For each of them, the authors have measured the average protein concentration, the average mRNA concentration and their lifetime. They interpret different sources of proteins, either from the protein production mechanism itself (i.e. the transcription of mRNAs and the translation of proteins) or from external factors, like the gene replication, division or the sharing of common resources in the production.

One can interpret the experimental results in the light of classical stochastic models of protein production [11, 65, 64] that predict the variability of protein production. But these classical models do not consider important phenomenon that may have an impact on the production: they usually do not consider the replication of the gene at some point in the cell cycle and also do not take into account the random assignment of each protein in either of the two daughter cells at division. Moreover, they only consider the number of mRNAs and proteins without any explicit notion of volume, and therefore the notion of concentration is unclear in this case. Since the experimental study only considers concentrations, it would be difficult to have quantitative comparisons between the variances predicted by these models and the ones experimentally measured.

We therefore propose here a more realistic model of protein production that takes into account all the basic features that can be expected for the production of a type of protein inside a cell cycle: the transcription, the translation, the gene replication and the division. We also explicitly represent the number of mRNAs and proteins in the growing volume of the bacteria, so that we can consider the concentration of each of these entities in the cell. We will then estimate the theoretical contribution of these features to the protein variance through mathematical analysis. We will also make a biological interpretation of these results with respect to experimental measures.

4.1. Model and Results

Our model presented here is adapted from a simple gene constitutive model that does not consider any gene regulation (it is a particular case of the model presented in [64]). Like this classical model, our model is gene centered (it represents the production of only one protein), and has two stages (transcription and translation) where both the number of mRNAs, $M$, and proteins, $P$, are explicitly represented. Each event (creation or degradation of an mRNA, creation of a protein, etc.) is supposed to occur at random times. The time intervals are considered as exponentially distributed. In addition to this classical framework, we introduce a notion of cell cycle represented by cell growth, division and gene replication. Periodically two events occurs: the replication at a deterministic time $\tau_{R}$ in the cell cycle that doubles the mRNA production rate and the division after a time $\tau_{D}$ starting a new cell cycle. We can then consider the evolution of one type of mRNA and protein through a cell lineage (see Figure 5a).

Overall, the model represents five aspects that intervene in protein production (see Figure 5b):

**mRNAs: **

Messenger-RNAs are transcribed at rate $\lambda_{1}$ before the replication and at rate $2\lambda_{1}$ after gene replication. When transcription happens, the number $M$ of mRNAs increases by $1$. As in the previous model, each mRNA has a lifetime of rate $\sigma_{1}$ (so the global mRNA degradation rate is $\sigma_{1}M$).

**Proteins: **

Each mRNA can be translated into a protein at rate $\lambda_{2}$ (so the global protein production rate is $\lambda_{2}M$). The number of proteins $P$ is increased by $1$. As the protein lifetime is usually of the order of several cell cycles, we consider that the proteins do not degrade.

**Gene replication: **

After a deterministic time $\tau_{R}$ after each division (with $\tau_{R}<\tau_{D}$), the gene replication occurs. The gene responsible for the mRNA transcription is replicated, hence doubling the mRNA transcription rate until next division.

**Division: **

Divisions occur periodically after deterministic time interval $\tau_{D}$. After each division, we only follow one of the two daughter cells (see Figure 5a) so that each mRNA and each protein have a probability $1/2$ to be in this next considered cell. The result is a binomial sampling: for instance, with $M$ mRNAs just before division, the number of mRNAs after division follows a binomial distribution ${\mathcal{B}}\left(n,p\right)$ with parameters $n=M$ and $p=1/2$. Moreover, as there is only one copy of the gene in the newborn cell, the mRNA transcription rate is anew set to $\lambda_{1}$ until the next gene replication.

**Volume growth: **

The volume $V(s)$ considers the entire volume of the cell at any moment $s$ of the cell cycle and it increases as the cell grows. One considers the volume growth of the cell as deterministic. In real life experiments, a bacteria volume globally grows exponentially (see [70]) and approximately doubles its volume at the time of division $\tau_{D}$. As a consequence, in our model, for a time $s$ in the cell cycle, then the volume grows as

[TABLE]

with $V_{0}$ the typical size of a cell at birth.

From this, if $M_{s}$ and $P_{s}$ denote respectively the number of mRNAs and proteins at time $s$, the concentrations are

[TABLE]

Our main results are the analytical expressions for the distribution of the number of mRNAs and the first moments of the number of proteins at any time $s$ of the cell cycle, when the system is at equilibrium (ie after many cell divisions).

Proposition 4.1**.**

At equilibrium, the mRNA number $M_{s}$ at time $s$ in the cell cycle follows a Poisson distribution of parameter

[TABLE]

The proof of this proposition uses the framework of Marked Point Poisson Processes to determine the mRNA distribution at any time of the cell cycle knowing the random variable $M_{0}$, the number of mRNA at birth. Then, the distribution of $M_{0}$ can be determined as the equilibrium stipulates that the distribution of $M$ is the same at time [math] and just after the next division. The proof will be given in [34] (in preparation).

Proposition 4.2**.**

At any time $s\in[0,\tau_{R}[$ before replication, the mean and the variance of $P_{s}$ are given as a function of the moments of $\left(P_{0},M_{0}\right)$ by:

[TABLE]

with $x_{0}$ defined in Proposition 4.1.

Moreover, we can show a similar result for $s\in[\tau_{R},\tau_{D}[$, a time after gene replication. Then, at equilibrium, we have that the distribution of the couple $\left(M,P\right)$ after division is the same as at the beginning of the cell cycle. This allows us to determine the explicit values for $\mathbb{E}\left[P_{0}\right]$, $\mathbb{V}\mathrm{ar}\left[P_{0}\right]$ and $\mathbb{C}\mathrm{ov}\left[P_{0},M_{0}\right]$. These complementary results and the proof of Proposition 4.2 will be available in [34] (in preparation).

4.2. Discussion

In order to have realistic values for the parameters $\lambda_{1}$, $\sigma_{1}$, $\lambda_{2}$ and $\tau_{R}$, we fix them to correspond to the average production of each protein measured in Taniguchi et al. (2010) [68]. It gives groups of parameters that represent a large spectrum of different genes. Then, using Propositions 4.1 and 4.2, as well as the additional results, we can compute the concentration for every gene at any time of the cell cycle.

An example of the evolution of one protein concentration is shown in Figure 6. It appears that the mean concentration at a given time $s$ of the cell cycle $\mathbb{E}\left[P_{s}/V(s)\right]$ is not constant during the cell cycle. The curve of $\mathbb{E}\left[P_{s}/V(s)\right]$ fluctuates around $2\%$ of the global average protein production (denoted by $\mu_{p}$ in the graph). Experimental measures of average protein expression during the cell cycle show similar results: for instance, the expression of genes at different positions on the chromosome have been measured in [69]; it shows a similar profile during the cell cycle and depicts a fluctuation also around $2\%$ of the global average (see Figure 1.d and Figure S6.b of the article).

Overall, for all the genes considered, the effect of the cell cycle (the fluctuation of $\mathbb{E}\left[P_{s}/V(s)\right]$ around $\mu_{p}$) is small compared to the standard deviation at each moment of the cell cycle $\sqrt{\mathbb{V}\mathrm{ar}\left[P_{s}/V(s)\right]}$ (the coloured areas in Figure 6). We have determined the range of parameters where, on the contrary, the fluctuations due to the cell cycle would be higher than the inherent variability of the system. It appears that such regime is possible only for highly active gene (high $\lambda_{1}$) and for mRNAs not very active (low $\lambda_{2}$) and which exist for short periods of time (high $\sigma_{1}$). Biologically, such regime seems unrealistic because the cost of mRNA production would be very high.

In Figure 7, we compare the experimental results obtained in the article of Taniguchi et al. (2010) [68] and our analytical expressions obtained in our model, for the global coefficient of variation (CV) of protein concentration as a function of the average concentration. In the experimental measures, it clearly appears two regimes for the CV of protein concentration: for low expressed proteins (mean protein concentration $<10$), the CV roughly scales inversely with the average concentration; for genes with higher protein production (mean protein concentration $>10$), the CV becomes independent of the average protein production level, the plateau is around $10^{-1}$. In our model, the global tendency of the CV approximately inversely scales the average protein concentration, and there is no lower bound limit, as in the experiments.

We have therefore shown through this model that both DNA replication and division can not explain the second regime observed in Figure 7a. Yet the authors of the experimental study propose another possible origin for the noise observed in highly produced proteins; they propose that it is due to fluctuations in the availability of RNA-polymerases and ribosomes. These two macro-molecules are shared among the productions of all the different proteins because they respectively participate to the transcription of mRNAs and the translation of proteins. Since they are present in limited quantities, fluctuations in their availability could indeed have repercussions on the protein noise. We will present in [34] an extension of the model presented here that takes into account this sharing of ribosomes and RNA-polymerases in the production of all proteins of the bacteria, and the impact it has on the protein noise.

Acknowledgements Bertrand Cloez would like to thank Romain Yvinec for his enthusiastic and stimulating organization of the ”Journées MAS 2016” and Florent Malrieu for many interesting discussions on the subject of his section. The research of Bertrand Cloez was partially supported by the Agence Nationale de la Recherche PIECE 12-JS01-0006-01 and by the Chaire ”Modélisation Mathématique et Biodiversité ” . The research of Aline Marguet was partially supported by the Chaire ”Modélisation Mathématique et Biodiversité”. Renaud Dessalles would like to thank Philippe Robert and Vincent Fromion for their support.