Adrien Mazoyer
email: adrien.mazoyer@univ-grenoble-alpes.fr

Time inhomogeneous mutation models with birth-date dependence

Keywords:

Time inhomogeneous mutation models with birth-date dependence

Introduction

Mutation models are probabilistic descriptions of the growth of a population of cells, in which scarse mutations randomly occur. The first objective of these models is to study the distribution of the number of mutant cells at the end of the growth process. The classic mutation models 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.

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]. 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 (see Yule [34, p. 35]; Athreya and Ney [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]; an explicit form of the asymptotic distribution is given by Ycart [START_REF] Ycart | Fluctuation analysis: can estimates be trusted?[END_REF]. General lifetimes have also been studied by Ycart [START_REF] Ycart | Fluctuation analysis: can estimates be trusted?[END_REF], but no explicit distribution is available appart 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 (Angerer [2, Sec. 3.1]; Dewanji et al. [START_REF] Dewanji | A generalized Luria-Delbrück model[END_REF], Komarova et al. [START_REF] Komarova | The fixed-size Luria-Delbrück model with a nonzero death rate[END_REF], and Ycart [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], or where the cell divisions are asymmetric [START_REF] Montgomery-Smith | Age-dependent branching processes and applications to the Luria-Delbrück experiment[END_REF].

In the mutation models cited above, the lifetimes of the cells are supposed to be i.i.d., which is quite unrealistic. Indeed, during an experiment, a colony of cell grows in an environment which contains a finite amount of resources. Consider two instants s 1 and s 2 such that s 1 s 2 , then a cell born at time s 1 will complete its lifetime faster than a cell born at time s 2 . 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 Tan [START_REF] Tan | A stochastic Gompertz birth-death process[END_REF], Tan and Piantadosi [START_REF] Tan | On stochastic growth processes with application to stochastic logistic growth[END_REF], and Lambert [START_REF] Lambert | The branching process with logistic growth[END_REF]. The independence of lifetimes for single type branching processes was questioned quite early [START_REF] Kendall | On the Choice of a Mathematical Model to Represent Normal Bacterial Growth[END_REF]. Experiments have evidenced correlation between a cell and its descendants, and between two sisters conditioning on their mother [START_REF] Wang | Robust growth of Escherichia Coli[END_REF]. The effects of these correlations and many models have been discussed since then: see Louhichi and Ycart [START_REF] Louhichi | Exponential growth of bifurcating processes with ancestral dependence[END_REF] and references cited therein.

When the lifetimes are i.i.d., using the theory of branching processes [START_REF] Bellman | On age-dependent binary branching processes[END_REF][START_REF] Athreya | Branching Processes[END_REF], the distribution of the total number of mutant cells converges as the initial number of cells tends to infinity to a heavy-tailed distribution. This convergence has also been proved using an analytic way by Bartlett [5,Sec. 4.31], but for a restrictive case where the fitness is set to 1. (equal growth rates for normal and mutant cells). 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 formulation without assumption on the growth model for the mutant clones. However, no result for the non i.i.d. lifetimes case has been stated until now. In particular, there is no convergence result for the distribution of final mutant counts.

The main objective of this paper is to extend classic mutation models to the case where the split instant of a cell depends on its birth date. Cells deaths and non exponential lifetimes are also taken into account with this approach. General modeling assumptions are described in Section 2. The main tool used in this paper is an extension of the Bellman-Harris integral equation [START_REF] Bellman | On age-dependent binary branching processes[END_REF], which is discussed in Section 3. General solutions are also provided. They are used in the asymptotic context (large observation time and small mutation probability) of mutation models in Section 4. Some examples among which Haldane model and a more general case (non-exponential growth) are provided. The convergence results are finally applied in Section 5 to the case where the hazard functions associated to the split instant distribution of normal and mutant cells are proportional. In particular, the Luria-Delbrück distribution with cell deaths, denoted here by LDD distribution [START_REF] Ycart | Fluctuation analysis with cell deaths[END_REF], is recovered. A computation algorithm is proposed in Section 6.

Hypotheses and models

In this section, the probabilistic model is defined as a tree-indexed process (see Pemantle [START_REF] Pemantle | Tree-indexed process[END_REF] and Benjamini and Peres [START_REF] Benjamini | Markov chains indexed by trees[END_REF] for general references). Denote by T the infinite complete binary tree and 0 its root. The vertices of T are interpreted as cells, and each vertex has exactly 2 descendants. If x is a vertex of T, the number of edges between the root 0 and x is denoted by |x|. The subtree of T composed by a vertex x and all its descent will be denoted by T x . If x and y are two vertices of T, x y is the order relation that holds if x is in the path from 0 to y; x ≺ y holds if x is strictly in the path from 0 to x; x ∧ y is the most recent common ancestor of x and y. The mother of the cell x is denoted by x: it is the cell such that x x and |x| = |x| -1. Each cell produces two cells upon completion of its lifetime. Then a cell x 0 which is not the root 0 has a sister cell x 1 . A vertex x 0 and its sister x 1 satisfy

x 0 = x 1 = x 0 ∧ x 1 .
The evolution of a clone stemming from a single cell at instant 0 will be modeled by a stochastic process (C x ) x∈T indexed by the binary tree T. For any x ∈ T, C x is a couple (B x , T x ) where B x describes the nature (dead, normal or mutant) of the cell x:

• B x = 0 if the cell x is dead; • B x = 1 if the cell x is normal; • B x = 2 if the cell x is mutant;
and T x the instant at which the cell x completes its lifetime. The instant T x will be called final instant of x. Denote by R + = R + ∪ {+∞} the extended real line, and by B(R + ) its Borel σ-field. From the above settings, the stochastic process (C x ) x∈T is defined on the measurable space (Ω, A), where Ω = {0, 1, 2} × R + and its σ-algebra A = P({0, 1, 2}) × B(R + ). The fact that T x can be infinite will be discussed below. There remains to define a probability distribution on that space. Recall that the birth date of the root is set to 0. Assume that its nature B 0 is known. Note that a dead cell has no descent. In other words, if B x = 0 then B y = 0 for any y ∈ T x . Let π, γ, δ be real numbers in (0 ; 1), respectively interpreted as probability of mutation, of dying for a normal cell, of dying for a mutant cell.

Consider a cell x 0 = 0 and its sister x 1 . Their nature B x 0 and B x 1 depend only on the nature B x 0 of the mother cell:

• if B x 0 = 0, then B x 0 = 0 and B x 1 = 0 with probability 1;

• if B x 0 = 1, then:
-B x 0 = 1 and B x 1 = 2 with probability π/2; -B x 0 = 2 and B x 1 = 1 with probability π/2; -B x 0 = 0 and B x 1 = 0 with probability γ;

-B x 0 = B x 1 = 1 with probability 1 -π -γ; • if B x 0 = 2, then:
-B x 0 = 0 and B x 1 = 0 with probability δ;

-B x 0 = B x 1 = 2 with probability 1 -δ.
In other words, upon completion of its lifetime, any normal cell produces one normal and one mutant cell with probability π (this event is called a mutation), dies with probability γ or produces two normal cells with probability 1 -π -γ. Upon completion of its lifetime, any mutant cell dies with probability δ or produces two mutant cells with probability 1-δ. Moreover, the events of death or mutation do not depend on the final instant of the cell.

For any cell x (such that B x = 0), its final instant T x depends on its nature B x and its birth date, i.e. on the split instant of its mother T x : if B x = 1 and T x = s, the cumulative distribution function (cdf) of T x is denoted by F ν (s, •); if B x = 2 and T x = s, the cdf of T x is denoted by F µ (s, •). These cdfs satisfy 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, for any cell x, the distribution of T x may have a positive mass at infinity. This hypothesis requires the notion of subprobability measure on a measurable space (Ω , A), i.e. a measure η on (Ω , A) such that η(Ω) 1. For more details about subprobability measures, see for example Nguyen [23, p. 170]. Consider a subprobability measure η on (R + , B(R + )). Then the limit of its cdf F (t) as t tends to infinity is smaller than 1. Let us define for any a, b ∈ R + the following measure η:

η(]a ; b[) = η(]a ; b[) ; η([0 ; b[) = η([0 ; b[) ; and η(]a ; + ∞]) = η(]a ; + ∞[) . Since the set {]a ; b[ ; [0 ; b[ ; ]a ; + ∞] | a, b ∈ R + } is a topological basis of R + ,
the Carathéodory extension's theorem can be applied to extend the measure η on B(R + ). From the measure η, the following probability measure η can be defined for any A ∈ B(R + ):

η(A) = η(A) + 1 -η(R + ) 1 A∈B(R + ) B(R + )
with the associated cdf:

F (t) = F (t)1 t∈[0 ; +∞) + 1 t=+∞ . (2.1)
Remark that if η is a probability measure on R + , then η is also a probability measure on R + . The stochastic process (C x ) x∈T = (B x , T x ) x∈T models the evolution of a clone stemming from a single cell at instant 0. Actually, the evolution of a clone stemming from a single cell y can be deduced from the above process. It consists of a stochastic process (C x ) x∈Ty indexed by the binary tree T y with the same modeling assumptions as above, conditionally to C y .

There remains to define the dependences between the cells. Consider a cell x 0 = 0 and its sister x 1 . The final instants T x 0 and T x 1 are assumed to be independent conditionally to C x 0 . By extension, the clones (C y ) y∈Tx 0 and (C y ) y∈Tx 1 are independent conditionally to C x 0 . Consider now two cells x = 0 and y = 0 and their common ancestor x ∧ y. Assume that their common ancestor is neither x nor y. Then only one of the daughter cells of x ∧ y is in the path from 0 to x, and its sister is in the path from 0 to y. Thus, according to the previous dependence assumption, the final instants T x and T y are independent conditionally to C x∧y . Therefore, the clones (C w ) w∈Tx and (C w ) w∈Ty are independent conditionally to C x∧y .

Thereafter, the root 0 is assumed to be a normal cell, i.e. B 0 = 1. The model can be summarized as follows:

• at time 0, a single normal cell is 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, •) defined on R + ;

• 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.

All the results will be given in terms of the bivariate probability generating functions (pgf) of the numbers of normal and mutant cells, and the pgf of the number of mutant cells. The next section is dedicated to the calculation of these functions.

Integral equations for probability generating functions

Denote by N (s, t) the couple of numbers of normal and mutant cells at time t in the clone stemming from a single normal cell born at time s. Its bivariate pgf is defined by

ϕ(y, z, s, t) = n,m 0 y n z m P [N (s, t) = (n, m)] . (3.1)
Note that ϕ(1, z, s, t) is the pgf of the number of mutant cells in the clone stemming from a normal cell born at time s. Denote by M (s, t) the number at time t of mutant cells in the clone stemming from a single mutant cell born at time s. Its pgf is defined by

ψ(z, s, t) = m 0 z m P [M (s, t) = m] . (3.2)
One way to study these pgfs is to apply the well known Bellman-Harris integral equation [START_REF] Bellman | On age-dependent binary branching processes[END_REF]. However, this equation has been justified only for the case where the lifetimes are i.i.d.. This section extends it to the case where the final instant of a cell depends on its birth date. These equations have already been stated for the case of i. • the final instant of a cell of type i > 0 born at time s is a random variable with cdf F i (s, •) defined on R + such that F i (s, t) = 0 if t s;

• upon completion of a cell of type i > 0:

for any j ∈ {1; . . . ; l}, a random number K i,j of cell of type j is produced;

with probability γ i , the cell dies out;

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

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

For any 1 i l, denote by χ i the pgf of (K i,j ) j=1,...,l defined by

χ i (Z) = γ i + k 1 ,...,k l 0 P [K i,1 = k 1 , . . . , K i,l = k l ] l j=1 z k j j
, for any Z = (z j ) j=1,...,l ∈ [0 ; 1] l . Denote by X i,j (s, t) the number at time t of cell of type j in the clone stemming from a cell of type i born at time s. Denote by ϕ i (Z, s, t) the pgf of (X i,j (s, t)) j=1,...,l defined by

ϕ i (Z, s, t) = k 1 ,...,k l 0 P [X i,1 (s, t) = k 1 , . . . , X i,l (s, t) = k l ] l j=1 z k j j
.

Assume the initial cell completes its lifetime at a given time u > s. For times t < u, the mother cell is still alone in the corresponding clone. For times t u, the number of cells of type j is equal to the sum of cells of type j in the clones stemming from the cells produced by this first division. Then, the number X i,j (s, t|u) of cell of type j in the considered clone, knowing that the initial cell (of type i) completes its lifetime at a time u, is given by

X i,j (s, t|u) =   l k=1 K i,k h=1 X (h) k,j (u, t)   1 t u + 1 t<u 1 i=j , where X (h)
i,j (u, t) are i.i.d. copies of the variable X i,j (u, t). Denote by ϕ i (Z, s, t|u) the pgf of (X i,j (s, t|u)) j=1,...,l defined by

ϕ i (Z, s, t|u) = k 1 ,...,k l 0 P [X i,1 (s, t|u) = k 1 , . . . , X i,l (s, t|u) = k l ] l j=1 z k j j .
Since the growths of the clones are mutually independent, the following equation is obtained:

ϕ i (Z, s, t|u) = χ i [ϕ 1 (Z, u, t), . . . , ϕ l (Z, u, t)] 1 t u + z i 1 t<u .
Integrating with respect to the distribution F i (s, •) removes the conditioning on the final instant of the cell. Then the following integral equation is obtained:

ϕ i (Z, s, t) = t s χ i [ϕ 1 (Z, u, t), . . . , ϕ l (Z, u, t)] dF i (s, u) + z i (1 -F i (s, t)) . (BHM)
Similarly as for the homogeneous case [START_REF] Bellman | On age-dependent binary branching processes[END_REF], an intuitive interpretation of (BHM) can be given. For a given time t > s, there are two possibilities: either the division of the cell has taken place after t with probability 1 -F i (s, t). Then there is still one cell of type i, and φ i (Z, s, t) = z i (second term in (BHM)); or the division of the cell has taken place in [u ; u + du] (where s < u t) with probability dF i (s, u). In that case, each cell of type j will start a clone, accounted for φ j (Z, u, t). Since the number of cells of each type is given by pgf χ j , this accounts for the integral term of (BHM).

Consider now the model described in the previous section. In that case, l = 2. Assume that the type 1 correspond to the normal type. Then, for any (y, z) ∈ [0 ; 1] 2 :

χ 1 (y, z) = γ + πyz + (1 -π -γ)z 2 , and χ 2 (y, z) = δ + (1 -δ)z 2 .
Hence, applying (BHM) to pgf (3.1) and (3.2) leads to the following integral equations:

ϕ(y, z, s, t) = t s γ + πϕ(y, z, u, t)ψ(z, u, t) + (1 -π -γ)ϕ(y, z, u, t) 2 dF ν (s, u) + y(1 -F ν (s, t)) , (3.3)
and

ψ(z, s, t) = t s δ + (1 -δ)ψ(z, u, t) 2 dF µ (s, u) + z(1 -F µ (s, t)) . (3.4)
Until now, there were no specific assumptions on F ν and F µ , except their definition domain and the fact that F ν (s, t) = F µ (s, t) = 0 if t s. Thereafter, in order to solve (3.3), some hypotheses on F ν are 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, •) as (2.1); (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) .

Remark that h is by definition 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) .

For the moment, no assumptions on F µ are required. However, some of the particular cases introduced here will assume that F µ satisfies (H) too. In that case, its related function defined in (H 4 ) will be denoted by h µ . From now on, assume that F ν satisfies (H), and denote by h ν its related function defined in (H 4 ). Thus 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 rate on an interval [s ; t] associated to F ν . The function λ ν can be interpreted as the instantaneous rate associated to F ν on R + . 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 = +∞ .
Since h ν (s, t) = h ν (0, t) -h ν (0, s), replacing F ν by its expression in (3.3) leads to:

ϕ(y, z, s, t)e -hν (0,s) = t s [γ + πϕ(y, z, u, t)ψ(z, u, t)

+(1 -π -γ)ϕ(y, z, u, t) 2 λ ν (u)e -hν (0,u) du + ye -hν (0,t) .
Taking the derivative with respect to s and dividing by e -hν (0,s) leads to the following Riccati equation:

∂ s ϕ(y, z, s, t) = -λ ν (s) [γ -(1 -πψ(z, s, t))ϕ(y, z, s, t) +(1 -π -γ)ϕ(y, z, s, t) 2 , ( R1
)
with the condition ϕ(y, z, t, t) = y. Riccati equations may have explicit solutions depending on the coefficients (see for example Kucera [START_REF] Kucera | A review of the matrix Riccati Equation[END_REF] and Harko et al. [START_REF] Harko | Analytical Solutions of the Riccati Equation with Coefficients Satisfying Integral or Differential Conditions with Arbitrary Functions[END_REF]). One of the cases where (R1) can be solved is when γ = 0. This case makes sense in the context of mutation models: the objective is to obtain an explicit pgf for the mutant counts only.

The number of divisions occurring in dying clones remains bounded and can be neglected. Thus, it can be considered that observed mutants only come from divisions in surviving clones.

In order to simplify expressions, define the following function:

I(z, s, t) = t s λ ν (u)ψ(z, u, t) du . (3.5)
If γ = 0, (R1) reduces to a Bernoulli equation of order 2. Then the change of variable ϕ = 1/ϕ leads to the solution ϕ(y, z, s, t).

Proposition 3.1. Assume γ = 0. The general solution of the Riccati equation (R1) is given by ϕ(y, z, s, t) = e πI(z,s,t)-hν (s,t) 1 y

-(1 -π) t s λ ν (u)e πI(z,u,t)-hν (u,t) du -1 . (3.6)
Another case where (R1) has an explicit solution is when γ = δ and F µ satisfies (H) such that h µ (s, t) = h ν (s, t) for any (s, t) in R 2 + . In that case, ψ is a particular solution, and the general solution of (R1) is explicit [START_REF] Harko | Analytical Solutions of the Riccati Equation with Coefficients Satisfying Integral or Differential Conditions with Arbitrary Functions[END_REF]. However, it corresponds to the case where mutant and normal cells have the same death probability and the same final instant distribution. Mutant and normal cells are then indistinguishable, which seems to be of less practical relevance. An explicit solution can also be obtained if 1 -π -γ = 0. However, γ has to be less than 0.5 (supercritical process), and π is typically in practice of order 10 -11 -10 -9 . Thus this case will not be studied here.

As a direct consequence of Proposition 3.1, setting y = 1 and s = 0 in (3.6) leads to the following result.

Corollary 3.1. Assume γ = 0. The mutant counts at time t starting with a single normal cell at time 0 follows the distribution with pgf: φ(z, t) = e πI(z,0,t)-hν (0,t) 1 -(1 -π) t 0 λ ν (u)e πI(z,u,t)-hν (u,t) du -1 .

(3.7)

From now on, γ will be set at 0. Observe that no assumptions on the cdf F µ are required. In particular, F µ does not necessarily satisfy (H). As long as the pgf ψ is known, the distribution of the mutant counts at a given time is explicit. As an example, consider the Haldane model: 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. In this case, λ ν ≡ λ and h ν is given by

h ν (s, t) = λ(t -s) ,
and the cdf F µ (s, •) is defined for any t s by

F µ (s, t) = 1 if t s + a , 0 else .
Then F µ does not satisfy (H): property (H 1 ) is not satisfied. However, the pgf ψ is easily identifiable. Assume that the lifetime of any mutant cells is equal to a. Considering a cell born at time s, let b i (z) be the pgf of the size of its clone in the interval [s+ia ; s+(i+1)a). Then b 0 (z) = z, and for all positive integer i,

b i (z) = δ + (1 -δ) (b i-1 (z)) 2 . (3.8)
Therefore the pgf of the size at time t of a clone starting at time s is:

ψ(z, s, t) = i 0 b i (z)1 t∈[s+ia ; s+(i+1)a) .
Functions I(z, s, t) and ϕ(z, t) can be explicited. This example will be continued in the next section where the asymptotic model is considered. Assume now that F µ satisfies (H). There exists a function λ µ , with the same properties as λ ν , such that F µ (s, •) is given by

F µ (s, t) =    1 -exp - t s λ µ (u) du 1 s t if t < +∞ , 1 if t = +∞ .
By the same reasoning as for ϕ, the following Riccati equation for ψ is obtained from (3.4):

∂ s ψ(z, s, t) = -λ µ (u) δ -ψ(z, s, t) + (1 -δ)ψ(z, s, t) 2 , ( R2
)
with the condition ψ(z, t, t) = z. The general solution of (R2) can be explicited without specific hypotheses:

the growth rate of normal cells to hat of mutants. Thereafter, the parameter ρ will be mentionned as the fitness. In that case, (3.5) can be explicited:

I(z, s, t) = ρ(1 -2δ) 1 -δ log (1 -2δ)e (1-δ)hµ(s,t) (1 -δ) (1 -z)e h * µ (s,t) + z -δ = h ν (s, t) + ρ(1 -2δ) 1 -δ log (1 -2δ) (1 -δ) (1 -z)e h * µ (s,t) + z -δ . (3.10)
In particular, let f be a continuous, non negative, and increasing function on R + . Let µ be defined for (s, t) in R 2 + by

h µ (s, t) = log f (t) f (s) .
Then λ µ is given by

λ µ (s) = f (s) f (s) ,
and the expectation of the size at a finite time t 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. Moreover, if δ = 0, plugging (3.10) in (3.7) and applying the change of variable w = f (u)/f (t) leads to:

φ(z, t) = f (t) f (0) -ρ 1 -z + z f (0) f (t) πρ × 1 -(1 -π) 1 f (0) f (t) ρw ρ-1 (1 -z + zw) -πρ dw -1 . (3.11)
For example, if f is defined for any t 0 by f (t) = e t , (3.11) is given by

φ(z, t) =e -ρt 1 -z + ze -t -πρ × 1 -(1 -π) t 0 ρe -ρu 1 -z + ze -u -πρ du -1
, which is the inverse of formula [START_REF] Embrechts | A limit theorem for tails of discrete infinitely divisible laws with applications to fluctuation theory[END_REF] of Bartlett [5,p. 155]. Functions with a carrying capacity, such that logistic or Gompertz functions, can also be considered and plugged into (3.11). Corollary 3.1 is used in the next section in a relevant asymptotic context to get the convergence in distribution of the mutant count when n normal cells are initially present.

Asymptotic for mutation models

In this section, the previous results are applied to mutation models. A convergence theorem for the final number of mutant cells is proved, generalizing the analytic method initiated by Bartlett [5,Sec. 4.31]. A mutation model consists of n independent copies of the model described in Section 2. Denote by φ n (z, t) the pgf for the mutant counts at time t starting with n normal cells at time 0. Because of the independence of the n initial cells, φ n (z, t) is the n-th power of (3.7). 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 that 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. Considering this asymptotic context, the main objective of this section is to establish the convergence as n tends to infinity of φ n (z, τ n ). Before stating the result, recall that γ = 0, and that F ν satisfies (H). Denote by h ν its related function defined by (H 4 ), and by λ ν the associated instantaneous rate. The limit of h ν (0, t) as t tends to infinity will be denoted by h ν,∞ . Note that the result exposed in this section does not require that F µ satisfies (H). The Haldane model described earlier will be considered as an example. Define now the function I as I(z, t) = 1 1 -e -hν (0,t) t 0 ψ(z, u, t)λ ν (u)e -hν (u,t) du .

(4.1)

Remark that, assuming the probability distribution function of a mutation instant is given by λ ν e -hν (u,t) 1 [0 ; t] , the function I could be interpreted as the pgf of the size at a given time t of any mutant clone. The main result of this paper is the following convergence theorem:

Theorem 4.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) = α .
Assume that the limit lim t→+∞ I(z, 0, t)e -hν (0,t) (

exists and is finite. As n tends to infinity, the pgf of the number of mutants at time τ n , starting with n normal cells at time 0, converges to

φ(z) = exp {-m(1 -I ∞ (z))} , (4.3)
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 , and m = α 1 -e -hν,∞ .
Observe that (4.3) is the pgf of a Poisson compound with parameter m. By analogy with the homogeneous case, this parameter could be interpreted as the mean number of mutations, assuming that the number of mutation occasions is almost surely equivalent to n e hν (0,τn) -1 as n tends to infinity. The main tool required to prove Theorem 4.1 is Lemma 4.1 below. Lemma 4.1. For any π ∈ [0 ; 1[, z ∈ [0 ; 1], t ∈ R + , and s ∈ [0 ; t], the following bound holds:

e ±πI(z,s,t) -(1 ± I(z, s, t)) π 2 e I(z,0,t) .

The proof uses a power series expansion of e ±πI(z,s,

Hence:

e ±πI(z,s,t) -(1 ± πI(z, s, t))

k 2 | ± π| k I(z, s, t) k k!
π 2 e I(z,s,t) π 2 e I(z,0,t) .

An analytic proof for the case where mutant and normal cells are exponentially i.i.d. with equal rates has been provided by Bartlett [5,Sec. 4.31]. This approach has been adapted to prove Theorem 4.1.

Theorem 4.1. Define the following two functions:

f 1 (z, u, t, π) = e πI(z,u,t) -(1 + πI(z, u, t)) and f 2 (z, t, π) = f 1 (z, 0, t, π) .
According to Lemma 4.1: |f 1 (z, u, t, π)| π 2 e I(z,0,t) and |f 2 (z, t, π)| π 2 e I(z,0,t) .

(4.4)

Then, the second factor in (3.7)

1 -(1 -π n ) τn 0
λ ν (u)e πnI(z,u,τn) e -hν (u,τn) du , can be written as:

1 -(1 -π n ) τn 0 λ ν (u)e -hν (u,τn) f 1 (z, u, τ n , π n ) du + 1 -e -hν (0,τn) +π n τn 0 I(z, u, τ n )λ ν (u)e -hν (u,τn) du =e -hν (0,τn) - τn 0 λ ν (u)e -hν (u,τn) f 1 (z, u, τ n , π n ) du + π n 1 -e -hν (0,τn) - τn 0 I(z, u, τ n )λ ν (u)e -hν (u,τn) du - τn 0 λ ν (u)e -hν (u,τn) f 1 (z, u, τ n , π n ) du + π 2 n τn 0 I(z, u, τ n )λ ν (u)e -hν (u,τn) du .
Let:

f 3 (z, τ n ) = τn 0 I(z, u, τ n )λ ν (u)e -hν (u,τn) du , and
f 4 (z, τ n , π n ) = τn 0 λ ν (u)e -hν (u,τn) f 1 (z, u, τ n , π n ) du .
Since the limit (4.2) exists and is finite, the limit of f 3 (z, t) as t tends to infinity exists and is finite. Consider now the following term:

e -πnI(z,0,τn) 1 -(1 -π n ) τn 0 λ ν (u)e πnI(z,u,τn) e -hν (u,τn) du .

It can be written as:

(f 2 (z, τ n , π n ) + 1 -π n I(z, 0, τ n )) e -hν (0,τn) -f 4 (z, τ n , π n ) +π n 1 -e -hν (0,τn) -f 3 (z, τ n ) + f 4 (z, τ n , π n ) + π 2 n f 3 (z, τ n ) =(f 2 (z, τ n , π n ) + 1) e -hν (0,τn) -f 4 (z, τ n , π n ) + π n (f 2 (z, τ n , π n ) + 1) 1 -e -hν (0,τn) -f 3 (z, τ n ) + f 4 (z, τ n , π n ) -I(z, 0, τ n ) e -hν (0,τn) -f 4 (z, τ n , π n ) + π 2 n (f 2 (z, τ n , π n ) + 1)f 3 (z, τ n ) -I(z, 0, τ n ) 1 -e -hν (0,τn) -f 3 (z, τ n ) + f 4 (z, τ n , π n ) -π 3 n I(z, 0, τ n )f 3 (z, τ n ) .
Multiplying by e hν (0,τn) :

1 φ(z, τ n ) =(f 2 (z, τ n , π n ) + 1) 1 -e hν (0,τn) f 4 (z, τ n , π n ) + π n (f 2 (z, τ n , π n ) + 1) e hν (0,τn) -1 -e hν (0,τn) f 3 (z, τ n ) +e hν (0,τn) f 4 (z, τ n , π n ) -I(z, 0, τ n ) 1 -e hν (0,τn) f 4 (z, τ n , π n ) + π 2 n (f 2 (z, τ n , π n ) + 1
)e hν (0,τn) f 3 (z, τ n ) -I(z, 0, τ n ) e hν (0,τn) -1 -e hν (0,τn) f 3 (z, τ n ) + e hν (0,τn) f 4 (z, τ n , π n ) -π 3 n I(z, 0, τ n )e hν (0,τn) f 3 (z, τ n ) . Remark now that according to inequality satisfied by f 1 in (4.4):

f 4 (z, τ n , π n ) π 2
n e I(z,0,τn) τn 0 λ ν (u)e -hν (u,τn) du π 2 n e I(z,0,τn) 1 -e -hν (0,τn) π 2 n e hν (0,τn) .

Then, since nπ n e hν (0,τn) tends to α as n tends to infinity, π n e hν (0,τn) tends to 0. Hence:

lim n→+∞ ne hν (0,τn) f 4 (z, τ n , π n ) = 0 .
Denote by o(π n , τ n ) any function such that no(π n , τ n ) tends to 0 as n tends to infinity. Then:

1 φ(z, τ n ) =f 2 + 1 + π n (f 2 + 1) e hν (0,τn) (1 -I(z, τ n )) -1 -I(z, 0, τ n ) + π 2 n (f 2 + 1)e hν (0,τn) I(z, τ n ) -I(z, 0, τ n ) e hν (0,τn) (1 -I(z, τ n )) -1 -π 3 n I(z, 0, τ n )e hν (0,τn) I(z, τ n ) + o(π n , τ n ) =1 + π n e hν (0,τn) 1 -1 -e -hν (0,τn) I(z, τ n ) -1 + o(π n , τ n ) .
Let (φ n ) n∈N be the sequence of functions defined by φ n (z, τ n ) = φ(z, τ n ) n . Then:

φ n (z, τ n ) = exp -n log 1 + π n e hν (0,τn) 1 -1 -e -hν (0,τn) I(z, τ n ) -1 + o(π n , τ n ) .
Since π n n is equivalent to αe -hν,∞ , the following limit is obtained:

lim n→+∞ φ n (z, τ n ) = exp {-m(1 -I ∞ (z))} , with m = α 1 -e -hν ∞ .
Observe that Theorem 4.1 holds whether F µ satisfies (H) or not. As an example, consider again the Haldane model exposed in Section 3: the final instants of normal cells are exponentially distributed with rate λ and the lifetimes of mutant cells are equal to a constant a. The function I(z, t) is then given by

I(z, t) = 1 1 -e -λt t 0 ψ(z, u, t)λe -λ(t-u) du = 1 1 -e -λt i 0 b i (z) t 0 1 [u+ia ; u+(i+1)a) (t)λe -λ(t-u) du = 1 1 -e -λt i 0 b i (z) e -λia -e -λt 1 t∈[ia ; (i+1)a[ + e -λia -e -λ(i+1)a 1 t∈[(i+1)a ; +∞[ ,
where the b i 's are given by (3.8). Hence the limit of I(z, t) as t tends to infinity:

I ∞ (z) = i 0 b i (z)e -λia 1 -e -λa . (4.5)
Remark that for δ = 0 and a = log(2), the result obtained by Ycart [START_REF] Ycart | Fluctuation analysis: can estimates be trusted?[END_REF] is recovered.

Then the pgf of the final number of mutants can be explicited applying Theorem 4.1.

The probabilities of final mutant counts under Haldane model with δ > 0 will be also explicited in Section 6.

Proportional hazard functions

Here, the assumptions on γ and F ν made in previous the section are still valid. Assume that F µ also satisfies (H). Denote by h µ its related function defined by (H 4 ), and by λ µ the associated instantaneous rate. The limit of h µ (0, t) as t tends to infinity will be denoted by h µ,∞ . Moreover, assume there exists ρ > 0 such that for any s 0:

λ ν (s) = ρ(1 -2δ)λ µ (s) . (H ρ )
An interpretation of assumption (H ρ ) was given at the end of Section 3. In this section, a convergence theorem is deduced from Theorem 4.1, and examples are discussed. In particular, the Luria-Delbrück model with cell deaths is recovered.

Consider first the case where:

h µ,∞ = lim t→+∞ h µ (0, t) < +∞ .
Under (H ρ ), the function I is given by (3.10). Then the limit lim t→∞ I(z, 0, t)e -hν (0,t) exists and is finite, and Theorem 4.1 can be applied. Consider now the case where h µ,∞ is infinite. Then:

I(z, 0, t) ∼ t→+∞ ρδ 1 -δ h * µ (0, t) .
Thus: lim t→+∞ e -hν (0,t) I(z, 0, t) = 0 , Then the following result can be deduced from Theorem 4.1:

Theorem 5.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) = α .
Under (H ρ ), as n tends to infinity, the pgf of the number of mutants at time τ n , starting with n normal cells at time 0, converges to

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

I ∞ (z) = 1 1 -e -ρh * µ,∞ 1 e -h * µ,∞ δ(1 -z) + w((1 -δ)z -δ) (1 -δ)(1 -z) + w((1 -δ)z -δ) ρw ρ-1 dw , with h µ,∞ = lim t→+∞ h µ (0, t) , h * µ,∞ = (1 -2δ)h µ,∞ , and m = α 1 -e -ρe -h * µ,∞ .
As a general application of Theorem 5.1, let f be a non negative and increasing function on R + , with finite limit f ∞ as t tends to infinity. Let µ be defined for (s, t) in R 2 + by

h µ (s, t) = log f (t) f (s) ,
Assume that hypothesis (H ρ ) is satisfied. Taking the limit as t tends to infinity of (3.10):

I ∞ (z) = ρ(1 -2δ) 1 -δ log       (1 -2δ) f ∞ f (0) 1-δ (1 -δ) (1 -z) f ∞ f (0) 1-2δ + z -δ      
, and:

I ∞ (z) = 1 1 -f * (0) f * ∞ ρ 1 f * (0) f * ∞ δ(1 -z) + w((1 -δ)z -δ) (1 -δ)(1 -z) + w((1 -δ)z -δ) ρw ρ-1 dw ,
where f * (t) = f 1-2δ for any positive t, and

f * ∞ = f 1-2δ ∞ .
Remark that only the ratio of f ∞ over f (0) has an influence on I ∞ (z). Another natural example is the classic case where the final instants of normal and mutant cells are both exponentially distributed, i.e.:

h µ (s, t) = λ(t -s) , and h ν (s, t) = ρ(1 -2δ)h µ (s, t) ,
where λ is a positive constant. Then the LDD distribution [START_REF] Ycart | Fluctuation analysis with cell deaths[END_REF] is recovered. Actually, if (H ρ ) is satisfied and h µ,∞ = +∞, the LDD distribution can be recovered from Theorem 5.1: Under (H ρ ), assume that h µ,∞ = +∞. As n tends to infinity, the distribution of the number of mutants at time τ n starting with n normal cells at time 0, converges to the distribution with pgf:

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

I ∞ (z) = 1 0 δ(1 -z) + w((1 -δ)z -δ) (1 -δ)(1 -z) + w((1 -δ)z -δ) ρw ρ-1 dw , and m = α 1 -e -ρh * µ,∞ = α .
In other words, the LDD distribution can be extended to the case where F ν (s, •) and F µ (s, •) are non-exponential distributions, as long as F ν (s, •) and F µ (s, •) are cdfs of true measures on R + and the associated hazard functions λ ν and λ µ are proportional.

Calculation algorithm

A probability computation algorithm for the distribution of the final mutant counts is described here under the hypotheses of Theorem 5.1. The pgf ψ is given by

ψ(z, s, t) = k 0 p k (s, t)z k ,
where the p k 's are defined in Proposition 3.3. Thus, (4.1) can be written as:

I(z, t) = k 0 r k (t)z k ,
where r k is defined for any k 0 by

r k (t) = t 0 p k (u, t)λ ν (u)e -hν (u,t) du . Hence: r 0 (t) = 1 1 -e -ρh * µ (0,t) 1 e -h * µ (0,t) δ -δw 1 -δ -δw ρw ρ-1 dw
and for all k > 0:

r k (t) = (1 -δ) k-1 (1 -2δ) 2 1 -e -ρh * µ (0,t) 1 e -h * µ (0,t) (1 -w) k-1 (1 -δ -δw) k+1 ρw ρ dw .
Then I ∞ (z) can be given by

I ∞ (z) = k 0 r k z k ,
where for any k 0:

r k = lim t→+∞ r k (t) .
In other words:

r 0 = 1 1 -e -ρh * µ,∞ 1 e -h * µ,∞ δ -δw 1 -δ -δw ρw ρ-1 dw ,
and for all k > 0:

r k = (1 -δ) k-1 (1 -2δ) 2 1 -e -ρh * µ,∞ 1
e -h * µ,∞

(1 -w) k-1 (1 -δ -δw) k+1 ρw ρ dw .

Moreover, (5.1) admits a series expansion for any z ∈ [0 ; 1]:

φ(z) = k 0 q k z k ,
where the q k 's can be easily expressed as function of the r k 's, using the following algorithm exposed by Embrechts and Hawkes [START_REF] Embrechts | A limit theorem for tails of discrete infinitely divisible laws with applications to fluctuation theory[END_REF]. Firstly: q 0 = φ(0) = e -m(1-r 0 ) .

The derivative of φ with respect to z is given by

dφ dz = m dI ∞ dz φ = m i 1
ir i z i-1 j 0 q j z j = m i 1,j 0 ir i q j z i+j-1 .

On the other hand:

dφ dz = k 1 kq k z k-1 .
Hence for any k > 0:

q k = m k i 1 ,j 0 i+j=k ir i q j = m k k i=1 ir i q k-i .
Naturally, if h µ,∞ = +∞, the probabilities of the LDD distribution [START_REF] Ycart | Fluctuation analysis with cell deaths[END_REF] are recovered.

Consider now the Haldane model. In that case, the pgf I is defined by (4.5). Considering a cell born at time s, let p Hence, the probabilities (r k ) k∈N can be explicited if the probabilities p (i) k k∈N can be explicited for any i 0. In practice, the Fast Fourier Transform can be used to identify the r (i) k 's. Then, the q k 's can be computed using the algorithm of Embrechts and Hawkes [START_REF] Embrechts | A limit theorem for tails of discrete infinitely divisible laws with applications to fluctuation theory[END_REF] described above.

Conclusion and perspectives

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. The main results are based on the theory of supercritical branching processes. It led to a family of distributions, modeling asymptotic number of mutants. These distributions depend on the expected number of mutations m, the death probability of mutant cells δ, and the final instant distributions F ν (s, •) and F µ (s, •) for normal and mutant cells born at a given time s. A convergence theorem for the final count of mutants has been proved for both cases where F ν (s, •) and F µ (s, •) are defined on R + or R + . The first case provides the possibility that a cell does not split or die before the end of the experiment. It enables to model more realistic growths, such as logistic growth. The particular case where the hazard functions λ ν and λ µ associated to F ν (s, •) and F µ (s, •) are proportional has been studied. Computation algorithm for probabilities has been described. Moreover, the LDD distribution is recovered when F ν (s, •) and F µ (s, •) are defined on R + and the associated hazard functions are proportional. The consequences for statistical inference and simulation must be developed. Since the R package flan (available on CRAN: https://cran.r-project.org/package=flan) provides tools for inference of mutation models for the case where final instants are i.i.d., an extension to the model proposed here is planned.

Corollary 5 . 1 .

 51 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) = α .

e

  of the size of its clone in the interval [s + ia ; s + (i + 1)a). In other words:b i (z) = k 0 p (i] k z k ,for any i 0, where the b i 's are given by(3.8). Therefore: -iaλ 1 -e -aλ p (i) k , and the probabilities (r k ) k∈N associated to pgf I ∞ are given by r k = i 0 e -iaλ 1 -e -aλ p (i) k .

  i.d. lifetimes (see for exampleKimmel and Axelrod [15, Chap. 5]). The case of multitype branching processes with non i.i.d. lifetimes is studied here using similar arguments. It leads to equation (BHM). Applying (BHM) to (3.1) gives equation(3.3). Description of this kind of processes with i.i.d. lifetimes can be found in Athreya and Ney[4, chap. 5] or Kimmel and Axelrod[15, chap. 6]. Consider a multitype branching process with a number l + 1 of cell types (incuding the "dead" type, indexed by 0). The model is the following:

Acknowledgements This research was supported by Laboratoire d'Excellence TOUCAN (Toulouse Cancer). The author is grateful to Bernard Ycart for comments on earlier drafts of the paper and to the anonymous referees for helpful remarks.

Proposition 3.2. The general solution of the Riccati equation (R2) is given by

where: h * µ (s, t) = (1 -2δ)h µ (s, t) . The proof of Proposition 3.2 simply consists in observing that 1 is a particular solution of (R2). Then the general solution (R2) is explicit [START_REF] Harko | Analytical Solutions of the Riccati Equation with Coefficients Satisfying Integral or Differential Conditions with Arbitrary Functions[END_REF] and given by (3.9).

Observe that if λ µ ≡ λ with λ a positive constant, then h µ (s, t) = λ(t -s) and Proposition 3.2 reduces to the example of Athreya and Ney [4, p. 109]. Moreover, if λ µ ≡ λ and δ = 0, (R2) reduces to a Bernoulli equation. Its solution is the pgf of the geometric distribution with parameter e -λ(t-s) , i.e. the pgf of a Yule process with parameter λ. (see Yule [34, p. 35] or Athreya and Ney [4, p. 109]). As a direct consequence of Proposition 3.2, the distribution of the number of mutant cells in a mutant clone started at time s can be explicited. Proposition 3.3. Denote by (p k (s, t)) k∈N the probabilities of the size at time t of a clone starting from a single mutant cell born at time s. Then:

and for k 1,

where:

1 -δ -δe -h * µ (s,t) , and:

In other words, a random variable with pgf ψ(z, s, t) is the following random mixture: either 0 with probability p 0 (s, t), or a geometric random variable with parameter P (s, t). Proposition 3.3. Writing ψ as the following rational function:

and

Then the probabilities p k (s, t) can be recursively identified:

and

Denote by P (s, t) the second term of the above product:

Then for k 2:

(1 -p 0 (s))P (s, t)

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

Observe that if λ µ ≡ λ with λ a positive constant, Proposition 3.3 reduces to formula (3.1) of [START_REF] Ycart | Fluctuation analysis with cell deaths[END_REF]. Moreover, the expectation of the size at a finite time t of a mutant clone started at time s is e h * µ (s,t) . Assume there exists ρ > 0 such that for any s 0:

The constant ρ can be interpreted as the instantaneous ratio of hazard functions. The assumption of proportional hazard functions is not new: in survival analysis, it is known as the Cox proportional-hazard regression model [START_REF] Cox | Regression Models and Life-Tables[END_REF], which is widely used. This assumption generalizes the notion of fitness defined in homogeneous mutation models as the ratio of


Emmanuel Jeannot
email: emmanuel.jeannot@inria.fr

Farouk Mansouri
email: farouk.mansouri@inria.fr

Guillaume Mercier
email: guillaume.mercier@bordeaux-inp.fr

A Hierarchical Model to Manage Hardware Topology in MPI Applications

Keywords: Hierarchy, Hardware Topology, Message Passing ACM Reference format:

The MPI standard is a major contribution in the landscape of parallel programming. Since its inception in the mid 90's it has ensured portability and performance for parallel applications on a wide spectrum of machines and architectures. With the advent of multicore machines, understanding and taking into account the underlying physical topology and memory hierarchy have become of paramount importance. The MPI standard in its current state, however, and despite recent evolutions is still unable to offer mechanisms to achieve this. In this paper, we detail several additions to the standard that give the user tools to address the hardware topology and data locality issues while improving application performance.

INTRODUCTION

Parallelizing or writing from scratch a parallel application is a very challenging task and this challenge has become even more bigger due to the current trend in processors design and supercomputers architecture. Indeed, the hardawre that the programmer has to tackle becomes more and more hierarchically organized. For instance, CPUs now feature various levels of memories that have different properties in terms of size, performance ansd even nature. As a consequence, a parallel application performance is likely to be impacted by the communication occuring between processes and by the way they access data. Thus, if a process accesses some data located in a memory bank physically far from the core it currently executes on, a penalty shall occur (NUMA effect). Also, if two processes share a cache level, they will communicate one with the other more efficiently. This is known as the data locality issue.

In order to better exploit the underlying hardware, applications need to take this locality phenomena into account. They need to get a better grasp of the underlying physical architecture they are running on. The current success of hwloc [START_REF] Broquedis | Hwloc: a Generic Framework for Managing Hardware Affinities in HPC Applications[END_REF] demonstrates the need for such a tool and the information it allows to gather. This knowledge may offer the possibility of optimize the code, to better exploit the memory hierarchy or the network topology as well as a global and comprehensive view of the hardware (a continuum of hierarchies between the network and the memory).

To this end, a relevent programming model (along with tools and libraries that implement it) can be of valuable help. The most obvious and natural choice would be to first look at what current parallel programming standards and libraries offer in this area.

In the case of the Message Passing Interface (MPI) [START_REF]MPI: A Message-Passing Interface Standard[END_REF], some mechanisms are already available that could make it possible for an application to exploit the hardware hierarchy to improve interprocess communication and locality. For instance, in the case of virtual topology management routines, such as MPI_Dist_graph_create or MPI_Graph_map, the reorder argument can be used to create a topology where processes are reorganized according to the underlying physical topology of the target architecture.

However, there are several issues with this approach: first, this argument might be used in this fashion, but that is not necessarily the case, which means that the expected behaviour is totally implementation-dependent (and thus not standard) hence likely to change from one implementation to the other or even worse from one implementation version to the other. An application cannot rely on a particular MPI implementation version to be able to use the specific features it needs. Then, in absence of dedicated and relevant mechanisms directly within the standard, an application is forced to use some side-effects of already available features, which constitutes an issue in terms of both interface expressiveness and usability.

In this paper, we present abstraction mechanisms that can help programmers to structure their applications based on physical topology criteria. Such structure can then be used to improve data locality or communication performance by taking into account information that would be otherwise unavailable to the underlying MPI implementation. We also present the implementation of this abstraction and mechanisms in the context of the MPI standard.

This paper is organized as follows: the current status of the MPI interface with regards to hardware topology management is discussed in Section 2 and our proposal is detailed in Section 3. Section results are analyzed in Section 5, while Section 6 concludes this paper and give potential future directions.

HARDWARE TOPOLOGY MANAGEMENT AND THE MPI STANDARD

In this section, we shall examine the current possibilities offered by the MPI standard to deal with the issue of hardware topology management in parallel applications. One key-characteristic of the MPI standard is its hardware-agnosticity. Indeed, it makes no assumptions about the hardware on which the application is going to be deployed and run. This behaviour ensures the portability of parallel programms using MPI. It is to be noted also that despite this independence from any hardware considerations in its programming model, the MPI standard and programming model does not prevent from accessing the hardware topology directly from an application.

Interactions with external tools

One way to tackle this issue is precisely to use an altogether different tool or interface, fully external to MPI. That is, one current practice is to deal with tools representing the hardware topology such as hwloc [START_REF] Broquedis | Hwloc: a Generic Framework for Managing Hardware Affinities in HPC Applications[END_REF], LibNuma [START_REF] Kleen | A NUMA API for Linux[END_REF], or Pthread sched [START_REF] Nichols | Pthreads Programming[END_REF]. These libraries could give a relevent representation of the hardware structure and components. However, they are rather low-level tools and need a good knowledge of the underlying hardware to be used correctly and efficiently. In addition, this practice increases the complexity of developing and supporting codes. Last, as these tools or interfaces are not standard, portability is not guaranteed.

Current status in the MPI standard and its implementations

The MPI library even offers some means to better understand the nature of the physical architecture in order to exploit it to its full potential. For instance, the extension of Remote Memory Operations with Shared Memory operations in the standard has allowed programmers to structure their applications to take into account the fact that processes are located on the same machine. This can be seen as an alternative to the use of multithreading when memory consumption is at stake [START_REF] Reinders | High Performance Parallelism Pearls[END_REF]Chapter 16]. In this particular case, the MPI standard acknowledges that some physical resource (e.g. memory) is shareable between processes and offers the tools to actually access this resource. However, despite the presence of some mechanisms in MPI to better understand and use the hardware, they address the issue only partially. Some MPI implementations offer mechanisms but since they are implementation-dependent they are typically nonstandard. Portability is therefore not guaranteed. The bottom line is: an application should not (or cannot) rely on a specific implementation nor on a specific version of an implementation (which is even worse). On the other hand, improving performances and scalability of applications is more efficient when locality is exposed during their design steps rather than only relying on MPI implementation optimizations. Thus, the MPI programming model as specified by the standard needs to offer high-level abstractions to take into consideration architecture topology at an early stage and before calling implementations features. The following paragraphs describe some of these mechanisms and their shortcomings.

Shared Memory constructs.

MPI_Comm_split_type can accept MPI_COMM_TYPE_SHARED as a value for its split_type argument. The outcome is a communicator that encompasses all processes in the original communicator that can create a shared memory region 1 By using the communicators produced by this function, it is possible to structure an application in order take into account the fact that memory is shared: for instance direct memory accesses are possible instead of relying on message-passing exchanges handled by the MPI library. The overheads of the library in case of shared-memory are eliminated. However, the various levels of cache, a major part of the hardware hierarchy, are left unexploited. As explained previously, this memory hierarchy is becoming more and more complex and performance gains are expected from a relevent exploitation of it.

2.2.2

Process topologies and reordering. The various process topologies available in MPI can help to structure an application but they are virtual topologies and quoting the standard itself: "The virtual topology can be exploited by the system in the assignment of processes to physical processor, if this helps to improve the communication performance on a given machine. How this mapping is done, however, is outside the scope of MPI." Some implementations use the reorder parameter of some functions (e.g. MPI_Graph_map, MPI_Dist_graph_create [START_REF] Träff | Implementing the MPI Process Topology Mechanism[END_REF], etc.) and take the opportunity to retrieve hardware information and make use of it [START_REF] Jeannot | Process Placement in Multicore Clusters: Algorithmic Issues and Practical Techniques[END_REF][START_REF] Mercier | Improving MPI Applications Performance on Multicore Clusters with Rank Reordering[END_REF][START_REF] Rashti | Multi-core and Network Aware MPI Topology Functions[END_REF]. Some others tailor the topology routine of MPI to a specific hardware [START_REF] Hatazaki | Rank Reordering Strategy for MPI Topology Creation Functions[END_REF][START_REF] Larsson | Implementing the MPI process topology mechanism[END_REF]. This, of course, falls into the non-standard category as it is implementation-dependent and is a side-effect of the function. It is not the primary goal of it and the fact that the underlying hardware is efficiently exploited can be seen as a bonus.

Process Managers and process mapping.

Process managers can also be of help when it comes to exploit the underlying hardware. Indeed, through their process mapping and binding options, they can allow the user to finely control the way the various MPI application processes are dispatched and executed [START_REF] Hursey | Locality-Aware Parallel Process Mapping for Multi-core HPC Systems[END_REF]. Thanks to an adequate placement policy enforced by both these mapping and binding parameters, it is possible to take into account the physical topology and reduce the communication costs for instance [START_REF] Brandfass | Rank Reordering for MPI Communication Optimization[END_REF][START_REF] Mercier | Towards an Efficient Process Placement Policy for MPI Applications in Multicore Environments[END_REF]. This is also used to improve collective communication performance [START_REF] Zhang | Process Mapping for MPI Collective Communications[END_REF], Unfortunately, these options are totally non-standard and even change from one version of a process manager to the other. This point is outside the scope of this paper, but standardizing these mapping and binding options/policies would ease the user when it comes to launch their application. However, there are tight interactions between mapping/binding polices and the mechanisms proposed in this paper: in absence of a (relevent) mapping/binding of processes, the performance improvements of taking into account the hardware topology are not expected to be as high as if an efficient policy was to be enforced.

MPI Sessions

. MPI sessions are a new concept that is currently being discussed by the MPI Forum. In its current state, the standard only allows to call MPI_Init and MPI_Finalize a single time in an application. This can raise some issues when multiple libraries that internally rely on MPI are used concurrently. Sessions can be seen a lightweight construct, even lighter than groups. A session encompasses MPI processes and some information can be attached to it, for instance about the application or the hardware. In this case, the sharing of hardware resources could be exploited by the application with several different MPI sessions.

PROPOSED EXTENSIONS TO THE MPI STANDARD

As exposed in the previous section, there are currently no means in the MPI standard to portably take into account the hardware topology at the application level. We think that it is important/necessary to offer high-level abstractions helping programmers to take care of locality and communication optimizations when they design their applications. This way we anticipate and facilitate an implementation work that optimizes communications according to the target architecture. Therefore, we propose to extend the MPI standard and detail in the rest of this section the relevant mechanisms and features needed to achieve these objectives.

Guidelines

Since its first version in 1994, the MPI standard has grown steadily in terms of number of available routines and functionalities. We therefore advocate for a minimal amount of changes and prefer to leverage existing mechanisms. We prefer not introducing new functions unless it is unavoidable and rather expand existing mechanisms. MPI is about communications and how they are managed. In this regard, exploiting the underlying topology boils down to be able to organize the various MPI application processes in a way that is both topologically-wise and performance-wise sensible. The same idea has already been used in the case of reordering: the processes that communicate a lot should be bound on two cores physically close to each other. Consequently, the sharing of caches is likely to decrease communication times and improve overall application performance.

The key idea is therefore to group processes into entities where a specific kind of resource is shared by all group 2 members (i.e. processes). MPI features a concept/construct that perfectly matches our needs: the communicator. As a consequence, we propose to make hardware topology information and structure available at the application level through a hierarchy of communicators. We want to help an application developer and guide him/her to build the most relevant communicator (hardware) topologically-wise, without any deep knowledge of the underlying architecture and regardless of the way the application processes are mapped and/or bound on the machine.

In this hierarchy, each communicator corresponds to a specific resource that is shared by all the processes belonging to it. For instance, if a process shares a L2 cache and a L3 cache with other processes, it will be part of the communicator encompassing all processes sharing the same L2 and part of the communicator encompassing all the processes sharing this level 3 cache. Creating communicators also allows the use of collective communications among processes that share a resource. Collective communication operations are a major feature of MPI. It gives the programmer some ability to structure his/her application and encourages him/her to improve the locality factor of the communications. 2 In the generic sense, we do not deal with the concept of MPI groups here.

Communicators creation

There exists a couple of functions in the MPI standard which create communicators:

• MPI_Comm_create is able to create such an object from an MPI Group. • MPI_Comm_dup duplicates a communicator taken as input argument of the function.

• MPI_Comm_split is able to create subcommunicators from the original one taken as an argument. This function yields k non-overlapping subcommunicators and the partition of the original communicator is determined by the color parameter (which can take MPI_UNDEFINED as a licit value).

Processes that call this function with the same value for color will belong to the same output (sub-)communicator. Since the idea behind our proposal is to create communicators based on the sharing of common resources, splitting some input communicator (MPI_COMM_WORLD being obviously a relevent but not mandatory candidate) is a natural fit for the goal we want to achieve: indeed the information about the sharing of the resource can be conveyed by the color argument.

However, the outcome is not likely to be the one expected: let us take for example the case of processes mapped onto different physical nodes (for the sake of simplicity, each core of each node executes its own process). Let us then assume that each node features several L3 caches. If we want to create as many communicators as the number of L3 caches in our configuration, we cannot provide a single color value, otherwise, all processes would end up belonging to the same communicator. Ideally, we would like to provide the same color value, but this value has to carry a different meaning in different processes.

MPI_Comm_split_type extension. Fortunately, this is exactly the behaviour of the MPI_Comm_split_type function.

MPI_Comm_split_type is part of the MPI-3 shared memory programming functionnalities and features the following prototype This function partitions the group associated with comm into disjoint subgroups, based on the type specified by the value assigned to the split_type parameter. A single value is currently defined in MPI 3.1: MPI_COMM_TYPE_SHARED. when this value is used, the input communicator is split into communicators, where each new communicator represents a shared-memory domain. That is, two processes belonging to the same subcommunicator are able to create a mutualy accessible shared-memory region. Obviously, there is no overlap between these new communicators. This function, along with the particular MPI_COMM_TYPE_SHARED value for the split_type parameter, already allows the user to better understand the way the processes are mapped onto the underlying hardware. It also gives the opportunity to take advantage of it since a different programming model (MPI-3 Shared Memory style) can be used in each new communicator and classical Message Passing between them.

The MPI standard stipulates that implementations may define their own values for the split_type parameter, in order to enforce specific behaviours. The flexibility granted by this approach is counterbalanced by its sheer lack of portability.

As a consequence, we propose to enrich the set of possible values for the split_type argument by adding a new one (MPI_COMM_TYPE _PHYSICAL_TOPOLOGY for instance 4 ). The info argument can be used to pass hints to the implementation about the way the split should be done.

Proprieties of the hierarchical communicators. A call to

MPI_Comm_split_type with this new value shall yield a communicator corresponding to the highest possible level in the hierarchy tree of the hardware topology. This newly produced communicator can then be used as an input argument in subsequent calls to MPI_Comm_split_type to produce other children communicators that correspond to deeper levels (see the example below for a practical use). Also, the newly produced communicators should retain the following proprieties:

• The last valid communicator produced in this fashion may be identical to MPI_COMM_SELF, but not necessarily. • Each recursively created new communicator should be a strict subset of the input communicator. That is, a call to MPI_Comm_compare(oldcomm,newcomm) should return MPI_UNEQUAL. This propriety ensures that we do not create unnecessary new communicators in case of redundancies of levels in the hardware topology. For instance, if a L3 cache and a L2 cache is shared between all processes, there is no need to create a communicator for both resources. • If no valid communicator is to be created, MPI_COMM_NULL should be returned.

These communicators calls will form a kind of hierarchy, mimicking the hardware one, as all new communicators are encompassed (so to speak) in their parent communicator. That is, if a process belongs to the communicator corresponding to the n-th level of the hierarchy, it also belongs to all communicators corresponding to levels 0 to n -1. It is to be noted that our abstraction does not make any distinction between the network and the nodes internal memory hierarchy. We give some means to organize an application according to the structure of the hardware, not its nature.

Creation of Roots Communicators.

One other usefull addition would be the ability to create at the qame time at each level of the hierarchy another communicator which includes all root processes of a hierarchical communicator. This forms another kind of hierarchy of its own and could ease the communication between all the levels of the hierarchy. This function could have the following prototype: 4 Or any suitable and meaningful name.

Interaction with process mapping/binding policies.

A special care should be taken regarding the current binding of process ranks. Indeed, the deepest level that shall be returned should correspond to the current process binding (e.g if a rank is bound to a L3 cache, no information below this level should be returned since it may use different L2 caches below when moving inside the binding. Any attempt to do so should return MPI_COMM_NULL. Moreover, unbound processes may move across an entire sharedmemory machine and therefore cannot belong to a communicator split deeper than the "machine" level: in such a case, the returned communicator shall be the same as the one returned by a call to MPI_Comm_split_type with the MPI_COMM_TYPE_SHARED value for the split_type parameter. 3.2.5 A practical example. We now detail a practical example of use of this new split_type value. Let us suppose that an MPI application is launched on several machines featuring the memory hierarchy depicted by Figure 1: each node is composed of two NumNode with a single socket (package) and 4 cores per socket. Each socket features its dedicated L3 cache and a L2 is shared between a pair of cores. At some point, we assume that the various processes of the application execute the following code: In this code snippet, NLEVELS is chosen approprietly so that there are enough elements in the newcomm array, but as we shall discuss later (see Section 3.4), our proposal does not make any assumption on the total number of levels in the hardware hierarchy nor on the nature of the hardware resource the communicator is supposed to represent for the processes.

A simple case. In our first case, we suppose that each process is bound to its own core. More precisely, process rank p i (in MPI_COMM_ WORLD) is bound to core number i. When the code shown previously is executed, the following communicators are created:

• As for the roots communicators we shall have the following situation (we suppose that rootscomm[i] corresponds to the hierarchy level newcomm[i]):

• rootscomm

[0] = {p 0 } • rootscomm[1] = {p 0 , p 4 } • rootscomm[2] = {p 0 , p 2 } and {p 4 , p 6 } • rootscomm[3] = {p 0 , p 1 }, {p 2 , p 3 }, {p 4 , p 5 } and {p 6 , p 7 }
A slightly more complicated case. Let us now suppose that the processes are not bound to a specific core: some processes are bound to cores, some others are bound on the L2 cache and some on a Nu-maNode. Figure 2 depicts such a case involving 8 processes. P0 and P1 are bound to their own cores, P2 and P3 are bound to a L2 cache, while P4, P5, P6 and P7 are bound to a NumaNode/Package/L3 cache. In this case, and according to the proprieties of the hierarchical communicators:

• The roots communicators would then be: This demonstrates that our proposal is flexible enough to accommodate the case of non-uniform binding policies within the same MPI application.

newcomm[0] = {p 0 , p 1 , p 2 , p 3 , p 4 , p 5 , p 6 , p 7 } • newcomm[1] = {p 0 , p 1 , p 2 ,
• rootscomm[0] = {p 0 } • rootscomm[1] = {p 0 , p 4 } • rootscomm[2] = {p 0 , p 2 }

Communicators characteristics query

The main aspect of our proposal deals with hierarchical communicators creation. However, we need to introduce other functions in order to make the use of these communicators more practical to the application developer. So far, with the proposed abstraction, users are able to leverage the structure of their hardware. However, more information might be needed, for instance in case of data distribution between the various communicators created at a certain level in the hierarchy.

Getting information for a hierarchical level.

A process is able to retrieve some information about a specific communicator it belongs to with a call to the follwing function: With:

IN comm: communicator (handle) OUT num_comms: number of siblings communicators (integer) OUT index: communicator index (integer) OUT type: type of communicator (string)

• num_comms is the number of communicators in the same level with the same parent communicator. For instance, in the case of a node similar to the one shown on Figure 1, a call to MPI_Comm_get_ hlevel_info with newcomm [START_REF] Broquedis | Hwloc: a Generic Framework for Managing Hardware Affinities in HPC Applications[END_REF] (that is, a L2 cache) would yield a value of 2 for num_comms, since there are 2 L2 caches per NumaNode (i.e. newcomm [START_REF] Brandfass | Rank Reordering for MPI Communication Optimization[END_REF], the parent communicator). • index is a kind of rank for each communicator which should be contiguously numbered and starting from 0. It is the rank of the communicator among all communicators created by its parent communicator.

• type is a string giving information about the kind of resource that the communicator represents. It should be unambiguous, like L2_Cache, L3_Cache or NumaNode.

All this information should be cached by the communicator in an info object attached to it containing a set of (key,value) proprely defined. This would require the use of the MPI_Comm_set_info and MPI_Comm_get_info functions.

3.3.2

Getting the minimal level. Another helpful feature would be the ability for a programmer to know the name (type) of the lowest level in the hardware hierarchy that is shared by some processes. To this end, we propose to add the following function: This function returns the name of the lowest level in the hierarchy shared by all the MPI processes which ranks in the communicator comm are listed in the rank array. If the calling process rank is not among the ranks listed in the array passed as an argument, the type returned should be "Unknown" or "Invalid".

Discussion

In this section, we discuss the design choices, advantages and potential drawbacks of our approach. We propose an abstraction based on existing MPI objects (the communicators), hierarchically modelling the hardware topology in order to improve performance and scalability. Communicators are often used in MPI applications and a well-understood concept to boot. A large majority of application developers are familiar with it (beyond MPI_COMM_WORLD). Therefore, using our proposal would require little effort, conceptually speaking. The communicators created do not feature a predetermined name, taken after the underlying resource it is supposed to represent. This ensures that there will no need to make any change nor modifications in the future in case new levels in the hardware hierarchy should appear. This also justifies why we do not specify a maximal depth for the hierarchy. By creating "recursively" new communicators, the user is able to get all needed objects until the bottom is reached and no new communicator can be produced. The user can of course chose the desirable depth by querying the name/type of the communicator created and deciding to go further or not. This choice also explains why we do not rely on several predetermined values for the split_type argument 5 (e.g. one value corresponding to a specific hierarchy level) because architectural changes in hardware would require modifications to the standard.

Currently, our design is based on a "recursive" call to MPI_Comm_ split_type in order to create the hierarchy. However, an alternate solution would be to introduce a function that create all communicators at once, returning an array for instance as well as its size (i.e. the hierarchy depth). Such a function is more simple to use, but 5 As implemented in Open MPI with the OMPI_COMM_TYPE_* values for instance. compells the user to create all the levels, even undesired ones. The chosen design allows for more control at the cost of some ease of use.

One drawback of our approach is the fact that it targets architectures which are hierarchically organized. It is the case for most machines, but there are some exceptions that our model does not currently address. And the same issue arises for network topologies that are not hierarchical, such as torus for instance. In such a case, should the split be made on just a particular dimension of the torus? The info argument could be used to pass this information. But if some form of hierarchy can be extracted from the network topology, we shall be able to exploit it with our mechanism.

Currently, the interface proposed features one new value for the split_type argument of the MPI_Comm_split_type function and three new routines. But one open question is whether or not we should offer more elements. Indeed, it could be interesting to create several hierarchies of communicators, based on different criteria with different split_type values. For instance, making a distinction between the network topology and the nodes internal topology could be a relevant idea. The info argument in MPI_Comm_split_type could also be used to split a communicator directly at a desired level, without creating the whole hierarchy.

As for the names (types) returned by the MPI_Comm_get_hlevel_ info function, they could be given or derived from the names used by an external tools such as hwloc. By doing so, portability could be enforced to some degree, or at least across MPI implementations that feature hwloc (which is the case of Open MPI, MPICH and MVAPICH for instance).

Last, its is possible to create the roots communicators with the current functions available in MPI. However, from a performance standpoint, an hwloc-based implementation (and thus a new function in the standard) is more efficient. It also removes the burden of implementing it by the user.

BENEFITS OF THE HIERARCHICAL MODEL FOR PARALLEL APPLICATIONS

In this section we propose to model the class of MPI applications that can leverage our hierarchical abstraction to enhance the performance of their communications. Let us suppose that an MPI If the communication operations C include collective operations, then the application can take advantage of the hierarchical model described by Algorithm 2, where L is the number of levels in the hierarchy allowed by the target architecture and generated automatically by our approach. Figure 4 illustrates the application communication model based on a three-level communicator hierarchy. Applying a hierarchical model can improve the communication times and thus the total time of the application. As a matter of fact, by exploiting more parallelism in the hierarchy levels, it is possible to enhance the scalability and performance of collective operations. In order to illustrate this claim we propose a short discussion to highlight the important features which can explain this benefit.

The collective communications are modeled with well-known works such as the Hockney [START_REF] Hockney | The Communication Challenge for MPP: Intel Paragon and Meiko CS-2[END_REF], the LogP [START_REF] Culler | LogP: Towards a Realistic Model of Parallel Computation[END_REF] or the PlogP [START_REF] Kielmann | Fast Measurement of LogP Parameters for Message Passing Platforms[END_REF] model. These models express the maximum time taken by a collective communication operation as a linear f unction depending on several variables such as:

• p: The number of processes performing the collective communication operation c. • m: The size of the data exchanged between processes.

• The collective algorithm used.

• α: The hardware latency that may be a function of m and p. • β: The hardware bandwidth that may be a function of m and p. For example, the processing time of the simple Flat Tree broadcast is expressed with the Hockney [START_REF] Hockney | The Communication Challenge for MPP: Intel Paragon and Meiko CS-2[END_REF] model as equal to :

(p -1)(α + mβ)

As shown by Algorithm 2 L levels of communications are present in the topology-aware hierarchical model according to the structure of the architecture. Thus, the time of collective operations will be represented as a sum of L functions depending on the previous cited features but where each level is characterised by a set of specific parameters: p i , α i , β i and m i . As a consequence, the Flat Tree broadcast example shall be expressed as:

L-1 i=0 (p i -1)(α i + mβ i )
The hardawre topology-aware hierarchical model could reduce the processing time of collective operations by exploiting the following points:

• Less number of processes per level: in most applications, it is enough to use a pyramidal structure in the hierarchical model of computation. In fact, as illustrated by Figure 4, the structure based only on the roots of hierarchical communicators in upper levels is enough to perform the communication. such a structure implies to use fewer processes p i by level, which reduces the processing time of the collective. • Parallelism by level: this point regards the exploitation of parallelism at each level of the hierarchy. Indeed, and except for the top level of the hierarchy, the communications inside the communicators of the same level are carried out in parallel. • Data locality and process affinity: This point highlights the advantage of using hardware-aware communicators. Indeed, taking into account the hardware affinity of processes placement in the machine, the communication between them is enhanced thanks to the cache optimizations. • Improved latency and bandwidth per level: this point concerns the physical latency and bandwidth which could be improved when the number of processes is small. In fact, because of the contention phenomenon, the higher the number of processes involved in the communication at the same time through an interconnect, the higher the latency and the lower the bandwidth delivered to each process. From all the above points, collective communication operations could be improved depending on the execution conditions. For instance, if we take the Flat Tree broadcast of a single message (m = 1) and operated in the simple context of 8 processes and the three-level hierarchy as described by Figure 4 then the result is:

The simple case: Time = (8 -1)(α + β)

The hierarchical case:

Time = (2 -1)(α 0 + β 0 ) + (2 -1)(α 1 + β 1 ) + (2 -1)(α 0 + β 1 )
If we compare both expressions, it is clear that the hierarchical approach enhances the simple collective with the factor of at least (7/3) times when latencies and bandwidths are equal in both the simple and the hierarchical cases. However, in the real situation, the simple execution could generate more contention on the interconnects of each level than the hierarchical execution. Thus, the factor of enhancement of collective time could be greater than (7/3).

EXPERIMENTAL RESULTS

In this section we present the experimentations we carried out their results to demonstrate the benefits of our proposal. We modified two benchmarks by introducing our hierarchical model of communicators and we compared their executions on three architectures. In all cases we used the roots communicators hierarchy generated by calling the primitive described in subsection 3.2.3. In addition, we used all cores of a targeted architecture and bound an MPI process on each of them.

Platforms and architectures

For our study, we used tree architectures: a network of 10 nodes (NTW10E5) and two SMP machines (SMPE12E5, SMP20E5) from the Plafrim platform [START_REF] Plafrim | Plate-forme Fédérative pour[END_REF]. The characteristics of these architectures are given in Table 1.

Collective Communications

The first benchmark is the broadcast and the reduce collective communication operations which we chose to test the proposed hierarchical communication model. There has been other works that aimed at exploiting hierarchy in the hardware in order to improve collective communications. In [START_REF] Zhu | Hierarchical Collectives in MPICH2[END_REF] a distinction is made between inter-node and intra-node communication, because shared-memory based communications are expected to be faster than their networkbased counterparts. In this case, the hierarchy was limited to two levels (intra vs. inter-node) and sometimes three (intra. vs intercluster). We generalized this approach and make no assumption about the number of levels. Moreover we are able to exploit the memory hierarchy inside the nodes of a cluster which is not addressed at all by these works. However, the aim here is not to rewrite collectives or to propose some new algorithms. Our goal is merely to experiment our abstraction in order to assess the potential gains achievable by optimizing communication and data locality. We compared two implementations:

• Native: this is the Open MPI implementation of the considered collective. • Hierarchical: this is a loop over the levels of the hierarchy calling the Open MPI version of the collective. Hence, we do not rewrite the collective but just call it through our hierarchy.

Figures 5, 6 and 7 show the enhancement factor of the hierarchicalbased implementation of six broadcast algorithms: Linear, Chain, Pipeline, Split Binary, Binary and Binomial. This factor is obtained by comparing the maximum time of processing several data sizes on the architectures described in Table 1 with the hierarchical broadcast and the native version of Open MPI broadcast. It is possible to note that the hierarchical approach we propose enhances almost all broadcast executions on the three architectures. In fact, the maxima achieved are roughly equal to 21x on NTW10E5, 11x on SMP20E5 and 22x on SMP12E5. This performance gain is due to two majors factors: first, our hierarchical approach allows to better exploit the parallelism and to reduce the complexity of broadcast algorithms. Second, the hardware topology-based communicators enhance the data locality and the hardware affinity between processes performing communications. Third, the communications are better pipelined over the network and the interconnect which reduce their total time.

Figures 8,9 and 10 present the same results for the Reduce collective with six algorithms: Linear, Chain, Pipeline, Binary, Binomial and In-order Binary. In this case also the enhancement factor is obtained by comparing the maximum time for executing the reduction collective of several data sizes on the architectures described in Table 1 with both hierarchical and native Open MPI reduce. Here, we note that the hierarchical approach we propose considerably enhances the first three algorithms: Linear, Chain, Pipeline. In fact, the maxima achieved are roughly equal to 39x on NTW10E5, 11x on SMP20E5 and 8x on SMP12E5. These performances are achieved by the hardware-aware and hierarchical decomposition of the communications. Indeed, the algorithms are better parallelized, the datalocality is enhanced and the communications are better pipelined. However, we note that the three last algorithms (e.g. Binary, Binomial and In-order Binary) are less improved than the others if not at all. This phenomenon is due to the structure of the algorithms since they are already based on a hierarchical, tree-based structure. Therefore, the algorithm complexity is not enhanced by our approach. The small performances obtained for these algorithms are only due to the hardware optimisations and could be more significant with a larger number of processes.

Hierarchical Matrix Multiplication

The second benchmark we carried out is a two-dimensional, hierarchical matrix multiplication [START_REF] Quintin | Hierarchical Parallel Matrix Multiplication on Large-Scale Distributed Memory Platforms[END_REF]. This application is an hierarchical extension of the SUMMA [START_REF] Van De Geijn | SUMMA: scalable universal matrix multiplication algorithm[END_REF] algorithm based on two levels of hierarchy. Its implementation is hardware oblivious and is based on exploring several hierarchy configurations i.e. different number of communicators on each two levels. Hence, the programmer need to manually specify the configuration at each execution by giving a configuration file as an input argument. With our approach, we abstract the manual specification of the hierarchy configuration. As a consequence we do not need to specify the configuration of the hierarchy in our implementation. It is based on the targeted architecture topology and automatically generated by using our proposed set of functions. The implementations we compare are:

• Simple: this is the implementation of the SUMMA algorithm based on broadcasting the blocks over rows and columns. • Hierarchical (xgroups): this is the hierarchical implementation of the SUMMA algorithm of matrix multiplication. The broadcasts of blocks are performed hierarchically over two levels. The hierarchy configuration is set by the user through a configuration file for each execution. • Hierarchical topological: this is the same implementation but leveraging our communicators hierarchy. The hierarchy configuration is based on the underlying hardware topology and automatically build using our primitives. We note that the Topological implementation (last bar) represent the minimum total time and achieve the better speed-ups of the simple implementation: 12x on NTW10E5, 5x on SMP20E7 and 6x on SMP12E5. Indeed, all implementations feature the same processing time (the blue part) but the communication times (the red part) is optimally reduced by the last implementation (topological). This is due to the optimal hardware configuration matching with the topology. In fact, the communications inside hardware topology-aware communicators is more efficient thanks to the data locality and the MPI processes affinity. In addition, pipelining the communication enhance processors occupation time and reduce the global time.

CONCLUSION AND FUTURE WORK

In this paper, we have presented an abstraction that can help the programmer to structure his/her application in order to take into account the hardware topology while he is designing it. We also presented how this model and abstraction could fit into an existing programming model/standard: the widely-used Message Passing Interface. At the expense of light changes and a few new features introduction, we have shown that performance improvements can be achieved in MPI implementations themselves, but more importantly in parallel applications directly. Indeed, we introduced our topology-based communicators in two collective communication operations (Broadcast and Reduce) and showed that in the cases where hierarchy could improve performance, taking into account the physical topology improves things even further. Tested on three different architectures, our approach enhances the performance of the tested collective communications by factors up to 21x, 11x and 22x for the broadcast and 39x, 11x and 8x for the reduce. It also reduces the total time of hierarchical matrix multiplication with the factors of 12x, 5x and 6x. Even if these gains of performance are considerable, the additional effort to use our approach is negligible. In fact, thanks to the proposed abstraction automatically generating an hierarchy of communicators, the user does not have to deal with the details of hardware characteristics or manipulating low-level tools. He only deals with high-level MPI objects: communicators.

The hardware-agnosticity nature of MPI might seem paradoxical with our goal, as we precisely seek to offer the programmer means to better understand and exploit the underlying hardware. We believe that the independence from hardware considerations is an important strength of the MPI standard. We intend to keep MPI hardwareagnostic but we also believe that giving the programmer more hints about this same hardware can be very beneficial performance-wise. The predicament is therefore to find the relevant level of abstraction for such a new functionality. Indeed, if an MPI application should be able to gather and use specific pieces of information about the hardware, this information should nevertheless be abstract enough to not be tailored for a particular class of hardware. This work is available as an external library that features all the functions presented in this paper. The code is downloadable at the following URL: https: //gforge.inria.fr/frs/download.php/file/36832/hsplit-rc-0.1.tar.gz or via a git clone from URL: https://scm.gforge.inria.fr/anonscm/git/ mpi-topology/mpi-topology.git We plan to further discuss this proposal at the MPI Forum. Our future works include the management of networks in our implementation. Currently, we only address a node internal topology (i.e. the memory hierarchy) and we do not create communicators corresponding the the various switch levels (for instance) that are present. This is only a matter of implementation as the concept we propose is ready yet for this aspect. As we use hwloc for memory hierarchies, we would like to support networks by integrating the Netloc [START_REF] Goglin | Netloc: Towards a Comprehensive View of the HPC System Topology[END_REF] software in our implementation. We also plan to address non-hierarchical topologies, especially regarding the network. Expressing the topology with a distance function seems a promising idea and we would like to explore it. Last, we believe that other objects in MPI implementations besides communicators could benefit from a knowledge of the underlying physical topology, for instance memory windows. Generalizing our approach could be of interest.

Figure 1 :

 1 Figure 1: A hierarchical node example.

  newcomm[0] = {p 0 , p 1 , p 2 , p 3 , p 4 , p 5 , p 6 , p 7 } • newcomm[1] = {p 0 , p 1 , p 2 , p 3 } (for NumaNode 0) and newcomm[1] = {p 4 , p 5 , p 6 , p 7 } (for NumaNode 1) • newcomm[2] = {p 2i , p 2i+1 )} (for each L2 number i) : 4 communicators • newcomm[3] = {p i } (for each Core number i) : 8 communicators

p 3 }

 3 (for NumaNode 0) and newcomm[1] = {p 4 , p 5 , p 6 , p 7 } (for NumaNode 1) • newcomm[2] = {p 0 , p 1 } (for L2 0) and newcomm[2] = {p 2 , p 3 } (for L2 1) • newcomm[3] = {p 0 } (for Core 0) and newcomm[3] = {p 1 } (for Core 1)

Figure 2 :

 2 Figure 2: A case of non-uniform binding policy for processes

Figure 3 :Algorithm 1 : 2 MPI

 312 Figure 3: Superstep of MPI application familly based on native collective communications

Figure 4 :

 4 Figure 4: Superstep of MPI application familly based on hierarchical collective communications

Algorithm 2 : 2 L 5 MPI _Collective_Communication_Operations(l); 6 MPI

 2256 MPI algorithm based on a hierarchical communication model 1 MPI _Initialization();(// Initialisation part) while L MPI _Comm_NU LL do ← MPI _Comm_hsplit_with_roots(); // Topological Communicators creation 3 foreach s← S -1..0 do // Algorithm super-steps 4 foreach l← L -1..0 do // Hierarchy levels _Processinд_Operations();

Figure 6 :Figure 7 :Figure 8 :Figure 9 :

 6789 Figure 6: Enhancement factor of hierarchical approach for Open MPI Broadcast implementations SMP20E7 (160 cores)

Figure 10 :Figure 11 :

 1011 Figure 10: Enhancement factor of hierarchical approach for Open MPI Reduce implementations on SMP12E5 (96 cores)

Figure 12 :Figure 13 :

 1213 Figure 12: Comparison of execution time of simple, topology-oblivious hierarchical and toplogical implementations on NTW10E5. 160 processes process 1 block of 8 doubles

  newcomm is the same communicator created by a call to MPI_Comm_ split_type with MPI_COMM_TYPE_PHYSICAL_TOPOLOGY as value for the split_type parameter. rootscomm is the communicator containing all processes that are roots in newcomm. A valid roots communicator can only be returned if the root process of oldcomm calls this function. MPI_COMM_NULL is otherwise returned by non-root processes.

	int MPI_Comm_hsplit_with_roots(MPI_Comm oldcomm,
	MPI_Info info,
	MPI_Comm *newcomm,
	MPI_Comm *rootscomm)
	With:
	IN oldcomm: communicator (handle)
	IN info:info argument (handle)
	IN newcomm: communicator (handle)
	OUT rootscomm: communicator (handle)

Table 1 :

 1 Characteristics of the used architectures

	Name	SMP12E5	SMP20E7	NTW10E5
	OS	Red Hat 4.8	SUSE Srv 11 Red Hat 4.8
	Kernel	3.10.0	2.6.32.46	3.10.0
	Nodes	1	1	10
	Cores per NUMA	8	8	6
	NUMA nodes	12	20	4
	Sockets	12	20	2
	NUMA groups	12	20	1
	Socket	E5-4620	E7-8837	E5-4620
	Clock rate	2600Mhz	2660Mhz	2600Mhz
	Hyper-Threading	Yes	No	No
	L1 cache	32K	32K	32K
	L2 cache	256K	32K	256K
	L3 cache	20480K	24576K	15360K
	Mem Interconnect NUMAlink6 NUMAlink5 QPI
	Node Interconnect N/A	N/A	InfiniBand
	Hierarchical levels 3	3	4
	GCC	5.1	5.1	5.1
	Open MPI	2.0.1	2.0.1	2.0.1
	Hwloc	2.0-git	2.0-git	2.0-git

It is obviously the case when they are located on the same machine, but could also be the case if the processes are on different machines linked by a network à la SCI.

Only the C version is shown.

Figure 5: Enhancement factor of hierarchical approach for Open MPI Broadcast implementations on NTW10E5 (240cores)

ACKNOWLEDGMENTS

This work has partially been supported by the PIA ELCI project of the French FSN and by the ANR MOEBUS project ANR-13-INFR-0001. Experiments presented in this paper were carried out using the PLAFRIM experimental testbed, being developed under the Inria PlaFRIM development action with support from Bordeaux INP, LABRI and IMB and other entities: Conseil Régional d'Aquitaine, Université de Bordeaux, CNRS and ANR in accordance to the programme d'Investissements d'Avenir (see https://www.plafrim.fr/). The authors would like to thank the MPI Forum for its feedback, especially Daniel Holmes.