Adrien Mazoyer

Fluctuation analysis on mutation models with birth-date dependence

published or not. The documents may come

Fluctuation analysis on mutation models with birth-date dependence

Introduction

Mutation models are probabilistic descriptions of the growth of a population of cells, in which scarce mutations randomly occur. Data are samples of integers, interpreted as final numbers of mutant cells. The frequent appearance in the data of very large mutant counts, usually called "jackpots", evidences heavy-tailed probability distributions. Mutation models have two objectives: study the distribution of the number of mutant cells at the end of the growth process; perform fluctuation analysis on data to estimate the probability for a mutant to appear at any division.

Any classic mutation model can be interpreted as the result of the three following ingredients [START_REF] Hamon | Statistics for the Luria-Delbrück distribution[END_REF]:

1. a random number of mutations occurring with small probability among a large number of cell divisions. Due to the law of small numbers, the number of mutations approximately follows a Poisson distribution. The expectation of that distribution is the product of the mutation probability by the total number of divisions;

2. from each mutation, a clone of mutant cells growing during a random time. Due to exponential growth, most mutations occur close to the end of the process, and the developing time of a random clone has exponential distribution. The rate of that distribution is the relative fitness, i.e., the ratio of the growth rate of normal cells to that of mutants;

3. the number of mutant cells that any clone developing for a given time will produce. The distribution of this number depends on the modeling assumptions, in particular the lifetimes of mutants.

This approach leads to a family of probability distributions which depend on the expected number of mutations and the relative fitness. One of the most used mutation models is the well known Luria-Delbrück model [START_REF] Luria | Mutations of bacteria from virus sensitivity to virus resistance[END_REF]. A review on the Luria-Delbrück distribution for the second half of the twentieth century can be found in [START_REF] Zheng | Progress of a half century in the study of the Luria-Delbrück distribution[END_REF]. Here we try to give a historical summary of previous works, including last decades. Mathematical descriptions were introduced by Lea and Coulson [START_REF] Lea | The distribution of the number of mutants in bacterial populations[END_REF], followed by Armitage [START_REF] Armitage | The statistical theory of bacterial populations subject to mutation[END_REF] and Bartlett [START_REF] Bartlett | An Introduction to Stochastic Processes, with Special Reference to Methods and applications[END_REF]. In that model, division times of mutant cells were supposed to be exponentially distributed. Thus a clone develops according to a Yule process [42, p. 35]; [4, p. 109], and its size at any given time follows a geometric distribution. The distribution of final mutant counts is also explicit when lifetimes of mutant cells are supposed to be constant. This latter model is called Haldane model by Sarkar [START_REF] Sarkar | Haldane's solution of the Luria-Delbrück distribution[END_REF]; a first practical algorithm for computing the Haldane distribution was proposed much later by Zheng [START_REF] Zheng | On Haldane's formulation of the Luria-Delbrück mutation model[END_REF]; an explicit form of the asymptotic distribution is finally given by Ycart [START_REF] Ycart | Fluctuation analysis: can estimates be trusted?[END_REF]. General lifetimes have also been studied in [START_REF] Ycart | Fluctuation analysis: can estimates be trusted?[END_REF], but no explicit distribution is available apart from the exponential and constant lifetimes. Other extensions of the Luria-Delbrück model take into account the case where cells have a certain probability to die rather than divide [2, sec. 3.1]; [START_REF] Dewanji | A generalized Luria-Delbrück model[END_REF][START_REF] Komarova | The fixed-size Luria-Delbrück model with a nonzero death rate[END_REF][START_REF] Ycart | Fluctuation analysis with cell deaths[END_REF], where final number of cells are random [START_REF] Angerer | An explicit representation of the Luria-Delbrück distribution[END_REF][START_REF] Komarova | The fixed-size Luria-Delbrück model with a nonzero death rate[END_REF][START_REF] Ycart | Unbiased estimates of mutation rates under fluctuating final counts[END_REF][START_REF] Zheng | A second look at the final number of cells in a fluctuation experiment[END_REF]. As mentioned above, the main statistical objective of mutation models is the estimation of the probability for a mutant cell to appear upon any given cell division. It is computed dividing estimate of the mean number of mutations by the mean final number of cells. Computing robust estimates is of crucial importance in medical applications, like cancer tumor relapse or multidrug resistance of Mycobacterium Tuberculosis for instance. Estimates are realizations of an estimator which is a random variable depending on the considered sample. A robust estimator satisfies two properties: consistency, and explicit asymptotic distribution. Thus confidence intervals and p-values can be computed.

Luria and Delbrück [START_REF] Luria | Mutations of bacteria from virus sensitivity to virus resistance[END_REF] have proposed two estimators. The first estimator, called p 0 estimator, is based on the relation between the probability to get a null count in the sample and the mean number of mutations: taking the negative logarithm of the relative frequency of zeros among the sample gives a robust estimate of the expected number of mutations [20, eq. (5)]. Remark that if the sample does not contain null counts, the method cannot be applied. The second estimator proposed in the same article is based on the relation between the mean number of mutants, the sample size, and the final mutant of cells [20, eq. (8)]. Since this estimator does not have expectation, it is not consistent and should not be used. A wide panel of estimation methods has been proposed since then [START_REF] Rosche | Determining mutation rates in bacterial populations[END_REF][START_REF] Foster | Methods for Determining Spontaneous Mutation Rates[END_REF]. Most of these methods are based on empirical median of the mutant count to reduce the heavy tail effects [18, eq. ( 25)], [13, eq. (6)]. Even if some median methods give good results in practice, the consistency property is not satisfied or cannot be checked: indeed the empirical median is not a robust estimator of the median for discrete distributions. Thus other methods which satisfy the properties of interest should be considered. Since the distribution of final numbers has an explicit form, the Maximum Likelihood (ML) seems to be an obvious optimal choice [START_REF] Ma | Analysis of the Luria-Delbrück distribution using discrete convolution powers[END_REF][START_REF] Stewart | Fluctuation Tests: How Reliable Are the Estimates of Mutation Rates?[END_REF][START_REF] Zheng | New algorithms for Luria-Delbrück fluctuation analysis[END_REF]. The computation of the likelihood and its derivatives can be numerically unstable in the sample contains large jackpots. One of the possible ways to reduce such tail effects is "Winsorization" of the sample [38, sec. 2.2], which consists in replacing any value of the sample that pass a certain bound by the bound itself. However, since very large numbers of mutants are not countable in practice, such cases where the Winsorization is required are not common. The last method exposed here uses the probability generating function [START_REF] Rémillard | Inference based on the empirical probability generating function for mixtures of Poisson distributions[END_REF][START_REF] Hamon | Statistics for the Luria-Delbrück distribution[END_REF], and is called Generating Function (GF) method. This method is comparable in terms of efficiency to the ML method. Moreover, this method has a good numerical stability and a negligible computing time. However, this method depends on tuning parameters. These parameters should be set according to the data, which is not possible in practice. Therefore, the GF method is implemented as a trained semi-parametric method. The p 0 , ML, and GF methods provide asymptotically normal estimators of the mean number of mutation. The estimation of the mutation probability can then be deduced dividing estimation of mean number of mutations by the mean final number of cells. The fluctuations of the final number of cells can also be taken into account [START_REF] Ycart | Unbiased estimates of mutation rates under fluctuating final counts[END_REF], in order to get a more accurate estimation. Sometimes, data are samples of couples of integers, interpreted as final numbers of mutants and final numbers of cells. In that case, the Maximum Likelihood can be used to estimate directly the mutation probability [START_REF] Ycart | Unbiased estimates of mutation rates under fluctuating final counts[END_REF]. Moreover, the relative fitness can also be estimated.

The lifetimes of the cells are supposed to be i.i.d. in the classic mutation models. Thus, the population grows exponentially or dies out [START_REF] Athreya | Branching Processes[END_REF][START_REF] Jagers | Stabilities and instabilities in population dynamics[END_REF]. However, the growth is in practice logistic [START_REF] Laird | Dynamics of tumor growth[END_REF]: it is exponential until an inflexion instant when the growth begins to slow and eventually tends to zero. Indeed, during an experiment, a colony of cells grows in an environment which contains a finite amount of resources. Then a cell born at a instant s 1 will complete its lifetime faster than a cell born at a instant s 2 s 1 . The Verhulst model [START_REF] Verhulst | Notice sur la loi que la population suit dans son accroissement[END_REF] is one of the most known deterministic growth model which takes into account this limitation. Logistic-type stochastic models are described by Allen [1, sec. 9.4.2], and mathematically studied by several authors among which [START_REF] Tan | A stochastic Gompertz birth-death process[END_REF][START_REF] Tan | On stochastic growth processes with application to stochastic logistic growth[END_REF][START_REF] Lambert | The branching process with logistic growth[END_REF]. Stewart et al. [START_REF] Stewart | Fluctuation analysis: the probability distribution of the number of mutants under different conditions[END_REF] proposed an approach to take into account the decreasing rate of division as the cells run out of resources. Houchmandzadeh [START_REF] Houchmandzadeh | General formulation of Luria-Delbrück distribution of the number of mutants[END_REF] described a discrete approach with a general growth model for mutant clones. A Luria-Delbrück model assuming that the replication number of any normal cell is limited whereas mutants are not has been exposed in [START_REF] Rodriguez-Brenes | Cellular replication limits in the Luria-Delbrück mutation model[END_REF]. However, none of these studies provide results for the non-i.i.d. lifetimes case, in particular on the distribution of final mutant count.

In a previous work [START_REF] Mazoyer | Time inhomogeneous mutation models with birth-date dependence[END_REF], an extension of the classic mutation models to the case where the split instant of a cell depends on its birth date has been proposed. The results on the asymptotic distribution of the mutant count were very similar to that of classic Luria-Delbrück model. Therefore the methods of estimation described above should be directly adapted to the model with birth-date dependence. As for the homogeneous case, the three methods provide consistent and asymptotically normal estimator for the parameters of interest. However, fast simulation cannot be deduced from the approach exposed in [START_REF] Mazoyer | Time inhomogeneous mutation models with birth-date dependence[END_REF]. Such algorithms are necessary to perform large scale simulation studies. Another approach of the model is proposed here: as for the homogeneous mutation models, the distribution of the final mutant count can still be interpreted as the result of three ingredients. As a direct consequence, a fast simulation algorithm can be deduced. The asymptotic results on the distribution of the mutant count of [START_REF] Mazoyer | Time inhomogeneous mutation models with birth-date dependence[END_REF] are recovered and extended to the case where the death of normal cells are taken into account.

General modeling assumptions are described in Section 2. The three-ingredients approach is exposed and used to prove the convergence in distribution of the final mutant count in the Section 3. Probability computation and simulation algorithms are exposed in Section 4. The case where the hazard functions associated to the split instant distribution of normal and mutant cells are proportional is studied. In particular, the Luria-Delbrück distribution with cell deaths [START_REF] Ycart | Fluctuation analysis with cell deaths[END_REF] is recovered. The Haldane model is also recovered, and extended to the case where mutant cell deaths are taken into account. The statistical question of estimation of the number of mutations and the relative fitness is studied in Section 5: assuming that the other parameters are known, the methods p 0 , ML, and GF methods can be extended to inhomogeneous models. Estimation biases induced by considering classic mutation models instead of model with birth-date dependence are illustrated with simulation experiments in Section 6. In particular, simulations seem to show that the bias will be in practice negligible, which encourages to continue to use the classic model for estimation.

Hypotheses and models

Notations and hypotheses are described in this section. A rigorous definition of the probabilistic model as a tree-indexed process has already been given in [23, sec. 2]. Thus, the dynamics are shortly described, and the modeling assumptions will be summarized at the end of this section.

Consider a normal cell born at a given instant s. At a random instant (called here a final instant) with cumulative distribution function (cdf) F ν (s, •), the cell produces one normal and one mutant cell with probability π (this event is called a mutation), two normal cells with probability 1 -π -γ, or dies with probability γ. Consider now a mutant cell, born at a given time s. At a random instant with cdf F µ (s, •), the mutant produces two mutant cells with probability 1 -δ or dies with probability δ. Starting from a single cell, whatever its nature, the set of all descendants constitutes a clone. Thus the clone size at a given time t denotes the number of cells alive at time t in the clone. Consider a given cell, the mutation or death events are independent from its final instant. Two cells are independent conditionally on their common ancestor. Therefore, the clones stemming from these cells are also independent conditionally on this ancestor. Remark that those dependence assumptions hold whatever the nature of the considered cells. At the beginning of the process, the population contains a given number n of normal cells and no mutants.

Some details about F ν and F µ are now given. The final instant of a cell born at a given time s 0 cannot be smaller than s. Thus, both cdf satisfies F ν (s, t) = 0 and F µ (s, t) = 0 for t s. Moreover, the total number of cells is in practice bounded by the carrying capacity. It corresponds to the maximum sustainable population: the closer to this bound the number of cells, the slower the growth of the population. In other words, some cells do not produce descendants before the end of the growth process. Thus, the distribution of the final instant of any cell may have a positive mass at infinity. The cdfs F ν (s, •) and F µ (s, •) are cdfs of subprobability measures on R + , i.e. the limit of F ν (s, t) and F µ (s, t) as t tends to infinity may be strictly smaller than 1 for any s in R + . For more details about subprobability measures, see for example [25, p. 170]. Thus, F ν (s, •) and F µ (s, •) are assumed to be cdfs on the extended real line R + = R + ∪ {+∞} for any s ∈ R + . Construction of such cdfs is described in [START_REF] Mazoyer | Time inhomogeneous mutation models with birth-date dependence[END_REF]. Additive assumptions on the cdf F ν (s, •) are now precised. For any s 0, let F (s, •) be a cdf on R + such that F (s, t) = 0 if t s. The cdf F will satisfy (H) if there exists a cdf of subprobability on R + , denoted by F (s, •), such that the following holds: (H 1 ) the cdf F is differentiable with respect to s and t, and decreasing in s;

(H 2 ) lim t→+∞ F (s, t) 1 for all s ∈ R + and F (s, t) = 0 if t s; (H 3 ) for any s 0, F (s, •) is deduced from F (s, •) by F (s, t) = F (s, t)1 t∈[0 ; +∞) + 1 t=+∞ ; (H 4 ) the function h defined for all (s, t) ∈ R 2 + by h(s, t) = -log 1 -F (s, t) ,
satisfies for any t s: h(s, t) = h(0, t) -h(0, s) .

By definition, h is positive, differentiable with respect to s and t, increasing in t, decreasing in s, and for any (s,t) ∈ R 2 + : h(s, t) lim t→+∞ h(0, t) .

There exists a positive, continuous, R + -valued function λ such that:

h(s, t) = t s λ(u)du .
The function h can be interpreted as the cumulative hazard rate on an interval [s ; t] associated to F . The function λ can be interpreted as the instantaneous hazard rate associated to F . The cdf F (s, •) is then defined on R + for any s ∈ R + by

F (s, t) =    1 -exp - t s λ(u)du 1 s t if t < +∞ , 1 if t = +∞ . (1)
Moreover, consider a positive, continuous, and R + -valued function λ. Then, the cdf deduced from λ by ( 1) satisfies (H). In particular, if λ(t) tends to 0 as t tends to infinity, then the limit of h(0, t) is finite. In that case, the limit of F (s, t) is smaller than 1 for any s 0. Notice also that if T (s) is a random variable with cdf F (s, •), such that F satisfies (H), then:

P [T (s) > u + t|T (s) > t] = e -h(s,u+t) e -h(s,u+t) = P [T (t) > u + t] .
This is quite similar to the memorylessness property of exponential distributions.

Condition (H 4 ) is quite restrictive. However, the mean growth of the population can be adjusted to any positive, continuous, increasing and R + -valued function [START_REF] Mazoyer | Time inhomogeneous mutation models with birth-date dependence[END_REF]. Remark that condition (H 4 ) is equivalent to satisfy:

F (s, t) = F (0, t) -F (0, s) 1 -F (0, s) .
Therefore, for any R + -valued cdf G, the cdf F defined by

F (s, t) = G(t) -G(s) 1 -G(s) ,
satisfies (H 4 ). From now on, F ν is assumed to satisfy (H). Its related function defined in (H 4 ) will be denoted by h ν . The limit of h ν (0, t) as t tends to infinity will be denoted by h ν,∞ . The instantaneous division rate associate to F ν will be denoted by λ ν . No additive assumptions on F µ are required yet. The modeling assumptions can then be summarized as follows:

• at time 0, n normal cells are present;

• the final instant of any cell depends on its nature and its birth date;

• the final instant of a normal cell born at time s is a random variable with cdf F ν (s, •)

which satisfies (H);

• upon completion of the lifetime of a normal cell:

with probability π one normal and one mutant cells are produced;

with probability γ the cell dies out;

with probability 1 -π -γ two normal cells are produced;

• the final instant of a mutant cell born at time s is a random variable with cdf F µ (s, •) defined on R + ;

• upon completion of the lifetime of a mutant cell:

with probability δ the cell dies out;

with probability 1 -δ two mutant cells are produced;

• for any cell, the events of death or mutation do not depend on its final instant;

• two cells, whatever their nature, are independent conditionally on their common ancestor;

• two clones are independent conditionally on the common ancestor of the two cells which started those clones.

"Three-ingredients" approach

The generalization of the "three-ingredients" interpretation exposed in the introduction is described in this section. Let (τ n ) n∈N be a sequence of observation instants, tending to infinity as n tends to infinity. Let (π n ) n∈N be a sequence of mutation probabilities, tending to 0 as n tends to infinity. Moreover, assume

lim n→+∞ π n ne hν (0,τn) = α ,
where α is some fixed positive real number. Remark that the constant α corresponds in the classic case to the mean number of mutations. Denote by Z n (t) the number of mutations in an interval [0 ; t], starting with n normal cells at time 0. The increasing sequence of mutation instants will be denoted by

T (n) i i∈N
.

Similarly to [START_REF] Hamon | Statistics for the Luria-Delbrück distribution[END_REF], a three-ingredients approach of model mutations can be given for the non-homogeneous birth date case.

Proposition 3.1. Assume γ = 0. Let π = (π n ) n∈N and τ = (τ n ) n∈N two sequences, and α a positive real such that:

lim n→+∞ π n = 0 , lim n→+∞ τ n = +∞ , lim n→+∞ π n ne hν (0,τn) = α .
Then:

(A (0) 1 )
As n tends to infinity, the distribution of the total number of mutations Z n (τ n ) tends to the Poisson distribution with parameter

m = α 1 -e -hν,∞ , where h ν,∞ = lim t→+∞ h ν (0, t) . (A (0)
2 ) As n tends to infinity, the joint distribution of the vector

T (n) 1 , . . . , T (n) k
of k fixed number of mutation instants in an interval [0 ; t] converges to the order statistics of a k sample of the distribution:

λ ν (u)e -hν (u,t) 1 -e -hν (0,t) 1 u∈[0 ; t] ,
i.e. the distribution λ ν (u)e -hν (u,t) truncated on [0 ; t].

In particular, from assertions (A (0)

2 ), the probability generating function (pgf) of the size at a given time t of any mutant clone is given by

I(z, t) = t 0 ψ(z, u, t) λ ν (u)e -hν (u,t) 1 -e -hν (0,t) du , (2)
where ψ(z, s, t) is the pgf of the size at time t of a clone stemming from a mutant born at time s.

Proposition 3.1.

Assertion (A (0) 1 ) Consider the binary branching process with a single initial cell and without mutations. Denote by N 1 (t) the total number of cells at instant t. For any t 0, consider the sequence (N n (t)) n∈N defined for any n > 0 by

N n (t) = n i=1 N (i) 1 (t) , (3)
where the N (i) 1 (t) are i.i.d. copies of N 1 (t). For any n > 0, N n (t) denotes the total number of cells living at time t in n independent copies of N 1 (t). According to Proposition 3.2. of [START_REF] Mazoyer | Time inhomogeneous mutation models with birth-date dependence[END_REF]:

E[N 1 (t)] = e hν (0,t) .
Let ε > 0. According to the law of large numbers, there exists for any t > 0 an integer n 0 (t) such that for any n > n 0 (t):

P N n (t) ne hν (0,t) -1 < 1 -e -hν (0,t) ε = 1 . (4)
Since:

N n (t) -n n (e hν (0,t) -1) -1 = 1 1 -e -hν (0,t) N n (t) ne hν (0,t) -e -hν (0,t) -1 = 1 1 -e -hν (0,t) N n (t) ne hν (0,t) -1 ,
for any n > n 0 (t):

P N n (t) -n n (e hν (0,t) -1) -1 < ε = 1 . (5)
Thus, the number of cell divisions during the time interval [0 ; t] is almost surely equivalent to n e hν (0,t) -1 . Let ϑ > 0. Since τ n tend to infinity, there exists n 1 ∈ N such that for any n > n 1 :

τ n > ϑ . (6)
Therefore, according to ( 5) and ( 6), for any n > max(n 0 (ϑ), n 1 ):

P N n (τ n ) -n n (e hν (0,τn) -1) -1 < ε = 1 . (7)
Thus, the total number of divisions N n (τ n )-n is almost surely equivalent to n e hν (0,τn) -1 . Since the mutant clones develop according to a different dynamic, the number of divisions of normal cells in a mutation model does not have the same distribution. However, since mutations are rare, the difference remains negligible: start with the n independent copies, and mark independently the division as potentially mutant with probability π n . Denote by (X n (t)) n∈N the sequence of the number of marked divisions at time t in the n copies. Since π n ne hν (0,τn) tends to α as n tends to infinity, the distribution of X n (τ n ) converges to the Poisson distribution with parameter m = α(1 -e -hν,∞ ). Thus the number of marked divisions remains bounded in probability. Consider now that the clones stemming from the marked divisions are mutant clones. If a division is marked, all the marked divisions occurring in the clone stemming from it will be ignored. In other words, the number of mutations Z n (τ n ) is smaller than the number of marked divisions X n (τ n ) with probability 1. Moreover, since X n (τ n ) is bounded in probability, the difference between Z n (τ n ) and X n (τ n ) is also bounded. Since this difference is bounded, the number of mutations occasions is equivalent in probability to n e hν (0,τn) -1 . In other words, for any ε > 0, there exists n 2 ∈ N such that for any n > n 2 :

P [|X n (τ n ) -Z n (τ n )| ε] = 0 .
Then the distribution of Z n (τ n ) tends to the Poisson distribution with parameter m. Hence (A

(0) 1 ). Assertion (A (0)
2 ) From ( 5), the number of mutation occasions in an interval [0 ; t] is equivalent in probability to n e hν (0,t) -1 . Then for any k ∈ N and any t 0:

lim n→+∞ P [Z n (t) = k] ι n (k, t) = 1 , with ι n (k, t) = n e hν (0,t) -1 k π k n (1 -π n ) n(e hν (0,t) -1) -k ,
where x is the only relative integer which satisfies for any x ∈ R:

x x < x + 1 ,
i.e. the integer part of x. Hence:

ι n (k, t) = n e hν (0,t) -1 n e hν (0,t) -1 -1 . . . n e hν (0,t) -1 -k + 1 k! × π n 1 -π n k exp n e hν (0,t) -1 log(1 -π n ) ∼ n→+∞ π n n e hν (0,t) -1 k k! exp -π n n e hν (0,t) -1 ∼ n→+∞ π n n e hν (0,t) -1 k k! exp -π n n e hν (0,t) -1 .
Moreover, by construction of {X n (t)} t 0 , the process {Z n (t)} t 0 satisfies the following properties:

1. {Z n (t)} t 0 is simple, i.e. for any t ∈ R + :

lim ∆t→0 1 ∆t P [Z n (t + ∆t) -Z n (t) > 1] = 0 ;
2. {Z n (t)} t 0 has independent increments; 3. Z n (0) = 0 with probability 1.

Denote by ξ n the intensity of the process {Z n (t)} t 0 defined for any t ∈ R + by

ξ n (t) = lim ∆t→0 1 ∆t P [Z n (t + ∆t) -Z n (t) = 1]
.

For any t = (t 1 , . . . , t k ) ∈ R k + , the event T (n) 1 = t 1 , . . . , T (n) k = t k will be denoted by T (n) = t For any k ∈ N and t ∈ R k + and conditionally to T (n) = t, the distribu- tion of the (k + 1)-th mutation instant T (n) k+1 is given by Proposition A.2 f T (n) k+1 | T (n) =t (t) = ξ n (t) t t k ξ n (u)du 1 0<t 1 <•••<t k <t .
The joint distribution of the first k mutation instants T (n) is also given by Proposition A.2:

f T (n) (t) = k i=1 ξ n (t i ) exp - k i=1 t i t i-1 ξ n (u)du 1 0<t 1 <•••<t k .
Consider now the inhomogeneous Poisson process {Y n (t)} t 0 with expectation:

m n (t) = nπ n e hν (0,t) -1 .
For any k ∈ N, denote by

S (n) = S (n) 1 , . . . , S (n) k
the k first occurrence instants of the process {Y n (t)} t 0 . Then, for any k ∈ N, t ∈ R + and t ∈ R k + , the three following assertions hold:

lim n→+∞ P [Z n (t) = k] P [Y n (t) = k] = 1 , (8)
and:

lim n→+∞ f T (n) (t) f S (n) (t) = 1 , (9)
and:

lim n→+∞ f T (n) k+1 | T (n) =t (t) f S (n) k+1 | S (n) =t (t) = 1 . (10)
Thus, for any k ∈ N, t ∈ R + , and conditionally to Z n (t) = k the distribution of the vector

T (n) is given for any t ∈ R k + by f (T (n) | Zn(t)=k) (t) = f T (n) (t)P Z n (t) = k T (n) = t P [Z n (t) = k] = f T (n) (t)P T (n) k+1 > t T (n) = t P [Z n (t) = k] .
Thus, according to ( 8), ( 9) and ( 10):

lim n→+∞ f (T (n) | Zn(t)=k) (t) f (S (n) | Yn(t)=k) (t) = 1 ,
for any k ∈ N. Conditionally to Y n (t) = k, the vector S (n) is distributed as the order statistics of k sample of the distribution (Proposition A.1):

m n (u) m n (t) 1 u∈[0 ; t] = λ ν (u)e hν (0,u) e hν (0,t)-1 1 u∈[0 ; t] = λ ν (u)e -hν (u,t) 1 -e -hν (0,t) 1 u∈[0 ; t] . Hence (A (0)
2 ). Remark that (A (0)

2 ) considers mutation instants, while similar results for homogeneous cases consider developing times of mutant clones. For example, according to Theorems 2.1 and 3.1 of [START_REF] Kuzcek | Almost sure limit results for the supercritical Bellman-Harris process[END_REF], the joint distribution of developing times of k fixed mutant clones converges to the product of k independent exponential distributions.

The asymptotic pgf of the total mutant counts can be explicited, as long as the pgf ψ of a clone size is known. Theorem 3.1. Assume γ = 0. Let π = (π n ) n∈N and τ = (τ n ) n∈N two sequences, and α a positive real such that:

lim n→+∞ π n = 0 , lim n→+∞ τ n = +∞ , lim n→+∞ π n ne hν (0,τn) = α .
As n tends to infinity, the pgf of the number of mutants at time τ n starting with n normal cells tends to the pgf:

φ(z) = exp {-m (1 -I(z))} , (11)
where

I(z) = lim t→+∞ I(z, t) = 1 1 -e -hν,∞ lim t→+∞ t 0 ψ(z, u, t)λ ν (u)e -hν (u,t) du . (12)
A first probabilistic interpretation of Theorem 3.1 is the generalization of the "threeingredients" description exposed in the introduction as follows:

1. a random number of mutations occurring with small probability among a large number of cell divisions. According to (A (0) 1 ), the number of mutations approximately follows a Poisson distribution with expectation m = α 1 -e -hν,∞ ; 2. each mutation appears at a random instant. According to (A (0)

2 ), any mutation instant asymptotically follows the distribution λ ν (u)e -hν (u,t) truncated on the interval [0 ; t]; 3. a mutant clone started at instant s develops according to a given process. Its size at time t follows the distribution with pgf ψ(z, s, t).

Remark that the second ingredient concerns the instants at which a given number of mutant clones are started, instead of their developing times. Observe that Theorem 3.1 holds whether if F µ satisfies (H) or not. The Haldane model can be considered as an example: the final instants of the normal cell are exponentially distributed with rate λ and the lifetimes of the mutant cells are equal to a constant a. Then the cdf F µ (s, •) is defined for any t s by

F µ (s, t) = 1 if t s + a , 0 else ,
and does not satisfy (H). Computation of pgfs and probabilites for Haldane model have already been done in [START_REF] Mazoyer | Time inhomogeneous mutation models with birth-date dependence[END_REF]. The pgf I is given by

I(z) = i 0 b i (z)e -λia 1 -e -λa . (13)
For any i 0, the pgf b i represents the size of a clone stemming from a mutant born at time s in the interval [s + ia ; s + (i + 1)a]. In other words, b 0 (z) = z, and for all i > 0:

b i (z) = δ + (1 -δ) (b i-1 (z)) 2 . ( 14
)
The probabilities of the final mutant count will be recalled in next section. In particular, the results obtained in [START_REF] Ycart | Fluctuation analysis: can estimates be trusted?[END_REF] are recovered when δ = 0 and a = log(2). Proposition 3.1 can be extended to the case γ 0.

Proposition 3.2. Let π = (π n ) n∈N et τ = (τ n ) n∈N
two sequences, and α a positive real such that:

lim n→+∞ π n = 0 , lim n→+∞ τ n = +∞ , lim n→+∞ π n nω(τ n )e h * ν (0,τn) = α , where h * ν (s, t) = (1 -2γ)h ν (s, t) ,
and

ω(t) = 1 -2γ 1 -γ -γe -h * ν (0,t)
, the probability that a clone stemming from a normal cell born at time 0 is not died out at time t. Then:

(A (γ) 1 )
As n tends to infinity, the distribution of the total number of mutations Z n (τ n ) tends to the Poisson distribution parameter:

m = α 1 -ω ∞ e h * ν,∞ -1 , where h * ν,∞ = lim t→+∞ h * ν (0, t) ,
and

ω ∞ = lim t→+∞ ω(t) ; (A (γ)
2 ) As n tends to infinity, the joint distribution of the vector

T (n) 1 , . . . , T (n) k
of k mutation instants in an interval [0 ; t] converges to that of the order statistics of a k sample of the distribution

(ω (u) + λ ν (u)(1 -2γ)ω(u)) e -h * ν (u,t) ω(t) -e -h * ν (0,t) 1 u∈[0 ; t] . Proposition 3.2. Assertion (A (γ)
1 ) Consider again the binary branching process with a single initial cell and without mutations: each cell can die out with probability γ instead of dividing into two cells. Denote by N 1 (t) the total number of cells at instant t. For any t 0, consider the sequence (N n (t)) n∈N defined by [START_REF] Armitage | The statistical theory of bacterial populations subject to mutation[END_REF]. For any t > 0, each copy of N 1 (t) may survive until t with probability ω(t). By the same reasoning as for assertion (A (0) 1 ), the number of divisions in the surviving clones occurring in [0 ; t] is almost surely equivalent to n ω(t)e h * ν (0,t) -1 . Moreover, the number of divisions in the proportion 1 -ω(t) of dying clones is bounded and can be neglected. Thus, the total number of divisions in the n copies is almost surely equivalent to n ω(τ n )e h * ν (0,τn) -1 . Let us mark independently and with probability π n the cells divisions. Consider (X n (t)) n∈N the sequence of marked divisions at time t in the n copies. By the same reasoning as for assertion (A

(0) 1 ), the distribution of the number of mutations Z n (τ n ) tends to the Poisson distribution with parameter m = α 1 -ω ∞ e h * ν,∞ -1 . Assertion (A (γ)
2 ) The number of mutation occasions in an interval [0 ; t] is equivalent in probability to n ω(t)e h * ν (0,t) -1 . Therefore:

lim n→+∞ P [Z n (t) = k] ι n (k, t) = 1 ,
where

ι n (k, t) = n ω(t)e h * ν (0,t) -1 k π k n (1 -π n ) n(ω(t)e h * ν (0,t) -1) -k .
Moreover:

ι n (k, t) ∼ n→+∞ π n n ω(t)e h * ν (0,t) -1 k k! exp -π n n ω(t)e h * ν (0,t) -1 ∼ n→+∞ π n n ω(t)e h * ν (0,t) -1 k k! exp -π n n ω(t)e h * ν (0,t) -1 .
For any t 0 and as n tends to infinity, the distribution of Z n (t) is equivalent to the Poisson distribution with parameter:

m n (t) = nπ n ω(t)e h * ν (0,t) -1 . Assertion (A (γ)
2 ) is then deduced by the same reasoning as for assertion (A (0)

2 ). Therefore, Theorem 3.1 can also be extended to the case where γ 0. Theorem 3.2. Let π = (π n ) n∈N and τ = (τ n ) n∈N two sequences, and α a positive real such that:

lim n→+∞ π n = 0 , lim n→+∞ τ n = +∞ , lim n→+∞ π n nω(τ n )e h * ν (0,τn) = α
As n tends to infinity, the pgf of the number of mutants at time τ n starting with n normal cells tends to the pgf [START_REF] Houchmandzadeh | General formulation of Luria-Delbrück distribution of the number of mutants[END_REF] with:

m = α 1 -ω ∞ e h * ν,∞ -1 ,
and

I(z) = 1 ω ∞ -e -h * ν,∞ lim t→+∞ t 0 ψ(z, u, t) (ω (u) + λ ν (u)(1 -2γ)ω(u)) e -h * ν (u,t) du . ( 15
)
From now on, assume that F µ satisfies (H) and denote by h µ its related function defined in (H 4 ). There exists a positive, continuous, R + -valued function λ µ such that:

h µ (s, t) = t s λ µ (u)du .
From Proposition 3.2 of [START_REF] Mazoyer | Time inhomogeneous mutation models with birth-date dependence[END_REF], the pgf of the size at time t of a clone stemming from a mutant cell born at time s is given by

ψ(z, s, t) = δ(1 -z) + e -h * µ (s,t) ((1 -δ)z -δ) (1 -δ)(1 -z) + e -h * µ (s,t) ((1 -δ)z -δ) , (16)
where, for any (s, t) ∈ R 2 + : h * µ (s, t) = (1 -2δ)h µ (s, t) . Assume also there exists ρ > 0 such that for any t 0:

(1 -2γ)λ ν (t) = ρ(1 -2δ)λ µ (t) . ( 17
)
The constant ρ can be interpreted as the instantaneous ratio of hazard functions λ ν and λ µ .

The assumption of proportional hazard functions is not new: in survival analysis, it is known as the Cox proportional-hazard regression model, which is widely used. Observe that, in the case where λ µ is a positive constant, the constant ρ corresponds to the relative fitness in Luria-Delbrück model. This designation will be kept from now on. Under the asymptotic context of Thereom 3.2, the distribution of final mutant count depends on the mean number of mutation m, the fitness ρ, the death parameters γ and δ, and the limit of h µ (0, t).

Theorem 3.3. Let π = (π n ) n∈N and τ = (τ n ) n∈N two sequences, and α a positive real such that:

lim n→+∞ π n = 0 , lim n→+∞ τ n = +∞ , lim n→+∞ π n nω(τ n )e h * ν (0,τn) = α .
Under [START_REF] Lambert | The branching process with logistic growth[END_REF] and as n tends to infinity, the pgf of the number of mutants at time τ n starting with n normal cells tends to the pgf [START_REF] Houchmandzadeh | General formulation of Luria-Delbrück distribution of the number of mutants[END_REF] with:

m = α 1 -ω ∞ e ρh * µ,∞ -1 ,
and

I(z) = (1 -2γ)ρ ω ∞ -e -ρh * µ,∞ 1 e -h * µ,∞ δ(1 -z) + v((1 -δ)z -δ) (1 -δ)(1 -z) + v((1 -δ)z -δ) Q(v)v ρ-1 dv , (18)
where

h µ,∞ = lim t→+∞ h µ (0, t) , h * µ,∞ = (1 -2δ)h µ,∞ , and
Q(v) = 1 -γ -2γve -ρh * µ,∞ 1 -γ -γve -ρh * µ,∞ 2 . ( 19
)
In particular, when γ = 0, (18) becomes:

I(z) = 1 1 -e -ρh * µ,∞ 1 e -h * µ,∞ δ(1 -z) + v((1 -δ)z -δ) (1 -δ)(1 -z) + v((1 -δ)z -δ) v ρ-1 dv , (20)
which corresponds to Theorem 5.1 of [START_REF] Mazoyer | Time inhomogeneous mutation models with birth-date dependence[END_REF].

Theorem 3.3. Plugging ( 16) in [START_REF] Kuzcek | Almost sure limit results for the supercritical Bellman-Harris process[END_REF] and changing e -hµ(s,t) into v leads to

I(z) = (1 -2γ)ρ ω ∞ -e -ρh * µ,∞ lim t→+∞ 1 e -h * µ (0,t) δ(1 -z) + v(1 -δ)z -δ) (1 -δ)(1 -z) + v(1 -δ)z -δ) 1 -γ -2γve -ρh * µ (0,t) 1 -γ -γve -ρh * µ (0,t) 2 v ρ-1 dv .
Since for any v ∈ [0 ; 1] and t 0:

1 -γ -2γve -ρh * µ (0,t) 1 -γ -γve -ρh * µ (0,t) 2 1 1 -γ ,
applying the dominated convergence theorem leads to the result.

A first example consists in considering a non-negative and increasing function f on R + , with finite limit f ∞ as t tends to infinity. Let h µ be defined for (s, t) in R 2 + by

h µ (s, t) = log f (t) f (s) . ( 21
)
Then the expected size of a mutant clone started at time s is

f (t) f (s) 1-2δ
. In other words, it is possible to fit the average trajectory of the development of the clones to any appropriate function of time defining h µ as [START_REF] Ma | Analysis of the Luria-Delbrück distribution using discrete convolution powers[END_REF]. Moreover, only the ratio of f ∞ over f (0) has an influence on I(z).

As another particular case, if h µ,∞ is infinite, the Luria-Delbrück with cells deaths distribution [START_REF] Ycart | Fluctuation analysis with cell deaths[END_REF] is recovered.

Corollary 3.1. Let π = (π n ) n∈N and τ = (τ n ) n∈N two sequences, and α a positive real such that:

lim n→+∞ π n = 0 , lim n→+∞ τ n = +∞ , lim n→+∞ π n nω(τ n )e h * ν (0,τn) = α .
Assume h µ,∞ = +∞. Under (17) and as n tends to infinity, the pgf of the number of mutants at time τ n starting with n normal cells tends to the pgf [START_REF] Houchmandzadeh | General formulation of Luria-Delbrück distribution of the number of mutants[END_REF] with m = α and:

I(z) = ρ 1 0 δ(1 -z) + v((1 -δ)z -δ) (1 -δ)(1 -z) + v((1 -δ)z -δ) v ρ-1 dv . ( 22
)
In other words, the Luria-Delbrück with cells deaths distribution can be extended to the case where F ν (s, •) and F µ (s, •) are non-exponential distributions, as long as they satisfy (H), are cdfs of true measures on R + , and such that the associated functions λ ν and λ µ are proportional. Moreover, the asymptotic distribution of the number of mutants does not depend on the γ parameter in this specific case. That was already notable in [START_REF] Ycart | Fluctuation analysis with cell deaths[END_REF]: the death parameter of normal cells (denoted by q in this reference) does not appear anymore when taking the asymptotic. From now on, the different mutation models will be referred to from the notations of Table 1.

IM M (m, γ, δ, F ν , F µ )
Inhomogeneous mutation models where the final instant of a normal (resp. mutant) cell born at time s has cdf F ν (s, •)

(resp. F µ (s, •)) (Theorem 3.2)
ILD(m, ρ, γ, δ, h µ,∞ ) IM M models where F µ satisfies (H) under ( 17) (Theorem 3.3)

LD(m, ρ, δ) ILD models where h µ,∞ = +∞ (Corollary 3.1)

H(m, ρ, δ) Haldane model
Table 1: Mutation models designations. The parameters m, ρ, γ, δ, and h µ,∞ respectively, denote the mean number of mutation, the fitness parameter, the probability of dying for a normal cell, the probability of dying for a mutant cell, the asymptotic cumulative division rate of mutants. F ν and F µ , respectively, denote the cdf of the final instant for normal cells, the cdf of final instant for mutant.

Probability calculations

Computation and simulation algorithms for the distribution of the final mutant counts are described here. Consider first the case of IM M models. The pgf (15) can be written as

I(z, t) = k 0 r k (t)z k ,
where r k is defined for any k ∈ N and t ∈ R + by

r k (t) = 1 ω(t) -e -h * ν (0,t) t 0 p k (u, t) (ω (u) + λ ν (u)(1 -2γ)ω(u)) e -h * ν (u,t) du .
and the p k (s, t)'s are probabilities of the size at time t of a mutant clone started at time s. These probabilities mainly depend on the model assumptions. For example, if the cdf F µ satisfies (H), the p k (s, t)'s are given by Proposition 3.3 of [START_REF] Mazoyer | Time inhomogeneous mutation models with birth-date dependence[END_REF]

p 0 (s, t) = δ(1 -e -h * µ (s,t) ) 1 -δ -δe -h * µ (s,t)
, and for k > 0:

p k (s, t) = (1 -p 0 (s, t))P (s, t)(1 -P (s, t)) k-1 ,
where:

P (s, t) = (1 -2δ)e -h * µ (s,t)
1 -δ -δe -h * µ (s,t) . Thus:

r 0 (t) = 1 ω(t) -e -h * ν (0,t) t 0 δ 1 -e -h * µ (u,t) 1 -δ -δe -h * µ (u,t) (ω (u) + λ ν (u)(1 -2γ)ω(u)) e -h * ν (u,t) du ,
and for all k > 0:

r k (t) = 1 ω(t) -e -h * ν (0,t) t 0 1 - δ 1 -e -h * µ (u,t) 1 -δ -δe -h * µ (u,t) P (u, t)(1 -P (u, t)) k-1 (ω (u) + λ ν (u)(1 -2γ)ω(u)) e -h * ν (u,t) du ,
Then, (15) can be given by

I(z) = k 0 r k z k ,
where for any k ∈ N, r k is the limit of r k (t) as t tends to infinity. In particular, the r k 's for ILD models are given by

r 0 = (1 -2γ)ρ ω ∞ -e -ρh * µ,∞ 1 e -h * µ,∞ δ -δv 1 -δ -δv v ρ-1 Q(v)dv ,
and for all k > 0:

r k = (1 -2γ)ρ(1 -2δ) 2 (1 -δ) k-1 ω ∞ -e -ρh * µ,∞ 1 e -h * µ,∞ (1 -v) k-1 (1 -δ -δv) k+1 Q(v)v ρ dv , (23)
where Q is given by [START_REF] Lehmann | Theory of Point Estimation[END_REF]. Remark that for γ = δ = 0, (23) leads to:

r k = ρ 1 -e -ρhµ,∞ 1 e -hµ,∞ (1 -w) k-1 w ρ dw = ρ 1 -e -ρhµ,∞ (B(ρ + 1, k) -B e -hµ,∞ (ρ + 1, k)) ,
where B x (θ, β) is the incomplete Beta function defined for any x ∈ (0 ; 1) by

B x (θ, β) = x 0 w θ-1 (1 -w) β-1 dw ,
and B(θ, β) is the complete Beta function, i.e. B(θ, β) = B 1 (θ, β).

As an other example, consider now H models. The pgf I is then given by [START_REF] Jones | Luria-Delbrück Fluctuation Experiments: Design and Analysis[END_REF]. For any i ∈ N, let us denote by p

(i) k k∈N
the probabilities associated to the pgf [START_REF] Komarova | The fixed-size Luria-Delbrück model with a nonzero death rate[END_REF]. Then pgf I can be rewritting as

I(z) = i 0 e -λia 1 -e -λa k 0 p (i) k z k = k 0 z k i 0 e -λia 1 -e -λa p (i) k .
The probabilities (r k ) k∈N are then given by

r k = i 0 e -λia 1 -e -λa p (i) k .
Then the probabilities (r k ) k∈N can be computed if the probabilities p

(i) k k∈N
are known for all i ∈ N. In practice, r (i) k k∈N can be identified using Fast Fourier Transform.

The probabilities (q k ) k∈N of the final mutant count can finally be explicited using the following algorithm [START_REF] Embrechts | A limit theorem for tails of discrete infinitely divisible laws with applications to fluctuation theory[END_REF] for Poisson compounds:

q 0 = e -m(1-r 0 ) , (24)
and for any k > 0:

q k = m k k i=1 ir i q k-i . (25)
The simulation algorithms are now introduced. Consider a inhomogeneous Poisson process {Y (t)} t 0 with a given expectation Λ(t). Denote by (S i ) i∈N the occurrence instants of the process {Y (t)} t 0 . Conditionally to

Y (t) = k, Λ(S 1 ) Λ(t) , . . . , Λ(S k ) Λ(t)
has same distribution as the order statistics of a k sample of the uniform distribution on (0 ; 1) (Proposition A.2). Consider first IM models with γ = 0. As a direct consequence of assertion (A (0)

2 ), the following corollary holds: Corollary 4.1. Assume γ = 0. Conditionally on Z n (t) = k and as n tends to infinity, the vector e hν (0,T 1 ) -1 e hν (0,t) -1 , . . . , e hν (0,T k ) -1 e hν (0,t) -1 converges in distribution to the order statistics of a k sample of the uniform distribution on (0 ; 1).. Then, a k sample of mutation instants T 1 , . . . , T k in an interval [0 ; t] is drawn by: 1. sample k uniform variables U 1 , . . . , U k ; 2. deduce the order statistics U (1) , . . . , U (k) ; 3. apply for any i = 1 . . . k:

T i = h -1 ν,0 log U (i) e hν (0,t) -1 + 1 ,
where for any s, h -1 ν,s is the function which satisfies:

h ν s, h -1 ν,s (u) = u and h -1 ν,s (h ν (s, t)) = t .
From (H 4 ), this function is well defined. For example, if h ν is defined as [START_REF] Ma | Analysis of the Luria-Delbrück distribution using discrete convolution powers[END_REF], then:

h -1 ν,s (u) = f -1 [e u f (s)] .
According to Corollary 4.1, a k sample of mutation instants T 1 , . . . , T k in an interval [0 ; t] is drawn by: 1. sample k uniform variables U 1 , . . . , U k ; 2. deduce the order statistics U (1) , . . . , U (k) ; 3. apply for any i = 1 . . . k:

T i = h -1 µ,0 1 ρ(1 -2δ) log U (i) e ρh * µ (0,t) -1 + 1 , (26)
where for any s, h -1 µ,s is the function which satisfies:

h µ s, h -1 µ,s (u) = u and h -1 µ,s (h µ (s, t)) = t .
Hence the following algorithm is used to compute sample of ILD distribution with γ = 0:

1. simulate a random number k of mutations, according to the Poisson distribution with parameter m;

2. compute the clone size stemming from each mutation:

(a) simulate the mutation instant s as above (with t = +∞); (b) compute p 0 (s, +∞) and P (s, +∞); (c) make the random choice:

• with probability p 0 (s, +∞), output 0 (extinction of the clone);

• with probability 1 -p 0 (s, +∞), output a geometric random number with parameter h µ (s, +∞).

3. sum the k clone sizes.

Remark that this algorithm can be used only if h µ,∞ is finite. Otherwise, the mutations instants cannot been computed applying (26) on uniform variables. However, from Corollary 3.1, the LD distribution is recovered when h µ,∞ = +∞. Thus the algorithm exposed by Ycart [START_REF] Ycart | Fluctuation analysis with cell deaths[END_REF] can be used. Remark that IM M models with γ 0 could also be considered, as a consequence of assertion (A (γ)

2 ); Corollary 4.2. Conditionally on Z n (t) = k, the vector ω(T 1 )e h * ν (0,T 1 ) -1 ω(t)e h * ν (0,t) -1 , . . . , ω(T k )e h * ν (0,T k ) -1 ω(t)e h * ν (0,t) -1 converges as n tends to infinity to the order statistics of a k sample of the uniform distribution on (0 ; 1).

Note that the computation of the inverse function of

g(t) = ω(t)e h * ν (0,t) -1 = 1 -2γ (1 -γ)e -h *
ν (0,t) -γe -2h * ν (0,t) -1 , has not been studied here. However, Newton-Raphson algorithm can be used to approximate it. Thus, as for the case where γ = 0, the simulation of a k sample of mutation instants is possible.

Parameter estimation

This section is dedicated to the adaptation of the p 0 , ML, and GF methods mentioned in the introduction. Their construction will be quickly recalled. The ILD and H models with γ = 0 are considered here. Remark that the case γ > 0 could be theoretically considered for the ML and GF methods. However, this case has not been implemented yet and is not considered in Section 6. The methods exposed in this section perform estimation only for m and ρ. Indeed, the fluctuations of the distribution of final mutant counts with respect to δ are very small [START_REF] Ycart | Fluctuation analysis with cell deaths[END_REF]. In practice, only the magnitude of δ can be measured [START_REF] Stewart | Aging and death in an organism that reproduces by morphologically symmetric division[END_REF][START_REF] Fontaine | Mutations in two global regulators lower individual mortality in Escherichia Coli[END_REF]. The same is also true for the parameter h µ,∞ : it is related to the death parameter δ as follows:

h µ,∞ = lim t→+∞ 1 1 -2δ log [E [M (0, t)]] ,
where M (s, t) is the size at time t of a mutant clone started at time s. In other words, estimating h µ,∞ is possible only if both estimates of the death parameter and the expected final size of a mutant clone started at time 0 are available. Thus, the identifiability of the model is hard in practice. In that sense, δ and h µ,∞ are assumed to be known, which is not realistic. In this section, (X 1 , . . . , X n ) will denote a sample of n i.i.d. random variables following ILD or H models. The p 0 method was the first method introduced by Luria and Delbrück [START_REF] Luria | Mutations of bacteria from virus sensitivity to virus resistance[END_REF] to estimate m for classic mutation models when δ = 0. It uses the fact that the probability of null counts is e -m . This relation remains true for ILD models, whether h µ,∞ is finite or not. Hence m can be estimated by m0 = -log (q 0 ) , where q0 is the relative frequency of zero among mutant counts. Thus, m0 is a consistent and asymptotically estimator of m. From the ∆-method (see for example [37, p. 79]), its asymptotic variance is given by v m0 = 1 -q 0 nq 0 .

An extension of the p 0 estimator to the case where δ > 0 has been described by Ycart [START_REF] Ycart | Fluctuation analysis with cell deaths[END_REF] for LD models (called Fixed Point estimator). It uses the fact that the probability of extinction of a supercriticial process is a fixed point of the pgf of the number of cells at time t, in a clone starting from one cell at time 0 [4, Lemma 1, p.141]. This probability is given by Theorem 1 of [4, Chap. I]:

δ * = δ 1 -δ .
Under ILD models, δ * is a fixed point of the pgf ( 16) and [START_REF] Luria | Mutations of bacteria from virus sensitivity to virus resistance[END_REF]. Therefore, a consistent and asymptotically normal estimator of m is then given by m0 =

-log φn (δ * )

1 -δ * ,
where φn denotes the empirical pgf of the sample:

φn (z) = 1 n n i=1 z X i . ( 27
)
For any z ∈ (0, 1), φn (z) is a consistent and asymptotically normal estimator of φ(z).

From the ∆-method, the asymptotic variance of m0 is given by

v m0 = 1 n(1 -δ * ) 2 φ(δ 2 * ) φ(δ * ) 2 -1 .
Remark that the p 0 estimator of m depends only on δ. In other words, the p 0 method does not depend on the growth model, and can be for any mutation model. However, this method does not diretly yield an estimator of ρ. If an estimate is desired, the Maximum Likelihood can be used for ρ only, setting m = m0 .

As the probabilities q k 's and their derivatives with respect to m and ρ are explicit for ILD and H models, the Maximum Likelihood method seems to be an obvious choice to estimate m and ρ. Moreover, an equivalent of r k for large values of k can be computed rewriting [START_REF] Mazoyer | Time inhomogeneous mutation models with birth-date dependence[END_REF] as:

r k = ρδ 2 1 -e -ρh * µ,∞ 1 0 (1 -v) k-1 (1 -δ * v) k+1 v ρ dv - e -h * µ,∞ 0 (1 -v) k-1 (1 -δ * v) k+1 v ρ dv = ρδ 2 1 -e -ρh * µ,∞ k -ρ-1 k 0 1 -w k k-1 1 -δ * w k k+1 w ρ dw - e -h * µ,∞ 0 (1 -v) k-1 (1 -δ * v) k+1 v ρ dv ,
where the constant δ is given by

δ = 1 -2δ 1 -δ .
The following equivalent is then obtained:

r k ∼ k→+∞ ρδ 2 1 -e -hν,∞ Γ(ρ + 1) k ρ+1 1 -δ 1 -2δ ρ+1 - e -h * µ,∞ 0 (1 -v) k-1 (1 -δ * v) k+1 v ρ dv , ( 28
)
where Γ is the Gamma function. Remark that for h µ,∞ = +∞, Formula (3.4) of [START_REF] Ycart | Fluctuation analysis with cell deaths[END_REF] is recovered. Consider first ILD models. The derivatives of the r k 's with respect to ρ are now given:

∂r 0 ∂ρ = ρ 1 -e -ρh * µ,∞ 1 e -h * µ,∞ δ -δv 1 -δ -δv v ρ-1 log(v)dv + 1 ρ - h * µ,∞ e -ρh * µ,∞ 1 -e -ρh * µ,∞ r 0 ,
and for any k > 0:

∂r k ∂ρ = ρ 1 -e -ρh * µ,∞ 1 e -h * µ,∞ (1 -v) k-1 (1 -δ -δv) k+1 v ρ log(v)dv + 1 ρ - h * µ,∞ e -ρh * µ,∞ 1 -e -ρh * µ,∞ r k .
From ( 28), an asymptotic equivalent for large k is given by

∂r k ∂ρ = ρδ 2 1 -e -ρh * µ,∞ Γ(ρ + 1) k ρ+1 1 -δ 1 -2δ ρ+1 (ρ + 1) -log kδ - e -h * µ,∞ 0 (1 -v) k-1 (1 -δ * v) k+1 log(v)v ρ dv + 1 ρ - h * µ,∞ e -ρh * µ,∞ 1 -e -ρh * µ,∞ r k ,
where is the Digamma function. Remark also that for δ = 0:

∂r k ∂ρ = ρ 1 -e -ρhµ,∞ [B(ρ + 1, k)( (ρ + 1) + (ρ + 1 + k)) - e -hµ,∞ 0 log(w)w ρ (1 -w) k-1 dw + 1 ρ - h µ,∞ e -ρhµ,∞ 1 -e -ρhµ,∞ r k
The gradient of the q k 's can then be deduced from the r k 's and their derivatives:

∂q 0 ∂m = -(1 -r 0 )q 0 and ∂q 0 ∂ρ = - ∂m ∂ρ (1 -r 0 ) -m ∂r 0 ∂ρ q 0 = -m h * µ,∞ e -ρh * µ,∞ 1 -e -ρh * µ,∞ (1 -r 0 ) - ∂r 0 ∂ρ q 0
The derivative of φ with respect to m is given by

∂φ ∂m (z) = -(1 -I(z))φ(z) = - j 0 q j z j - i 0 r i z i j 0 q j z j = - j 0 q j z j - i,j 0 r i q j z i+j
On the other hand:

∂φ ∂m (z) = k 0 ∂q k ∂m z k .
Hence for any k > 0:

∂q k ∂m = -     q k - i,j 0 i+j=k r i q j     = k i=1 r i q k-i -q k . ( 29
)
Similarly, the derivative of φ with respect to ρ is given by

∂φ ∂ρ (z) = - ∂m ∂ρ (1 -I(z)) -m ∂h ∂ρ (z) φ(z) = -m h * µ,∞ e -ρh * µ,∞ 1 -e -ρh * µ,∞ (1 -I(z)) - ∂h ∂ρ (z) φ(z) = -m h * µ,∞ e -ρh * µ,∞ 1 -e -ρh * µ,∞ 1 - i 0 z i r i - i 0 ∂r i ∂ρ z i j 0 q j z j = -m h * µ,∞ e -ρh * µ,∞ 1 -e -ρh * µ,∞ - i 0 z i h * µ,∞ e -ρh * µ,∞ 1 -e -ρh * µ,∞ r i + ∂r i ∂ρ j 0 q j z j
Hence for any k > 0:

∂q k ∂ρ = -m h * µ,∞ e -ρh * µ,∞ 1 -e -ρh * µ,∞ q k - k i=1 q k-i h * µ,∞ e -ρh * µ,∞ 1 -e -ρh * µ,∞ r i + ∂r i ∂ρ . ( 30
)
Define the log-likelihood by

(X 1 , . . . , X n ) = n i=1 log(q X i ) = M j=1 log(q j ) n i=1 1 X i =j , (31)
where M = max j X j is the sample maximum. The couple ( mML , ρML ) maximizing ( 31) is consistent and asymptotically normal [START_REF] Lehmann | Theory of Point Estimation[END_REF]Theo. 5.1,Chap.6]. The asymptotic variances of mML and ρML are given by

v mML = I 2,2 det(I) and v ρML = I 1,1 det(I)
, where I = (I i,j ) i,j∈{1,2} is the following information matrix:

I = M j=0           ∂q j ∂m 1 q j 2 ∂q j ∂m ∂q j ∂ρ 1 q 2 j ∂q j ∂m ∂q j ∂ρ 1 q 2 j ∂q j ∂ρ 1 q j 2      n i=1 1 X i =j     
.

Using the algorithms ( 25), ( 29) and ( 30), the log-likelihood and its derivatives can be calculated iteratively. However, the formulas must be applied to vectors as large as M . Therefore, as for the homogeneous case [START_REF] Hamon | Statistics for the Luria-Delbrück distribution[END_REF], the procedure can be very long and numerically unstable. In practice, this instability problem can be avoided using Winsorization [38, sec. 2.2.]: any value of the sample that pass a certain bound is replaced by the bound itself. All information above the bound is lost. In extreme cases where the sample minimum is greater than the bound, i.e. for large m and/or small ρ, irrelevant results will be returned. However, these difficulties should not be met in practice. Indeed, very large numbers of mutants are not countable and do not appear in data sets.

As mentioned previously, the optimization of the likelihood with respect to δ and h µ,∞ is hard in practice. However, the derivatives of the r k 's and the q k 's with respect to δ and h µ,∞ can also be iteratively computed. Consider first the derivatives with respect to δ:

∂r 0 ∂δ = 2ρh µ,∞ e -ρh * µ,∞ 1 -e -ρh * µ,∞ r 0 - δ 1 -e h * µ,∞ ρe (ρ-1)h * µ,∞ 1 -δ -δe h * µ,∞ + 1 1 -e -ρh * µ,∞ 1 e -h * µ,∞ (1 -v)ρv ρ-1 (1 -δ -δv) 2 dv , and
∂q 0 ∂δ = - ∂m ∂δ (1 -r 0 ) -m ∂r 0 ∂δ q 0 = -m 2ρh µ,∞ e -ρh * µ,∞ 1 -e -ρh * µ,∞ (1 -r 0 ) - ∂r 0 ∂δ q 0 .
For all k > 0:

∂r k ∂δ = 1 1 -e -ρh * µ,∞ r k -2δ(1 -δ) -(1 -2δ)(k -1) (1 -2δ)(1 -δ) + 2ρh * µ,∞ e -ρh * µ,∞ + (1 -2δ)(1 -δ) k-1 - 1 -e -h * µ,∞ k-1 1 -δ -δe -h * µ,∞ 2h * µ,∞ ρe -ρh * µ,∞ + 1 e -h * µ,∞ (1 -v) k-1 (1 + v)(k + 1) (1 -δ -δv) k+1 ρv ρ dv .
Since the derivative of φ with respect to δ is given by

∂φ ∂δ (z) = - ∂m ∂δ (1 -I(z)) -m ∂h ∂δ (z) φ(z) = -m 2ρh µ,∞ e -ρh * µ,∞ 1 -e -ρh * µ,∞ (1 -I(z)) - ∂h ∂δ (z) φ(z) = -m 2ρh µ,∞ e -ρh * µ,∞ 1 -e -ρh * µ,∞ - i 0 z i ∂r i ∂δ + 2ρh µ,∞ e -ρh * µ,∞ 1 -e ρh * µ,∞ r i j 0 q j z j .
Hence, for any k > 0:

∂q k ∂δ = -m 2ρh µ,∞ e -ρh * µ,∞ 1 -e -ρh * µ,∞ q k - k i=1 q k-i ∂r i ∂δ - 2h µ,∞ e -ρh * µ,∞ 1 -e -ρh * µ,∞ r i .
The derivatives of the r k 's and the q k 's with respect to h * µ,∞ can be computed as follows:

∂r 0 ∂h * µ,∞ = ρe -ρh * µ,∞ 1 -e -ρh * µ,∞ -r 0 + δ 1 -e -h * µ,∞ 1 -δ -δe -h * µ,∞
, thus:

∂q 0 ∂h * µ,∞ = - ∂m ∂h * µ,∞ (1 -r 0 ) -m ∂r 0 ∂h * µ,∞ q 0 =m ρe -ρh * µ,∞ 1 -e -ρh * µ,∞ (1 -r 0 ) + ∂r 0 ∂h * µ,∞ q 0 .
And for any k > 0:

∂r k ∂h * µ,∞ = ρe -ρh * µ,∞ 1 -e -ρh * µ,∞ -r k + e -h * µ,∞ 1 -e -h * µ,∞ k-1 1 -δ -δe -h * µ,∞ k+1
Since the derivative of φ with respect to h * µ,∞ is given by

∂φ ∂h * µ,∞ (z) = - ∂m ∂h * µ,∞ (1 -I(z)) -m ∂I(z) ∂h * µ,∞ φ(z) =m ρe -ρh * µ,∞ 1 -e -ρh * µ,∞ (1 -I(z)) + ∂I(z) ∂h * µ,∞ φ(z) =m ρe -ρh * µ,∞ 1 -e -ρh * µ,∞ + i 0 z i ∂r i ∂h * µ,∞ - ρe -ρh * µ,∞ 1 -e -ρh * µ,∞ r i j 0 q j z j .
Hence, for any k > 0:

∂q k ∂h * µ,∞ = m ρe -ρh * µ,∞ 1 -e -ρh * µ,∞ q k + k i=1 q k-i ∂r i ∂h * µ,∞ - ρe -ρh * µ,∞ 1 -e -ρh * µ,∞ r i .
Consider now H models. In that case, the derivatives of the r k 's with respect to ρ and δ are given by . For any z ∈ (0 ; 1):

∂r k ∂ρ = i 0 r (i) k -iae -naρ + (i + 1)ae -(i+1)aρ = i 0 r (n) k ae -iaρ (i + 1)e -aρ -n , and
∂r k ∂δ = i 0 ∂p (i) k ∂δ e -iaρ -e -(i+1)aρ + p (i) k iρ (1 -δ) e -iaρ - (i + 1)ρ (1 -δ) e -(i+1)aρ = 1 -e -aρ i 0 ∂p (i) k ∂δ e -iaρ + 2ρe -a i 0 p (i) k e -iaρ i -(i + 1)e -aρ .

The derivatives ∂p

∂b i ∂δ (z) = k 0 ∂p (i) k ∂δ z k = 1 -(b i-1 (z)) 2 + 2(1 -δ) ∂b i-1 ∂δ (z)b i-1 (z) , and
∂b 0 ∂δ (z) = 0 .
The sequence ∂p

(i) n ∂δ k 0
can then be deduced for any integer i from the sequence of polynoms ∂b i ∂δ i 0 using Fast Fourier Transform. However, each polynom b i has 2 i coefficients. The computation of these coefficients is very long. The last method exposed here relies on the generating function φ of final mutant count. Indeed, the parameter of a Poisson compound can be easily estimated [START_REF] Rémillard | Inference based on the empirical probability generating function for mixtures of Poisson distributions[END_REF][START_REF] Marcheselli | Parameter estimation for the discrete stable family[END_REF]. The GF method proposed in [10, Section 4] uses this approach to estimate m and ρ under LD models when δ = 0. This method can also be used for H models [START_REF] Ycart | Fluctuation analysis: can estimates be trusted?[END_REF] or LD models with δ > 0 [START_REF] Ycart | Fluctuation analysis with cell deaths[END_REF]. It is extended here to the ILD models. The pgf I and its derivative with respect to ρ are required in this method. Thus they are implemented as the following numerically stable expressions :

I(z) = δ * + z * (1 -δ * ) 1 -e -ρh * µ,∞ 1 e -h * µ,∞ ρv ρ 1 + z * v dv , and
∂I(z) ∂ρ = z * (1 -δ * ) 1 -e -ρh * µ,∞ 1 - ρh * µ,∞ e -ρh * µ,∞ 1 -e -ρh * µ,∞ 1 e -h * µ,∞ v ρ 1 + z * v dv + 1 e -h * µ,∞ ρv ρ 1 + z * v log(v)dv ,
where the constant δ * and z * are given by

δ * = δ 1 -δ , and z * = z -δ * 1 -z .
Consider z 1 , z 2 , z 3 in (0 ; 1). The GF estimators of m and ρ are the following:

mGF (z 3 ) = log φn (z 3 ) I ρGF (z 1 ,z 2 ) (z 3 ) -1 and ρGF (z 1 , z 2 ) = g -1 (ŷ n ) ,
where I x is the pgf [START_REF] Luria | Mutations of bacteria from virus sensitivity to virus resistance[END_REF] for ILD models or [START_REF] Jones | Luria-Delbrück Fluctuation Experiments: Design and Analysis[END_REF] for H models setting ρ = x, φn is the empricial pgf [START_REF] Rodriguez-Brenes | Cellular replication limits in the Luria-Delbrück mutation model[END_REF] and:

g(x) = I x (z 1 ) -1 I x (z 2 ) -1 and ŷn = log φn (z 1 ) log φn (z 2 )
.

From Theorem 3.4 of [START_REF] Rémillard | Inference based on the empirical probability generating function for mixtures of Poisson distributions[END_REF] and ∆-method, it can be proved that the estimator couple ( mGF , ρGF ) is strongly consistent and asymptotically normal, with explicit asymptotic variance.

Proposition 5.1. Let z 1 , z 2 , z 3 in (0 ; 1), two by two distinct. Consider the random vector

√ n ( φn (z 1 ), φn (z 2 ), φn (z 3 )) -(φ(z 1 ), φ(z 2 ), φ(z 3 )) ,
and its asymptotic covariance matrix C = (c(z i , z j )) i,j=1,2,3 given by

c(z i , z j ) = φ(z i z j ) -φ(z i )φ(z j ) .
Then, the couple of random variables:

√ n (( mGF , ρGF ) -(m, ρ)) (32)
converges in distribution to the bivariate centered normal distribution with covariance matrix A t CA, where A = (a i,j ) i=1,2,3 j=1,2 is a 3 × 2 matrix with:

a 1,1 = ma 1,2 I(z 3 ) -1 ∂I(z 3 ) ∂ρ ; a 1,2 = I(z 2 ) -1 mφ(z 1 ) ∂I(z 1 ) ∂ρ (I(z 2 ) -1) -∂I(z 2 ) ∂ρ (I(z 1 ) -1) a 2,1 = ma 2,2 I(z 3 ) -1 ∂I(z 3 ) ∂ρ ; a 2,2 = I(z 1 ) -1 mφ(z 2 ) ∂I(z 2 ) ∂ρ (I(z 1 ) -1) -∂I(z 1 ) ∂ρ (I(z 2 ) -1) a 3,1 = 1 φ(z 3 )(I(z 3 ) -1) ; a 3,2 = 0
The proof of Proposition 5.1 has already been exposed by Hamon and Ycart [START_REF] Hamon | Statistics for the Luria-Delbrück distribution[END_REF]. By definition, the GF estimators depend on the arbitrary values of z 1 , z 2 and z 3 . They can be seen as tuning parameters and should be appropriately chosen according to the sample to minimize the asymptotic variances of Proposition 5.1. Since the unknown values of m and ρ influence also the variance, optimal values for z 1 , z 2 and z 3 cannot be easily identified. However, the fluctuations of the variances with respect to these three parameters are quite small. Their values have been set in [10, p. 1262] using simulation experiments, such that the tuning parameters minimizes the mean squared error. In that sense, the GF method is a semi-empirical method, since the estimators have been trained by simulated data. Another issue in this method concerns the estimation of ρ which requires to identify the zero of the monotone function g(ρ) -ŷn . In practice, the research domain is bounded. This can be a problem if the sample does not contain jackpot, i.e. for large theoretical values of ρ. However in that case, it should be considered that a mutation model may not be adapted to the data.

The fact is that the GF estimators are in practice comparable in precision to ML estimators, with a negligible computing time [24, Fig. 12]. For that reason, ML optimization can be initialized by GF estimates, to significantly improve computing time. This choice nevertheless has an inconvenient in practice: if the GF method is not able to estimate ρ, there is no chance to estimate it with the ML method. These simulation studies have been implemented in R [START_REF] Development | R: A Language and Environment for Statistical Computing[END_REF], using the R package flan [START_REF] Mazoyer | flan: An R package for inference on mutation models[END_REF], which is available on CRAN (https://cran.r-project.org/package=flan).

Assume that h µ is defined as [START_REF] Ma | Analysis of the Luria-Delbrück distribution using discrete convolution powers[END_REF], where f is solution of the logistic equation:

f (t) = f (t) 1 - f (t) f ∞ ,
where f ∞ denotes the finite limit of f as t tends to infinity, i.e. the carrying capacity. Thus:

f (t) = f ∞ 1 + f∞ f (0) -1 e -t
.

The associated cdf F µ is then given for any (s, t) ∈ R + × R + by

F µ (s, t) = 1 - f (s) f (t) 1 t∈[s ; +∞) + 1 t=+∞ , (33)
and the asymptotic cumulative division rate h µ,∞ by

h µ,∞ = log f ∞ f (0) .
Recall that the choice of the function f is a matter only for the simulation: mutant count distribution depends only on the ratio h µ,∞ (Theorem 3.3). In that sense, f (0) is set here to 1. Simulation experiments have been made along the following lines for each of the 24 sets of parameters m = (0.5, 2, 4, 8), ρ = (0.8, 1, 1.2) and h µ,∞ = (log(100), log(10 4 )):

1. draw 10000 samples of size 100, under inhomogeneous model with cdf (33); 2. for each sample, compute ML and GF estimates of m and ρ under ILD, LD and H models;

3. for each model, observe the empirical distribution of relative bias θ/θ, where θ is the estimator of the true value θ.

Note that the sample size considered here is not realistic: in practice, the sample size rarely exceeds 20. However, the objective of this section is to observe the bias induced by estimating under LD or H models when data are realizations of ILD model. Therefore, large size samples are used in order to reduce the parameter identifiability issues. Similar simulation studies with more realistic sample sizes have also been performed and will be exposed at the end of this section. model(right boxplots), according to the value of h µ,∞ . Empirical distributions of the estimators mGF and ρGF under same assumptions are shown in Figures 2 and4. In the four figures, red lines mark theoretical values, blue lines mark relative biases of 0.9 and 1.1.

According to Figure 1, ML method correctly estimates small values of m under LD model. Positive biases, which seem to increase as ρ decreases, are observed for larger m. The fitness ρ is overestimated by the ML method under LD model. The bias seems to increase as m increases and/or ρ decreases. For the extreme case where m = 8 and ρ = 0.8, a wide proportion of the estimates have a relative bias larger than 20%. However, using ML method under H model seems to provide overestimated values of m. The estimates of ρ seem also to be biased when ρ = 1: the bias is positive when ρ < 1, negative when ρ > 1. For larger value of h µ,∞ (Figure 3), there is almost no distinction between LD and ILD models when using ML method, whatever the theoretical value of m and ρ. However, the behaviour of the ML estimates under H model is harder to interpret. On the other hand, the GF method provides similar results, with slightly larger variances (see Figures 2 and4).

Similar simulations have also been performed using ML method on samples of size 20. Related visual results are exposed in Figures 5 and6. These figures illustrate the identifiability issues: the estimation of ρ is not possible if the sample does not contain jackpot, i.e. for small m and/or large ρ.

Two facts have to be noticed. First, even if the whole population of cells has a logistic growth, it is possible in practice to stop the experiment before the inflexion instant. Thus, the assumption of an exponential growth could be considered and the LD distribution could be used to perform estimations. Moreover, the initial number of cells is in practice of order 10 3 -10 4 , the final number of cells of order 10 8 -10 9 . According to Figures 3 and6, the practical value of h µ,∞ is then such that the bias induced by considering an exponential growth instead of logistical growth should be negligible. Therefore, it seems that the LD model should be preferred since the identifiability issues have less impact for this model.

Conclusion

An extension for the classic mutation models to the case where the final instant of a cell depends on its birth date has been proposed. Results are based on the decomposition as three ingredients of any mutation model. This approach led to a family of distributions for the asymptotic mutant count. These distributions depend on the expected number of mutations m, the death probabilities γ and δ for normal and mutant cells, and the final instant cdf F ν (s, •) and F µ (s, •) for normal and mutant cells born at a given time s. The previous results obtained with an analytic approach are recovered and generalized to the case where γ > 0. Computation of probabilities and simulation algorithms have been described recovered. The latter has also been extended to the case where δ > 0. The consequence for statistical inference has been studied: robust estimation methods have been extended to the ILD models, in the specific cases where the parameters δ and h µ,∞ are known, and γ = 0; biases induced by considering classic mutation models instead of birth-date dependence model have been studied with simulation experiments. Considering the order of h µ,∞ in practice, the biases on m and ρ induced in practice by considering the LD distribution instead of the ILD distribution seems negligible. Thus, still considering the LD distribution for estimating parameters seems worthy. The R package flan which has been used for simulation study is available on CRAN (https://cran.r-project. org/package=flan). Note that flan does not take into account the death of normal cells (i.e. γ > 0) for now. The extension to the case γ > 0 is being developed and should appear on CRAN. In the same time, a web-tool based on flan and R package shiny is also being developed. It allows to use some features of flan (in particular hypothesis testing and simulation) without any installation or knowledge in R. A first version is already available at https://toltex-shiny.u-ga.fr/RodaShiny/ShinyFlan/.

A Point processes

The main properties of point processes used in this paper are shortly exposed in this appendix. Consider a sequence (T i ) i∈N . Each T i represents the occurring time of a given event, such as a mutation. Let us denote by {N (t)} t 0 the associated point process, defined for any t ∈ R + by N (t) = max {i ∈ N ; T i t} , which represents the cumulated number of events occurring in [0 ; t]. Assume that the following holds:

1. T 0 = 0 and N (0) = 0 with probability 1; 2. The paths {N (t)} t 0 are cadlag and increasing ; 3. the counting process {N (t)} t 0 is simple, i.e. for any t ∈ R + :

This intensity represents the probability of occurring after a time t, conditionally to the number of events before t and their respective appearance time. Thus ξ is a random variable depending on N (t), T 1 , . . . , T N (t) . To simplify the redaction, it will be sometimes denoted here by ξ(t, N (t)).

2. The conditional intensity defined by

ξ(t, n) = E[ξ(t, N (t)) | N (t) = n] (35)
The function ξ(t, n) is deterministic. It represents the probability of occurring after a time t, knowing that n events have occurred before t.

The distribution of the event number occurring before a given time can then be explicited.

Proposition A.1. For any t ∈ R + , the distribution of the event number occurring before a given time t is given by

P [N (t) = 0] = exp - t 0 ξ(u, 0)du , (36)
and for any n > 0,

P [N (t) = n] = 0<t 1 <•••<tn n i=1 ξ(t i , i -1) exp - n i=0 t i+1 t i ξ(u, i)du dt 1 . . . dt n , (37
) with convention t 0 = 0 et t k+1 = t.

Consider now the distribution of the occurring instants. For any vector t (n) = (t 1 , . . . , t n ) ∈ R n + , the event (T 1 = t 1 , . . . , T n = t n ) will be denoted by T (n) = t (n) . Proposition A.2. The probability distribution function (pdf ) of the instant T 1 of the first event is given by f T 1 (t) = ξ(t, 0) exp -t 0 ξ(u, 0)du .

(

) 38
Conditionally to T (n) = t (n) , the pdf of the instant T n+1 of the n + 1-th event is given by

f (Tn+1|T (n) =t (n) ) (t) = ξ(t, n, t (n) ) exp - t tn ξ(u, n, t (n) )du 1 0<t 1 <•••<tn . ( 39
)
The joint pdf of the vector T (n) is given by

f T (n) (t (n) ) = n i=1 ξ(t i , i -1, t (i-1) ) exp - n i=1 t i t i-1 ξ(u, i -1, t (i-1) )du 1 0<t 1 <•••<tn . (40)

5 ρ

 5 Boxplots of Figures 1 and 3 show the empirical distributions of the estimators mML and ρML obtained under LD model (left boxplots), H model (center boxplots), and ILD m = 0.

Figure 1 :

 1 Figure 1: ML estimates under LD, H, and ILD models on data drawn with ILD model (h µ,∞ = log (100)). For each of the 12 sets of parameters m = (0.5, 2, 4, 8) (rows) and ρ = (0.8, 1, 1.2) (columns), 10 4 samples of size 100 of the ILD(m, ρ, 0, 0, h µ,∞ ) distribution were simulated with cdf (33) and h µ,∞ = log(100). For each column, the first three boxplots represent the distribution of the 10 4 ratio mML /m obtained under LD model (left), H model (center), and ILD model (right) ; the last three boxplots represent the distribution of the 10 4 ratio ρML /ρ obtained under LD model (left), H model (center), and ILD model (right). Red horizontal lines mark theoretical value. Blue horizontal lines mark relative biases of 0.9 and 1.1.

Figure 2 :

 2 Figure 2: GF estimates under LD, H, and ILD models on data drawn with ILD model (h µ,∞ = log (100)). For each of the 12 sets of parameters m = (0.5, 2, 4, 8) (rows) and ρ = (0.8, 1, 1.2) (columns), 10 4 samples of size 100 of the ILD(m, ρ, 0, 0, h µ,∞ ) distribution were simulated with cdf (33) and h µ,∞ = log(100). For each column, the first three boxplots represent the distribution of the 10 4 ratio mGF /m obtained under LD model (left), H model (center), and ILD model (right) ; the last three boxplots represent the distribution of the 10 4 ratio ρGF /ρ obtained under LD model (left), H model (center), and ILD model (right). Red horizontal lines mark theoretical value. Blue horizontal lines mark relative biases of 0.9 and 1.1.

Figure 3 :

 3 Figure 3: ML estimates under LD, H, and ILD models on data drawn with ILD model (h µ,∞ = log (10 4 )). For each of the 12 sets of parameters m = (0.5, 2, 4, 8) (rows) and ρ = (0.8, 1, 1.2) (columns), 10 4 samples of size 100 of the ILD(m, ρ, 0, 0, h µ,∞ ) distribution were simulated with cdf (33) and h µ,∞ = log(10 4 ). For each column, the first three boxplots represent the distribution of the 10 4 ratio mML /m obtained under LD model (left), H model (center), and ILD model (right) ; the last three boxplots represent the distribution of the 10 4 ratio ρML /ρ obtained under LD model (left), H model (center), and ILD model (right). Red horizontal lines mark theoretical values. Blue horizontal lines mark relative biases of 0.9 and 1.1.

Figure 4 :

 4 Figure 4: GF estimates under LD, H, and ILD models on data drawn with ILD model (h µ,∞ = log (10 4 )). For each of the 12 sets of parameters m = (0.5, 2, 4, 8) (rows) and ρ = (0.8, 1, 1.2) (columns), 10 4 samples of size 100 of the ILD(m, ρ, 0, 0, h µ,∞ ) distribution were simulated with cdf (33) and h µ,∞ = log(10 4 ). For each column, the first three boxplots represent the distribution of the 10 4 ratio mGF /m obtained under LD model (left), H model (center), and ILD model (right); the last three boxplots represent the distribution of the 10 4 ratio ρGF /ρ obtained under LD model (left), H model (center), and ILD model (right). Red horizontal lines mark theoretical values. Blue horizontal lines mark relative biases of 0.9 and 1.1.

5 ρ

 5 . The Luria-Delbrück distribution with cell deaths and the Haldane model are m = 0.

Figure 5 :

 5 Figure 5: ML estimates under LD, H, and ILD models on data drawn with ILD model (h µ,∞ = log (100)). For each of the 12 sets of parameters m = (0.5, 2, 4, 8) (rows) and ρ = (0.8, 1, 1.2) (columns), 10 4 samples of size 20 of the ILD(m, ρ, 0, 0, h µ,∞ ) distribution were simulated with cdf (33) and h µ,∞ = log(10 4 ). For each column, the first three boxplots represent the distribution of the 10 4 ratio mML /m obtained under LD model (left), H model (center), and ILD model (right) ; the last three boxplots represent the distribution of the 10 4 ratio ρML /ρ obtained under LD model (left), H model (center), and ILD model (right). Red horizontal lines mark theoretical values. Blue horizontal lines mark relative biases of 0.9 and 1.1.

Figure 6 :

 6 Figure 6: ML estimates under LD, H, and ILD models on data drawn with ILD model (h µ,∞ = log (10 4 )). For each of the 12 sets of parameters m = (0.5, 2, 4, 8) (rows) and ρ = (0.8, 1, 1.2) (columns), 10 4 samples of size 20 of the ILD(m, ρ, 0, 0, h µ,∞ ) distribution were simulated with cdf (33) and h µ,∞ = log(10 4 ). For each column, the first three boxplots represent the distribution of the 10 4 ratio mML /m obtained under LD model (left), H model (center), and ILD model (right) ; the last three boxplots represent the distribution of the 10 4 ratio ρML /ρ obtained under LD model (left), H model (center), and ILD model (right). Red horizontal lines mark theoretical values. Blue horizontal lines mark relative biases of 0.9 and 1.1.

1 .

 1 (t + ∆t) -N (t) > 1] = 0 .Consider now the two following functions: The intensity of the process {N (t)} t 0 defined by ξ(t, N (t)) ≡ ξ(t, N (t), T 1 , . . . , T N (t) t + ∆t) -N (t) = 1 | N (t), T 1 , . . . , T N (t) .

Simulation studyIf the model used for the estimation does not correspond to the theoretical model, the estimates can be biased. Simulation experiments have been performed to observe the bias induced by estimating under LD or H models when data are realizations of ILD model. The death parameters γ and δ are assumed to be zero.

Assume now that the process {N (t)} t 0 is a time inhomogeneous Poisson process. In that case, the intensity [START_REF] Tan | A stochastic Gompertz birth-death process[END_REF] is a deterministic time function:

Note that the increments are independent. The conditional intensity [START_REF] Tan | On stochastic growth processes with application to stochastic logistic growth[END_REF] is given by

According to Proposition A.1, the distribution of the event number occurring before a given time can be explicited. The direct consequence of Proposition A.2 provides the distribution of the occurring instant for a given number of events.

Proposition A.3. For any t ∈ R + , N (t) follows the Poisson distribution with expectation

The pdf of the instant T 1 of the first event is given by

Conditionally to T (n) = t (n) , the pgf of the (n + 1)-th instant T n+1 is given by

The joint pdf of the vector T (n) is given by

The inhomogeneous Poisson processes satisfy a wide panel of usefull properties. One of them concerns the joint distribution of T (n) , conditionally to N (t) = n.

Corollary A.1. Conditionally to N (t) = n, the joint pdf of the vector T (n) is the same as the order statistics of a n sample of the distribution ξ(u) m(t)

Corollary A.1 can also be written as Corollary A.2.

Corollary A.2. Conditionally to N (t) = n, the joint pdf of the vector m(T 1 ) m(t) , . . . , m(T n ) m(t) is the same as the order statistics of a n sample of the uniform distribution on (0 ; 1).