How does variability in cells aging and growth rates influence the malthus parameter?

. The aim of this study is to compare the growth speed of diﬀerent cell populations measured by their Malthus parameter. We focus on both the age-structured and size-structured equations. A ﬁrst population (of reference) is composed of cells all aging or growing at the same rate ¯ v . A second population (with variability) is composed of cells each aging or growing at a rate v drawn according to a non-degenerated distribution ρ with mean ¯ v . In a ﬁrst part, analytical answers – based on the study of an eigenproblem – are provided for the age-structured model. In a second part, numerical answers – based on stochastic simulations – are derived for the size-structured model. It appears numerically that the population with variability proliferates more slowly than the population of reference (for experimentally plausible division rates). The decrease in the Malthus parameter we measure, around 2% for distributions ρ with realistic coeﬃcients of variations around 15-20%, is determinant since it controls the exponential growth of the whole population.


Introduction
Recent biological studies draw attention to the question of variability between cells.We refer to the study of Kiviet et al. published in 2014 [13].A cell in a controlled culture grows at a constant rate v > 0, but this rate can differ from one individual to another.The biological question we address here states as follows.How does individual variability in the growth rate influence the growth speed of the population?The growth speed of the population is measured by the Malthus parameter we define thereafter, also called in the literature fitness.Even if the variability in the growth rate among cells is small, with a distribution of coefficient of variation around 10%, and even if its influence on the Malthus parameter would be still smaller, such an influence may become determinant since it characterises the exponential growth speed of the population.
1.1.1.Paradigmatic age or size equations.Structured models have been successfully used to describe the evolution of a population of cells over the past decades, we refer to Metz and Dickmann [16], the textbook of Perthame [22] and references therein.We focus in this study on classical structuring variables which are age (understood in a broad sense as a physiological age -it may not be the time elapsed since birth) or size.The concentration n(t, a, x) of cells of physiological age a ≥ 0 and size x > 0 at time t ≥ 0 satisfies (1) ∂x g x (a, x)n(t, a, x) + ∂ ∂a g a (a, x)n(t, a, x) + γ(x, a)n(t, a, x) = 0, g a (a = 0, x)n(t, a = 0, x) = 4 ∞ 0 γ(a, 2x)n(t, a, 2x)da, g x (a, x = 0)n(t, a, x = 0) = 0, n(t = 0, a, x) = n in (a, x), in a weak sense.The mechanism at work here can be described as a mass balance.The concentration of cells n(t, a, x) evolves through two terms of transport which stand for the growth and the aging, and with one term of fragmentation: -Transport terms.Both size and age evolve in a deterministic way, the growth speed and aging speed being respectively g x and g a .Let x t and a t be the size and physiological age at time t of a cell born with characteristics (0, x b ) = (a t=0 , x t=0 ).Then their evolution is given by dx t /dt = g x (a t , x t ) and da t /dt = g a (a t , x t ).-Fragmentation terms.We assume that the division is perfectly symmetric: an individual of size x divides into two individuals of size x/2 (and age 0).A cell of current size x and current physiological age a divides with probability γ(a, x)dt between t and t + dt.
The boundary condition for a = 0 ensures that the quantity of new cells of characteristics (0, x) is exactly twice the number of cells of characteristics (s, 2x) for s > 0 that have just divided.The factor 4 is thus the product of a factor 2 arising from the birth of two new cells at each division and another factor 2 coming from the 2x size of dividing cells.
Equation 1 also encloses two paradigmatic equations.1) If g x and γ do not depend on size, integrating (1) in x from zero to infinity, we obtain the age-structured equation.The age-structured equation is a classical equation and we refer to Rubinov [27] and to Perthame [22] (Section 3.9.1)for a complete study.2) If g x and γ do not depend on age, integrating (1) in a from zero to infinity, we obtain the size-structured equation, as introduced by Metz and Dickmann [16].We refer to Mischler and Scher [19] and to references therein for the study of this equation.

Introducing individual variability.
To take into account variability between cells, we extend the previous framework adding an individual feature -an aging rate or a growth rate -as structuring variable.Let V be a compact set of (0, ∞) where the individual feature takes its values and let ρ(v, dv ) be a Markov kernel with support in V × V (i.e.satisfying V ρ(v, dv ) = 1 for all v ∈ V).In this study we focus on two cases of special interest: Model (A+V).The age-structured model with variability, extending the classical setting (see Section 2).The individual feature v we introduce represents here an aging rate.The concentration n(t, a, v) of cells of physiological age a ≥ 0 and aging rate v ∈ V at time t ≥ 0 evolves as    ∂ ∂t n(t, a, v) + ∂ ∂a g a (a, v)n(t, a, v) + γ(a, v)n(t, a, v) = 0, g a (a = 0, v )n(t, a = 0, v ) = 2 S γ(a, v)n(t, a, v)ρ(v, dv )dvda, n(t = 0, a, v) = n in (a, v), in a weak sense, abusing slightly notation here, S being the state space [0, ∞) × V.

Model (S+V).
The size-structured model with variability, as already introduced in Doumic, Hoffmann, Krell and Robert [10] (see Section 3).The concentration n(t, x, v) of cells of size x ≥ 0 and growth rate v ∈ V at time t ≥ 0 evolves as still in a weak sense.
The underlying mechanism of the two previous equations is similar to the one described previously for (1): -Transport terms.The individual feature v does not evolve through time, it is given once and for all at birth.Only the physiological age or the size evolve deterministically through the transport terms in g a or g x .-Fragmentation terms.A cell of current physiological age a or current size x, and feature v, divides with probability γ(a, v)dt or γ(x, v)dt between t and t + dt.It gives birth to two new cells, each of one having age 0 or size x/2, and feature v with probability ρ(v, dv ).
We stress again that the age-structured and size-structured equations were extensively studied but mainly without variability.Here the novelty is the structuring variable v we add to take into account an individual feature.
1.2.1.Mathematical formulation of our problem.The initial biological question can now be reformulated mathematically.As we will see, the Malthus parameter is defined as the dominant eigenvalue of Model (A+V) or (S+V).We denote it by λ γ,ρ to stress the dependence not only on the division rate γ but also on the variability Markov kernel ρ(v, dv ).The long-time behaviour of the solution to Model (A+V) or (S+V) is expected to be n(t, a, v) ≈ e λγ,ρt N γ,ρ (a, v) or n(t, x, v) ≈ e λγ,ρt N γ,ρ (x, v), with N γ,ρ a stationary profile, where we see that the growth speed of the system is governed by λ γ,ρ .Our aim is to compare the growth speed of the two following populations: 1) Population with variability.Each cell has its own aging or growth rate drawn according to a non-degenerated density ρ(dv ) such that V v ρ(dv ) = v (if one neglects heredity in the transmission of the individual feature).This population grows at speed λ γ,ρ .
2) Population of reference (without variability).All cells age or grow at the same rate v.This population grows at speed1 λ γ,v .
In other words, our aim is to compare λ γ,ρ to λ γ,v .1.2.2.Existing studies on the fitness.Before going ahead let us describe existing studies on the variations of the Malthus parameter.We sum up below two studies dealing with the size-structured model (without variability).We also mention a corpus of articles by Clairambault, Gaubert, Lepoutre, Michel and Perthame (see [5,11] and references therein) comparing the eigenvalue of a partial differential equation with time dependent periodic coefficients (birth and death rates) to the eigenvalue of the equation with time-averaged coefficients.They mainly focus on the agestructured system for the cell division cycle.Evolution equations with time periodic coefficients require specific techniques (Floquet's theory).
Asymmetry of the division.The influence of asymmetry in the division is investigated by Michel in [17] for the model described by −σ , with n(t, x = 0) = 0, where the size of each cell evolves linearly (at speed g x ≡ 1) and a cell of size x divides at a rate B(x) into two cells of sizes σx and (1 − σ)x, for some parameter σ ∈ (0, 1).The case of reference is the symmetric division case corresponding to σ = 1/2.If we denote by λ B,σ the Malthus parameter in this model, the aim is to compare λ B,σ to λ B,1/2 .Depending on the form of the division rate B, the asymmetry division is beneficial or not for the growth of the overall cell population.Two special cases are highlighted in [17], 1) qualitatively, if cells divide at high sizes (the support of B is far away from zero), then λ B,σ < λ B,1/2 which means that asymmetry slows the growth speed of the population, 2) qualitatively, if cells divides early (the support of B contains zero and B decreases), then λ B,σ > λ B,1/2 which means that asymmetry creates a gain speeding up the growth of the population.We refer to Theorems 2.1 and 2.2 of [17] for precise statements of the assumptions on B. The method developed in [17] is an interesting approach to investigate our question as regards variability.
Influence of the growth rate.The influence of the individual growth rate is investigated by Calvez et al. in [4] for the model described by ), where the size of each cell evolves exponentially at rate v and assuming that division is symmetric for simplicity.On this very simple model, the Malthus parameter is exactly equal to the common individual growth rate v, hence it seems that the faster the cells grow, the faster the overall cell population grows.This is not the case in general: if the growth speed would be g(x) instead of x, then increasing g by a factor may have the effect of diminishing the Malthus parameter.A plausible example: if around infinity g(x) is equivalent to x ν (up to a constant) with ν < 1 and if B(x) vanishes at infinity, being equivalent to x −γ (up to a constant) with 0 < γ < 1 − ν, then the Malthus parameter vanishes when inflating g by a multiplicative factor.This was proved in [4] (see Theorem 1 and Proposition 1 in Appendix 2, see also Figure 2(a)).Note that in the present study, we inflate the growth speed by a factor that depends on each individual.1.3.Main results and outline.Our main results are summed up in Table 1.Our analysis begins with Model (A+V) in Section 2, neglecting heredity in the transmission of the aging rate by picking ρ(v, dv ) = ρ(v )dv .We analytically prove that the Malthus parameter, well-defined by Theorem 3 (we give a proof for the sake of completeness in Section 7), can increase or decrease when introducing variability in the aging rate, depending on the form of the division rate γ (Theorem 4).We also compute the perturbation at order two of the Malthus parameter λ γ,ρα when ρ α converges in distribution to a Dirac mass as α → 0 (Theorem 5).
The paradigmatic age-structured model still being a toy model for cellular division (at least for E. coli, see Robert et al. [23]), we turn to Model (S+V) in Section 3. A numerical study, based on stochastic simulations, is carried out to give preliminary answers.For the division rate γ chosen as a power law with some lag (see Table 1), for ρ(•) a truncated Gaussian distribution with mean v, we infer that λ γ,ρ < λ γ,v .Admittedly this conclusion agrees with biological wisdom but the strength of our methodology is to quantify such a decrease.We evaluate the magnitude of the decrease around 2% when the variability distribution has a realistic coefficient of variation around 15%-20%.Such a decrease is far from being negligible since the Malthus parameter governs the exponential growth of the whole population.In addition, we observe a monotonous relationship: the Malthus parameter decreases when there is more and more variability in the growth rate.
In the aging rate In the growth rate Variations λγ,ρ λγ,ρ (Figure 1) (Figure 2) Comments Analytic result: -Theorems 4 and 5 Table 1.Variations of the Malthus parameter compared to the reference value when introducing variability between cells, for an experimentally realistic division rate.Note: λ γ,ρ means that λ γ,ρ < λ γ,v for a non-degenerated probability distribution ρ(•) with mean v (truncated Gaussian), and so on.
We provide perspectives in Section 4. Finally Section 5 is devoted to the proofs of our main results of Section 2. Section 6 gives supplementary figures and tables to complete the numerical study of Section 3.

The age-structured model with variability
In this section we study the age-structured model with variability (A+V).Lebowitz and Rubinow [14] and Rotenberg [24] already enriched the age-structured model by an additional feature.In [14], the model is structured by the age and the lifetime, called generation time, which is inherited from the mother cell to the daughter cells.The death or disappearance of the cells occur at a constant rate.In [24], the model is structured by a maturity (our physiological age) and by a maturity velocity (our aging rate).Rotenberg's model is further studied by Mischler, Perthame and Ryzhik [18], establishing the existence of a steady state and the long-time behaviour of the solution.
Rotenberg's model is close to the model we are interested in, but different.From a certain viewpoint, the model introduced by Rotenberg [24] and studied in [18] is more general than Model (A+V): there are distinct death and birth rates and maturity velocity can change at any moment during the life of an individual.However, in our Model (A+V) we allow for a general evolution of the aging rate (through the aging speed g a ) and, most importantly, we allow for heredity in the transmission of the aging rate (through the Markov kernel ρ) in Section 2.1 which follows.
2.1.1.Main assumptions.We use the notation S = [0, ∞)×V for the state space of the physiological age and the aging rate.We require the following two assumptions throughout this section.
Assumption 1 (Aging speed g a and division rate γ).The following conditions are fulfilled: (a) Both (a, v) ; γ(a, v) and (a, v) ; g a (a, v) are continuous.
(d) For any v ∈ V, we have Assumption 2 (Markov kernel ρ).There exists a continuous and bounded ρ : These two assumptions enable us to define the Malthus parameter studying the direct and adjoint eigenproblems below.In Assumption 1, we require (a) since we look for continuous eigenvectors.The requirement (b) on the aging speed g a authorises the case g a (a, v) = v we are interested in (Section 2.2).Note that the non-negativity of g a means that cells can only age, but not rejuvenate.The boundedness of the aging speed is technical.The condition (c) is needed in our proof for the uniqueness of the eigenelements.Finally note that (d) is not very restrictive and it can be read as follows: for any v ∈ V, a ; γ(a,v) ga(a,v) is a hazard rate and the associated density is bounded in a (and in v, which belongs to a compact set of (0, ∞)).Assumption 2 on ρ is rather strong but it simplifies the study of the eigenproblem.
These assumptions have the advantage of leading to a simple proof (in Section 5).The recent results of Mischler and Scher [19] should enable us to weaken it.We stress that our objective here is not the extensive study of the eigenproblem, but lies beyond with the study of the variations of the eigenvalue.
2.1.2.The direct and adjoint eigenproblems.We now precisely define the Malthus parameter.To that end, let us introduce the direct eigenproblem, (2) Let C b V) be the set of functions f : V → R which are bounded and continuous and let C 1 b R + ) be the set of functions f : R + → R which are bounded and continuously differentiable.For f : Theorem 3. Work under Assumptions 1 and 2. There exists a unique solution (λ γ,ρ , N γ,ρ , φ γ,ρ ) to the direct and adjoint eigenproblems (2) and (3) . The unique λ γ,ρ defined in such a way is what we call the Malthus parameter (or fitness) and let us now study its variations (with respect to ρ, for γ fixed).

Influence of variability on the Malthus parameter.
A preliminary remark first: let us consider the case of a constant division rate, γ(a, v) = c > 0 for any (a, v) ∈ S. Since γ is constant, we expect that variability in the aging rate has no influence on the Malthus parameter.When ρ(v, dv ) = ρ(v )dv , one can easily check that λ γ,ρ = c which is independent of ρ.

Model specifications.
From now on we consider Model (A+V) with the aging speed set to (4) g a (a, v) = v which means that the physiological age is proportional to the time elapsed since birth, up to a factor v which may change from an individual to another.
Cells divides at a rate γ(a, v)dt = g a (a, v)B(a)dt = B(a)da, since da = g a (a, v)dt.Thus one has to see the division rate B(a) as a rate per unit of physiological age and γ(a, v) = g a (a, v)B(a) as a rate per unit of time.We mimic here the choice made by S. Taheri-Araghi et al. [29] in a more general model -choice relying on biological evidence.

2.2.2.
Variations.In this section, contrarily to the previous one, we neglect heredity in the transmission of the aging rate assuming that (6) ρ(v, dv ) = ρ(v )dv for some continuous and bounded ρ : In this framework2 , the eigenvectors solution to the eigenproblem (2)-( 3) are explicit and an implicit relation uniquely defines the eigenvalue (see Lemma 7 below, Section 5).
From now on we denote λ γ,ρ by λ B,ρ (to stress Specification ( 5)) and λ B,v stands for λ B,δv where δ v is the Dirac mass at point If necessary, set B(a) = 0 for a ≥ a max , and define One wants to compare 1) the Malthus parameters λ B,ρ solution to controlling the growth speed of the population with variability ; to 2) the Malthus parameters λ B,v solution to controlling the growth speed of the population of reference.
Theorem 4. Consider Model (A+V) with Specifications (4), ( 5) and (6).Assuming in addition that First note that Theorem 4 is valid independently of the density ρ on V (with mean v).The first derivative of the density f B being (for B differentiable), a condition of sign on B − B 2 is nothing but imposing f B to be decreasing (statement (i)) or increasing (statement (ii)) .
One canonical example is the following.Assume that B is constant equal to b > 0 so that γ(a, v) = vb for any (a, v) ∈ S.Then, f B is the exponential density of mean b −1 : it is nonincreasing and statement (i) of Theorem 4 applies.
We give the proof of Theorem 4 now: relying mainly on Jensen's inequality, we obtained a straightforward proof.
Proof of Theorem 4. Let λ B,v be such that and λ B,ρ be such that We know that 1) both H v and H ρ are continuous on [0, ∞), 2) both H v and H ρ are decreasing on [0, ∞), 3) when λ goes to infinity, H v (λ) → 0 and Then, by the intermediate values theorem, we know there exists a unique positive λ B,v such that H v (λ B,v ) = 1 and a unique positive λ B,ρ such that H ρ (λ B,ρ ) = 1.Recall that we want to compare λ B,ρ to λ B,v .We claim that 5) if f B ≤ 0, then for any λ > 0. Together with points 1) to 4), this enables us to conclude that λ B,ρ < λ B,v .
Let us prove 5).For any λ > 0, introduce Then H v and H ρ can be written Jensen's inequality, since V vρ(v)dv = v, and 5) immediately follows.
We now compute the second derivative of h λ , Integrating by parts, recalling that we choose B differentiable, we get for any v ∈ V when f B < 0. This ends the proof for the case f B ≤ 0 and the case f B > 0 is treated in the same way (we need f B to be bounded for the integration by parts).
2.2.3.Small perturbations.The density of the aging rates is thought of as a deformation around an average aging rate v.Given a baseline density ρ : V → [0, ∞) continuous and bounded with mean v, we define, for α ∈ (0, 1], (11) so that Vα ρ α (v)dv = 1 and the mean value is constant, i.e.Vα vρ α (v)dv = v.The density ρ α converges in distribution to the Dirac mass δ v as α → 0. Thus we want to investigate the behaviour of the Malthus parameter λ B,ρα -which is a perturbation of λ B,v -for small value of α.
Theorem 5. Consider Model (A+V) with Specifications (4), ( 5) and ρ(v, dv ) = ρ α (v )dv .Then α ; λ B,ρα is twice differentiable and the Malthus parameter λ B,ρα defined by (9) satisfies with λ B,v defined by (10) and where Note that the perturbation plays at order 2, but not at order 1 since the mean v is preserved (we prove dλ B,ρα dα α=0 = 0 in Proposition 8 below, Section 5).The amplitude of the perturbation at order 2 depends on the baseline aging rate density ρ only through its means v and its variance σ 2 .The proof heavily relies on the explicit expressions of the eigenelements solution to the eigenproblem (2)-(3), available when neglecting heredity in the transmission of the aging rate (see Lemma 7 below, Section 5).Remark 6.If B is differentiable, integrating by parts we get and the variation study of α ; λ B,ρα around α = 0 follows as in the proof of Theorem 4 above.
The class of B for which Theorem 4 applies is unfortunately quite restrictive.In view of applications, a monotonous density f B seems not realistic since we expect to observe a mode in the density of the physiological ages at division (exactly as for lifetimes).One can see the distribution of lifetimes, also called doubling times or generation times, in [21] (Figure 1.D), and in [29] for various medium (Figure S10).Given a more realistic division rate B, in this simple model, one can numerically solve Equations ( 9) and (10) in order to compare λ B,ρ to λ B,v .
A plausible division rate would be of the form B(a) = (a − 1) β 1 {a≥1} .For such a division rate, Theorem 4 does not apply since f B is non-monotonous.Figure 1 shows curves α ; λ B,ρα for different values of β and we see that all these curves are non-increasing, which means that variability in the aging rate slows down the growth speed of the overall cell population, for such division rates B and the tested baseline variability kernel (a truncated Gaussian distribution).

The size-structured model with variability
In this section we study the size-structured model with variability (S+V).A precise definition of the Malthus parameter λ γ,ρ is based on the associated eigenproblem.We refer to Doumic and Gabriel [9] (size-structured model without variability) and references therein for the study of the eigenproblem.We carry out a numerical study of our initial question.
3.1.Stochastic approach.Our aim is to approximate numerically the Malthus parameter and to do so we use a stochastic approach.One could use deterministic methods and approximate the Malthus parameter via a discretisation of the PDE.The stochastic approach we propose enables us to build confidence intervals for λ γ,ρ (in order to check if the variation compared to a reference value is significant or not).In addition, the computational cost of our method should be invariant to the dimension of the model and it could be used as well for more complex models.The growth rate of one cell is drawn according to where u − is the parent of u.We choose g x (x, v) = vx which means the size of each sizes evolves exponentially.Such an assumption is realistic, we refer to the analysis of single cells growth by Schaechter et al. [25] and by Robert et al. [23] (Figure 2).Then the size at time t of cell u is Stochastically, cells divide according to the rate γ(x, v)dt, where x b e vt = x is the size after a time t for a cell born with size x b and growth rate v. Cells divide into two equal parts, ( 12) From the two previous relations, we readily obtain that (ξ u , u ∈ T) is a bifurcating Markov chain (see [12] for a definition) with explicit transition This stochastic modeling and the previous deterministic modeling match.Define the set of living cells at time t, (13) Define the measure n(t, dxdv) as the expectation of the empirical measure at time t over smooth test functions f : where (ξ t u , τ u ) is the size at time t of u together with its growth rate, constant through time.Then n(t, dxdv) satisfies in the weak sense of measures the partial differential equation of Model (S+V).We refer to Theorem 1 of Doumic et al. [10] for such a result.

Stochastic tools.
Relying on a stochastic procedure, we numerically evaluate λ γ,ρ for given division rates γ and growth rate densities ρ.Our estimates of λ γ,ρ are based on trees observed between time 0 and time T .The theoretical framework legitimating our approach is thus a continuous time setting as in Bansaye et al. [2] or Cloez [6].
Estimation of the Malthus parameter.The value of the Malthus parameter is encoded in the time evolution of the tree.To build an estimator one can exploit the asymptotic behaviour of empirical means at time t, as studied in Cloez [6], (14) e −λγ,ρt when t is large, in probability, where c γ,ρ (f ) is a constant and W γ,ρ a non-degenerated random variable.Now define the sum of all sizes of living cells at time t, called biomass in the biological literature.If we observe (ξ ) and (ξ T u , u ∈ ∂T T ) we can define an estimator of the Malthus parameter as ( 16) Approximation (14) ensures that this estimator is close to the true value for large T .Note that a large collection of estimators based on ( 14) can be built, choosing different test functions, and we discuss this in Section 6 (see Figure 3).We see that the estimator ( 16) requires the observation of the whole population at two different large times, chosen for simplicity as T /2 and T .Thus in order to numerically approximate λ B,ρ by ( 16) for given parameters B and ρ one need to simulate continuous time trees up to a large time T .
Simulation of a continuous time tree up to time T .Let us be given a division rate B and a density ρ on V. To simulate a full tree up to a time T , we have to keep the birth and division dates.Set ξ ∅ = x ∅ some given initial value and b ∅ = 0.For any cell u − ∈ T, given its birth time b u − and its size at birth ξ u − , we compute (ζ u − , τ u − ) and its division time d u − in the following way, 1) Draw its growth rate τ u − according to ρ.
2) Draw ξ u given ξ u − = x according to the transition P γ (x, y)dy.
3) Recalling ( 12), compute its lifetime ζ u − by ( 17) Set its division time  Division rate.Mimicking again S. Taheri-Araghi et al. [29], we work with a division rate per unit of time γ chosen as (18) γ(x, v) = vxB(x) and that is why we call B a division rate per unit of size.As previously, from now on we denote the Malthus parameter λ γ,ρ by λ B,ρ to emphasize we work under Specification (18).We choose B of the general form for x 0 ≥ 0 and β ≥ 0. This power form in β with some lag x 0 is inspired by experimental data.For instance we refer to [23], Supplementary Figures S1 and S2, where the division rate is estimated for the age-structured and size-structured models.Adjusting β and x 0 is an easy way to mimic qualitatively all these curves.The estimation is not accurate for large sizes or ages since observations are lacking and we opt not to truncate B making it constant after some threshold but to let it grow to infinity.For such a choice of γ and B, one has to draw ξ u given ξ u − according to (20) 2B(2s)ds 1 y≥x/2 dy, in Step 2) of the previous algorithm 3 .In our numerical tests, we pick x 0 = 1 and β = 2.
Parametrisation of the variability.We suppose here that the growth rate of the daughter cell is independent of its parent, by picking a variability kernel ρ(v, dv ) independent of v. Given a non-degenerated baseline density ρ : V → [0, ∞) with mean v, we define, for α ∈ (0, 1], The coefficient of variation, denoted by CV indexed by the probability distribution considered, is defined as the quotient of its standard error and mean.Then CV ρ = η/v and CV ρα = αCV ρ .Set ϕ(x) = (2π) −1/2 e −x 2 /2 the density of the standard gaussian and Φ(x) = x −∞ ϕ(y)dy its cumulative distribution function.We pick for the baseline density ρ, That is to say, ρ is the truncation on [v min , v max ] = V of a Gaussian distribution with mean v and standard deviation σ η .We set v min = 0, v max = 2 so that v = 1 and we choose σ η = 0.70 so that CV ρ = 50% (using the known formulae of the moments of a truncated Gaussian distribution).Once ρ is fixed, it enables us to define the collection ρ α , α ∈ (0, 1] .For α ∈ (0, 1] the support of the density ρ α is V α = [1 − α, 1 + α] and CV ρα = α/2 for such a choice of ρ.The larger the coefficient of variation CV ρα , the more variability in the growth rate. Estimation of the curve CV ρα ; λ B,ρα .For a given α ∈ (0, 1), we simulate M = 50 continuous time trees up to a large time T , picking x ∅ = 2, v ∅ = 1, with division rate B = B x0=1,β=2 and growth rate density ρ(v, dv ) = ρ α (v )dv .
We obtain a collection ( λ T,m , m = 1, . . ., M ) of M estimators of the Malthus parameter computed according to (16).We denote the mean of these M estimators by λ (α) T .It enable us to obtain a reconstruction of the curve CV ρα ; λ B,ρα .
3 Note that for the choice (19) of B, we can easily do it inverting the cumulative distribution function.
The question now states as follows: what can we say about λ B,ρα compared to λ B,v = v the average growth rate, which is the Malthus parameter of a population without variability in the growth rate?3.2.2.Main result: the Malthus parameter decreases when introducing variability.T , together with a confidence interval for each α, such that it contains 95% of the M estimators.The aim is to get confidence intervals of good precision and to that end we choose large enough times T .The mean numbers of living individuals at time T is at least of magnitude 50 000, see Supplementary Table 2 in Section 6. (Such an amount of data, in a full tree case, would correspond to the observation of the first 16 generations.)We first observe that λ B,ρα is significantly lower than the value of reference V vρ α (v)dv = v = 1, for CV ρα exceeding 10%.In addition, the curve is significantly decreasing as CV ρα increases, i.e. as there is more and more variability in the growth rate.Division rate change.The result is robust when changing the division rate B. In particular, we have explored several couples of parameters (x 0 , β) (recall Parametrisation (19)).The conclusion remains the same: the Malthus parameter significantly decreases when there is more and more variability in the growth rate.Supplementary Table 3 in Section 6 displays the results for x 0 = 1 and β = 8.Asymmetric division.So far we have assumed symmetry in the division.However a slight asymmetry can arise, see Marr et al. [15] or Soifer et al. [26] (Figure 8).Asymmetry is defined as the ratio of the size at birth of one daughter cell over the size at division of its parent.The distribution of this ratio shows a mode at 0.50 and has a coefficient of variation of about 4%.In our numerical study, we amplify asymmetry since we aim at measuring an effect of asymmetry on the curve CV ρα ; λ B,ρα .We immediately extend the model presented in Section 3.1 to allow for asymmetry: we assume that a cell of size x divides into two cells of sizes ux and (1 − u)x with u uniformly drawn on the interval [ε, 1 − ε], independently of everything else, for 0 ≤ ε < 1/2.A dividing cell does not produce an arbitrarily small cell, that is why we chose ε = 0. However we fix ε quite close to 0 for our study.Supplementary Table 4 in Section 6 collects the results for ε = 0.10.Once again there is no significant change compared to Figure 2. Thus asymmetry seems to have no significant impact on the Malthus parameter, in presence of variability.
Linear growth.Consider now that the size at time t of cell u is ξ t u = ξ u + τ u t for t ∈ [b u , d u ).In order to simulate a continuous time tree, we still use the transition defined by (20) in Step 2) of our algorithm and one has only to replace (17) by We still observe a penalisation due to variability in terms of the overall cell population growth (see Supplementary Table 5 in Section 6).For coefficient of variations around 15%-20% in the growth rates, the decrease of the Malthus parameter is estimated at 1-1.5%.
Unit size versus unit time division rate.Notice that the model studied in Doumic et al. [10] corresponds to the choice γ(x, v) = B(x) (there B is a rate per unit of time) instead of γ(x, v) = vxB(x) (here B is a rate per unit of size).This is fundamentally different.As one can see in (20), with the choice γ(x, v) = vxB(x), the size at birth of a cell actually does not depend on the growth rate of its parent whereas it would be the case with the choice γ(x, v) = B(x), since B(2s) vs ds 1 {y≥x/2} , (see [10], Equation ( 11)), and this is the main difference.In order to simulate a continuous time tree up to a given time, at Step 2) of our algorithm, we draw ξ u given ξ u − and τ u (simulated in Step 1)) according to the previous equation (we use a rejection sampling algorithm for this Step 2)).Whatever the specification is, our results concerning the Malthus parameter remain unchanged (see Supplementary Table 6 in Section 6): we observe a decrease when there is more and more variability (the decrease seems slightly higher in the case γ(x, v) = vxB(x) than in the case γ(x, v) = B(x)).
Three main conclusions regarding Model (S+V) are in order.For a unit size division rate B experimentally plausible, when there is variability in the growth rate among cells, 1) The Malthus parameter is lower than the value of reference computed assuming all cells grow at the mean growth rate.
2) The variation is of magnitude 2% for experimentally realistic coefficients of variation in the growth rates distribution, around 15-20%. 3) In addition the Malthus parameter is monotonous: it decreases when there is more and more variability.
These conclusions are robust as argued above changing the unit size B, introducing asymmetry in the division, assuming the individual growth of each cell is linear instead of exponential or even considering a unit time division rate B. We stress that our conclusion 1) coincides with conventional wisdom in biology and our methodology has the advantage of bringing some quantification through 2).

Discussion
The scope of the perspectives is large and includes both theoretical, with analytical and statistical aspects, and experimental issues.We mention here some open questions.
Eigenproblem.In order to prove existence and uniqueness of the eigenelements in Model (A+V), we would like to find minimal assumptions on the aging speed g a , on the division rate γ and on the variability kernel ρ (based on the general results of Mischler and Scher [19] for instance).
Without heredity -ρ(v, v ) = ρ(v ).For Model (A+V), can we build more general classes of B (including experimentally more plausible B) in order to discriminate between the two cases λ B,ρ > λ B,v and λ B,ρ < λ B,v ?For Model (S+V), one can ask if the these two cases are possible.Is there some plausible B such that λ B,ρ > λ B,v ?We would like to build general classes of B to discriminate between the two cases (specifying a density ρ if needed, a truncated Gaussian for instance), and to compute the perturbations of λ B,ρα around λ B,v for ρ α tending to δ v as α → 0 (using the same kind of tools as to prove Theorem 5, see also Michel [17]).
With heredity -general ρ(v, v ).We would like to take into account heredity in the transmission of the aging or growth rate, considering a general Markov kernel ρ(v, v ).In particular, for Model (A+V), a natural question is: how does heredity in the transmission of the aging rate influence the results of Theorems 4 and 5?
Malthus parameter.Beyond estimation, to conduct statistically reliable tests on the Malthus parameter would be of great interest, especially for experimental issues.
Alternative models.Some other models successfully describe the division of E. coli.We refer to Amir [1] and Taheri-Araghi et al. [29].We wonder what can be said on the Malthus parameter in these two models.Preliminary answers are given in Olivier [20] (Chapter 4).

Proof of Theorem 5
As a preliminary, note that neglecting heredity in the transmission of the aging rate enables us to obtain explicit expressions for the eigenvectors N γ,ρ and φ γ,ρ , and an implicit relation which uniquely defines the Malthus parameter λ γ,ρ .Lemma 7. Work under Assumption 1. Assume in addition that ρ(v, dv ) = ρ(v )dv for some continuous and bounded ρ : The unique solution N γ,ρ such that 3) are respectively given by, for any (a, v) ∈ S, with κ, κ normalizing constants such that S N γ,ρ = 1 and S N γ,ρ φ γ,ρ = 1.
For a proof, one can easily check that N γ,ρ and φ γ,ρ defined in such a way satisfy respectively (2) and ( 3).The uniqueness is guaranteed by Theorem 3 (see a proof in the appendix, Section 7).
First note that since ρ α converges in distribution to the Dirac mass at point v as α → 0, we get the convergence of λ B,ρα to λ B,v as α → 0 using the characterisations ( 9) and (10).In a first step, we aim at computing the first derivative of α ; λ B,ρα .
Step 1.Let α ∈ (0, 1] be fixed.For 0 < ε < α, we claim that Indeed, operators A α and A * α are dual, then [17], Lemma 3.2 and Equation (3.11)).Using the definition (23) of A α , we get which we insert in (24) to obtain and we let ε go to zero to study the differentiability on the left of α ; λ α at point α.In the same way we compute ε −1 (λ α+ε − λ α ) to study the differentiability on the right of α ; λ α at point α.
Step 2. The adjoint eigenvector φ α has an explicit expression we exploit now.For β ∈ (0, 1], recalling the explicit expression of φ β given by Lemma 7, with κ β defined in Lemma 7. Thus (25) becomes Before going ahead in computations, recall Specifications (4), g a (a, v) = v, and ( . The previous equality boils down to ( 26) with f B defined by ( 8) and κα = κ α κ α equal to Step 3. In the case of the kernel (11), we can explicitly compute the derivative with respect to α.
After a change of variables (setting α −1 (v − v(1 − α)) as new variable), (26) reads Inverting the derivative in α and the integral and computing the derivative with respect to α, we get the announced result, recalling the definition (27) of κα .
As a corollary of Proposition 8, for all division rate B, when α goes to zero, the first derivative is null, (28) dλ B,ρα we picked a baseline density ρ with mean v).So we compute the second derivative when α converges to zero.
Step 1.We first treat the second term of Gathering (31) and (32) enables us to compute the right-hand side of (30).
Step 2. We now check that the first term of converges to zero as α → 0. One readily checks that Robustness of our results: division rate change.Number of cells versus biomass.Recall Approximation (14).Observing (ξ t u , τ u ), u ∈ ∂T t , or only a component of it for all living cells, at two different times, T /2 and T for instance, with T large enough, one can estimate the Malthus parameter, with a free choice for the smooth test function f .One common choice is f (x, v) = x (the empirical mean defines in this case the mean size of the cells) and it lead us to the estimator (16).Another choice would be f ≡ 1 (the empirical mean defining here the mean number of cells) and it lead us to the estimator where |∂T T | stands for the cardinality of the set (13) of living particles at time T .It is interesting to compare these two estimators.For a fixed T , the estimator ( 16) is better than (33) in the sense that its standard deviation is smaller, as pointed out in Supplementary Figure 3. However the two estimators perform equivalently for large T .Green upper curve: estimation by (33) via the number of cells.

Appendix
Theorem 3 concerns the eigenproblem of the age-structured model with variability.Contrary to Lemma 7, we now work in the general case where ρ is a Markov kernel.
Proof of Theorem 3. We first study the direct eigenproblem (2), then we turn to the adjoint eigenproblem (3).At last we prove uniqueness of the eigenelements.The methodology we use is inspired by the one of [8].In this proof we denote by C b (V) = X the Banach space of bounded continuous functions f : V → R, equipped with the supremum norm, f X = sup v∈V |f (v)|, for V compact set of (0, ∞).
Direct eigenproblem.We split the proof into four steps.
Step 1.Since N γ,ρ satisfies (2), we immediately deduce that, for any (a, v) ∈ S, Using the boundary condition leads us to Note that, for every v ∈ V, a ; f γ/ga (a, v) is a density.Equation (35) leads us to define an operator G λ : X → X by for any λ ≥ 0. In the next steps, we look for a solution (λ, f ) to the equation G λ (f ) = f .We work on the space X of continuous functions since we want to apply the Krein-Rutman theorem [7] (then the interior of the positive cone is the set of positive functions).
Step 2. We introduce a so-called regularized operator, which is strictly positive.For a fixed ε > 0, set where |V| stands for the Lebesgue measure of the compact set V ⊂ (0, ∞) and define the operator for any λ ≥ 0. We claim that G λ,ε is 1) strictly positive on X (i.e. for any f ∈ X non-negative and different from the null function, (G λ,ε f )(v ) > 0 for any v ∈ V), 2) a linear mapping from X into itself, 3) continuous and 4) compact.Thus we are now in position to apply the Krein-Rutman theorem (we use Theorem 6.5 of [22]).For any λ ≥ 0 there exist a unique µ λ,ε > 0 and a unique positive U It just remains to prove the four claimed properties. 1) is precisely achieved thanks to the regularisation ρ ε of ρ by (37).2) The linearity is obvious and for f ∈ X , we have Therefore by the Ascoli-Arzelà theorem for any λ ≥ 0 the family G λ,ε (f ), f ∈ X is compact in X .
Step 4. In this last step, the aim is to let ε tend to zero.Let λ ε be uniquely defined by (39) and denote by U λε,ε = U ε ∈ X the associated positive eigenvector such that U ε X = 1.On the one hand, the family, (U ε , 0 < ε < 1) is compact in X (recall that U ε = G λε,ε U ε and use again the Ascoli-Arzelà theorem).On the other hand, the family (λ ε , 0 < ε < 1) is bounded, recalling (39) and (41).Thus we can extract a subsequence, still denoted by (λ ε , U ε ), converging to ( λ, Ū ) in R × X with λ ≥ 0 and Ū ∈ X positive such that Ū X = 1.Since Adjoint eigenproblem.The proof follows the same steps as in the direct eigenproblem.
Step 4. We let ε go to zero as previously and we find ( λ * , H) with λ * ≥ 0 and H non-negative, φ * X = 1, such that G * λ * H = H, which means we have found a solution to (42).Recalling (43), we set Uniqueness of the eigenelements.We successively prove the uniqueness of the eigenvalue, of the direct eigenvector and of the adjoint eigenvector.
The proof of Theorem 3 is now complete.

3. 1 . 1 .
Stochastic modeling.We closely follow here the description given in Doumic et al.[10].Introduce the infinite genealogical treeT = ∞ m=0 {0, 1} m .Informally we may view T as a population of cells (where the initial individual is denoted by ∅).Each node u ∈ T, identified with a cell of the population, has a mark(b u , ζ u , ξ u , τ u ),where b u is the birth time, ζ u the lifetime, ξ u the size at birth and τ u the individual feature of cell u.We introduce d u the division time of u, d u = b u + ζ u .

3 . 2 .
and do not simulate further the descendants of u − .(ii) If d u − < T , keep (b u , ξ u ) and go back to Step 1).Influence of variability on the Malthus parameter.

3. 2 . 1 .
Numerical protocol.We first specify the division rate γ and the growth rate Markov kernel ρ.

Table 4 . 1 . 6 .
Model (S+V).Division rate γ(x, v) = vxB(x) with B(x) = (x − 1) 2 1 {x≥1} .Asymmetric division (a cell of size x splits into two cells of size ux and (1 − u)x for u uniformly drawn on [0.1, 0.9]).Estimation of the Malthus parameter λ B,ρα (mean and 95% confidence interval based on M = 50 Monte Carlo continuous time trees simulated up to time T ) with respect to the coefficient of variation of the growth rates density ρ α with mean v = 1.Reference (all cells grow at a rate v = 1):λ B,v = Table Model (S+V).Division rate γ(x, v) = B(x) with B(x) = (x − 1) 2 1 {x≥1}.Estimation of the Malthus parameter λ B,ρα (mean and 95% confidence interval based on M = 50 Monte Carlo continuous time trees simulated up to time T ) with respect to the coefficient of variation of the growth rates density ρ α with mean v = 1.Reference (all cells grow at a rate v = 1): λ B,v = 1.

Figure 3 .
Figure 3. Model (S+V).Standard deviation of two estimators of the Malthus parameter as T increases (based on M = 50 Monte Carlo continuous time trees simulated up to time T ), for ρ α=0.3 and division rate γ(x, v) = vxB(x) with B(x) = (x − 1) 2 1 {x≥1} .Blue lower curve: estimation by (16) via the biomass.Green upper curve: estimation by (33) via the number of cells.

Table 3 .
Model (S+V).Division rate γ(x, v) = vxB(x) with B(x) = (x − 1) 8 1 {x≥1} .Estimation of the Malthus parameter λ B,ρα (mean and 95% confidence interval based on M = 50 Monte Carlo continuous time trees simulated up to time T ) with respect to the coefficient of variation of the growth rates density ρ α with mean v = 1.Reference (all cells grow at a rate v = 1): λ B,v = 1.