Pauline Lafitte-Godillon

Kilian Raschel
email: kilian.raschel@lmpt.univ-tours.fr

Viet Chi Tran

Extinction probabilities for a distylous plant population modeled by an inhomogeneous random walk on the positive quadrant

Keywords: Inhomogeneous random walk on the positive quadrant, boundary absorption, transport equation, method of characteristics, self-incompatibility in flower populations, extinction in diploid population with sexual reproduction AMS: 60G50, 60J80, 35Q92, 92D25

In this paper, we study a flower population in which self-reproduction is not permitted. Individuals are diploid, that is, each cell contains two sets of chromosomes, and distylous, that is, two alleles, A and a, can be found at the considered locus S. Pollen and ovules of flowers with the same genotype at locus S cannot mate. This prevents the pollen of a given flower to fecundate its own stigmata. Only genotypes AA and Aa can be maintained in the population, so that the latter can be described by a random walk in the positive quadrant whose components are the number of individuals of each genotype. This random walk is not homogeneous and its transitions depend on the location of the process. We are interested in the computation of the extinction probabilities, as extinction happens when one of the axis is reached by the process. These extinction probabilities, which depend on the initial condition, satisfy a doubly-indexed recurrence equation that cannot be solved directly. Our contribution is twofold : on the one hand, we obtain an explicit, though intricate, solution through the study of the PDE solved by the associated generating function. On the other hand, we provide numerical results comparing stochastic and deterministic approximations of the extinction probabilities.

Introduction

We consider the model of flower population without pollen limitation introduced in Billiard and Tran [START_REF] Billiard | A general stochastic model for sporophytic self-incompatibility[END_REF]. The flower reproduction is sexual: plants produce pollen that may fecundate the stigmata of other plants. We are interested in self-incompatible reproduction, where an individual can reproduce only with compatible partners. In particular, self-incompatible reproduction prevents the fecundation of a plant's stigmata by its own pollen. Each plant is diploid and characterized by the two alleles that it carries at the locus S, which decide on the possible types of partners with whom the plant may reproduce (as it encodes the recognition proteins present on the pollen and stigmata of the plant). We consider the distyle case with only two possible types for the alleles, A or a. The plants thus have genotypes AA, Aa or aa. The only interesting case is when A is dominant over a (see [START_REF] Billiard | A general stochastic model for sporophytic self-incompatibility[END_REF]), and we restrict to this case in this work. Then, the phenotype, i.e. the type of proteins carried by the pollen and stigmata, of individuals with genotypes AA (resp. Aa and aa) is A (resp. A and a). Only pollen and stigmata with different proteins can give viable seeds, i.e. pollen of a plant of phenotype A can only fecundate stigmata of a plant of phenotype a and vice-versa. It can be seen that seeds AA cannot be created, since the genotype of individuals of phenotype a is necessarily aa that combine only with individuals of phenotype A that have genotypes AA or Aa, therefore we can consider without restriction populations consisting only of individuals of genotypes Aa and aa. Each viable seed is then necessarily of genotype Aa or aa with probability 1/2. It is assumed that ovules are produced in continuous time at rate r > 0 and that each ovule is fecundated to give a seed, provided there exists compatible pollen in the population. The lifetime of each individual follows an exponential distribution with mean 1/d, where d > 0. In all the article, we consider

r > d (1.1)
which, we will see, is the interesting case. (b) Transition probabilities of the embedded random walk, that we denote by (X t , Y t ) t∈N (here and throughout, N is the set {0, 1, 2, . . .}), with an abuse of notation.

Let us denote by X t and Y t the number of individuals of genotype Aa (phenotype A) and aa (phenotype a) at time t ∈ R + . The process (X t , Y t ) t∈R + is a pure-jump Markov process with transitions represented in Fig. 1(a). A stochastic differential equation (SDE) representation of (X t , Y t ) t∈R + is given in [START_REF] Billiard | A general stochastic model for sporophytic self-incompatibility[END_REF]. Here we forget the continuous-time process, and we are interested in the embedded discretetime Markov chain, which we denote, with an abuse of notation, by (X t , Y t ) t∈N , and with transitions represented in Fig. 1(b):

P i,j [(X 1 , Y 1 ) = (i -1, j)] = d i (r + d)(i + j) , P i,j [(X 1 , Y 1 ) = (i + 1, j)] = r 2(r + d) , P i,j [(X 1 , Y 1 ) = (i, j -1)] = d j (r + d)(i + j) , P i,j [(X 1 , Y 1 ) = (i, j + 1)] = r 2(r + d)
,

where P i,j means that the process starts with the initial condition (X 0 , Y 0 ) = (i, j). The main and profound difficulty is that this random walk is not homogeneous in space, while techniques developed in the literature for random walks on positive quadrants mostly focus on the homogeneous case (see e.g. Fayolle et al. [START_REF] Fayolle | Random walks in the quarter-plane[END_REF], Klein Haneveld and Pittenger [START_REF] Klein Haneveld | Escape time for a random walk from an orthant[END_REF], Kurkova and Raschel [START_REF] Kurkova | Random walks in Z 2 + with non-zero drift absorbed at the axes[END_REF], Walraevens, van Leeuwaarden and Boxma [START_REF] Walraevens | Power series approximations for two-class generalized processor sharing systems[END_REF]). We introduce a generating function (1.4) that satisfies here a partial differential equation (PDE) of a new type that we solve. Although the particularity of the problem is exploited, these techniques and the links between probability and PDEs may be extended to carry out general studies of inhomogeneous random walks in cones. The introduction of PDEs through generating functions had been already used by Feller [START_REF] Feller | An introduction to probability theory and its applications[END_REF] for a trunking problem with an inhomogeneous random walk in dimension 1. To our knowledge, the case of inhomogeneous random walks in the cone with absorbing boundaries has been left open. In [START_REF] Rege | Queue-Length Distribution for the Discriminatory Processor-Sharing Queue Operations[END_REF], the discriminatory processor-sharing queue is considered but boundaries are not absorbing and the overall arrival rate is constant, which is not the case in our model. When one of the phenotype A or a disappears, reproduction becomes impossible and the extinction of the system occurs. We are interested in the probability of extinction of (X t , Y t ) t∈N (or, equivalently, in that of (X t , Y t ) t∈R + ). Let us introduce the first time at which one of the two types gets extinct:

τ 0 = inf{t ∈ N : X t = 0 or Y t = 0}. (1.2)
For i, j ∈ N, let us denote by

p i,j = P i,j [τ 0 < ∞] (1.3)
the absorption probabilities, and by

P (x, y) = i,j 1 p i,j x i y j (1.4)
their generating function. By symmetry arguments, we have, for all i, j ∈ N,

p i,j = p j,i . (1.5)
Moreover, for any i, j ∈ N such that i = 0 or j = 0, we have

p i,j = 1. (1.6)
In Section 2, we will see that the p i,j 's satisfy the Dirichlet problem associated with the following doubly-indexed recurrence equation

q i,j = di (r + d)(i + j) q i-1,j + dj (r + d)(i + j) q i,j-1 + r 2(r + d) q i,j+1 + r 2(r + d) q i+1,j (1.7)
and with the boundary condition (1.6). This problem does not admit simple solutions. There is no uniqueness of solutions to this problem. Note that the constant sequence equal to 1 is a solution. However, we are interested in solutions that tend to 0 as i or j tends to infinity, since, [START_REF] Billiard | A general stochastic model for sporophytic self-incompatibility[END_REF] (see Proposition 2.2 in this paper), estimates for p i,j were obtained through probabilistic coupling techniques; they show that in the case (1.1) we consider, p i,j is strictly less than 1. In fact, the p i,j 's correspond to the smallest positive solution of the Dirichlet problem, and are completely determined if we give the probabilities (p i,1 ) i 1 . We conclude the section with more precise estimates of the absorption probabilities p i,j as the initial state (i, j) goes to infinity along one axis (Proposition 2.3). These new estimates rely on Proposition 2.2 and on comparisons with one-dimensional random walks. In Section 3, we consider the generating function P (x, y) associated with the p i,j 's and show that it satisfies a PDE, that has one and only one solution, that is computed (Proposition 3.5) explicitly with a dependence on the (p i,1 ) i 1 , prompting us to use the name "Green's function".

Kf (i, j) = (f (i + 1, j) + f (i, j + 1)) r 2(r + d) + f (i -1, j) di (r + d)(i + j) + f (i, j -1) dj (r + d)(i + j) .
Following classical proofs (e.g. [START_REF] Baldi | Martingales and Markov chains[END_REF][START_REF] Revuz | Markov chains[END_REF]), the extinction probabilities (p i,j ) i,j∈N satisfy the equation:

∀i, j ∈ N * , f (i, j) = Kf (i, j) and ∀i, j ∈ N, f (i, 0) = f (0, j) = 1. (2.1)
The constant solution equal to 1 is a solution to (2.1). Let us prove that (p i,j ) i,j∈N is the smallest positive solution to (2.1). Let f be another positive solution. Let us consider M t = f (X inf{t,τ 0 } , Y inf{t,τ 0 } ), with τ 0 defined in (1.2). Denoting by (G t ) t∈N the filtration of (M t ) t∈N , we have:

E M t+1 | G t = E M t+1 1l τ 0 t + M t+1 1l τ 0 >t | G t = E M t 1l τ 0 t + f (X t+1 , Y t+1 )1l τ 0 >t | G t = M t 1l τ 0 t + 1l τ 0 >t E f (X t+1 , Y t+1 ) | G t = M t 1l τ 0 t + 1l τ 0 >t Kf (X t , Y t ) = M t 1l τ 0 t + 1l τ 0 >t f (X t , Y t ) = M t .
Hence (M t ) t∈N is a martingale, which converges on {τ 0 < ∞} to f (X τ 0 , Y τ 0 ) = 1 (see the boundary condition in (2.1)). Thus by using the positivity of f and Fatou's lemma, we obtain that for every i, j ∈ N:

f (i, j) = E i,j M 0 = lim t→∞ E i,j M t E lim inf t→∞ M t 1l τ 0 <∞ = E i,j 1l τ 0 <∞ = p i,j .
This concludes the proof of Point (i).

Let us now consider Point (ii). Assume that the probabilities (p i,1 ) i 1 are given, and let us prove, by recursion, that every p i,j can be computed. By symmetry, we only need to prove that this is the case for i j. Assume (Hrec j): for j ∈ N * all the p k, 's for j and k can be computed from the p i,1 's and let us prove that we can determine the p i,j+1 's for i j + 1. From (1.7) we get:

p i,j+1 = 2(r + d) r p i,j - 2 d i r(i + j) p i-1,j - 2 d j r(i + j) p i,j-1 -p i+1,j . (2.2)
All the terms in the r.h.s. of (2.2) are known by (Hrec j), and hence p i,j+1 can be computed for any i j + 1. This concludes the recursion.

The following result shows that there is almost sure extinction in the case r d. In the interesting case r > d, it also shows that there is a nontrivial solution to the Dirichlet problem (1.6)-(1.7).

Proposition 2.2 (Proposition 9 of [START_REF] Billiard | A general stochastic model for sporophytic self-incompatibility[END_REF]). We have the following regimes given the parameters r and d: (ii) If r > d(> 0), then there is a strictly positive survival probability. Denoting by (i, j) the initial condition, we have:

(i) If r d,
d r i+j p i,j d r i + d r j - d r i+j .
In Point (ii), only bounds, and no explicit formula, are available for the extinction probability p i,j . The purpose of this article is to address (1.7) by considering the Green's function P (x, y) introduced in (1.4).

Asymptotic behavior of the absorption probability as the initial state goes to infinity along one axis

In this part, using the result of Proposition 2.2, we provide more precise estimates of the asymptotic behavior of the absorption probability p 1,j = p j,1 when j → ∞. In particular, these estimates will be very useful when we tackle the deterministic numerical simulations (see Section 4.2).

Proposition 2.3. If j → ∞, then p 1,j = p j,1 = 2d r 1 j - 2d(r 2 + dr + 2d 2 ) r 2 (r + d) 1 j 2 + O 1 j 3 . (2.3)
Proof. In addition to τ 0 , defined in (1.2), we introduce

S = inf{t ∈ N : Y t = 0}, T = inf{t ∈ N : X t = 0},
the hitting times of the horizontal axis and vertical axis, respectively. Note that we have τ 0 = inf{S, T }. Let f : N → N be a function such that f (j) < j for any j 1. In the sequel, we will choose f (j) = j , with ∈ (0, 1) and where . denotes the integer part). We obviously have the identity:

p 1,j = P (1,j) [τ 0 < ∞] = P (1,j) [τ 0 f (j)] + P (1,j) [f (j) < τ 0 < ∞]. (2.4)
To prove Proposition 2.3, we shall give estimates for both terms in the r.h.s. of (2.4).

First step:

Study of P (1,j) [τ 0 f (j)].
Since f (j) < j, it is impossible, starting from (1, j), to reach the horizontal axis before time f (j), and we have

P (1,j) [τ 0 f (j)] = P (1,j) [T f (j)].
In order to compute the latter probability, we introduce two one-dimensional random walks on N, namely X -and X + , which are killed at 0, and which have the jumps

P i [X ± 1 = i -1] = q ± i , P i [X ± 1 = i + 1] = p ± i , P i [X ± 1 = i] = r ± i , q ± i + p ± i + r ± i = 1,
where

q ± i = di (r + d)(i + j ∓ f (j)) , p ± i = r 2(r + d) . (2.5)
Both X -and X + are (inhomogeneous) birth-and-death processes on N. These random walks are implicitely parameterized by j.

If T ± = inf{t ∈ N : X ± t = 0}, then P 1 [T -f (j)] P (1,j) [T f (j)] P 1 [T + f (j)].
(2.6)

The quantities P 1 [T ± f (j)] are computable: we shall prove that

P 1 [T ± f (j)] = 2d r 1 (j ∓ f (j)) - 2d(r 2 + dr + 2d 2 ) r 2 (r + d) 1 (j ∓ f (j)) 2 + O 1 (j ∓ f (j)) 3 .
(2.7)

The main idea for proving (2.7) is that the q ± i being very small as j → ∞, the only paths which will significantly contribute to the probability P 1 [T ± f (j)] are the ones with very few jumps to the right. Let us define Λ ± t (p) = {the chain X ± makes exactly p jumps to the right between 0 and t} = {there exist 0

q 1 < • • • < q p t -1 such that X ± q 1 +1 -X ± q 1 = • • • = X ± qp+1 -X ± qp = 1}.
We are entitled to write

P 1 [T ± f (j)] = P 1 [T ± f (j), Λ ± f (j) (0)] + P 1 [T ± f (j), Λ ± f (j) (1)] + P 1 [T ± f (j), ∪ p 2 Λ ± f (j) (p)], (2.8)
and we now separately analyze the three terms in the right-hand side of (2.8). First:

P 1 [T ± f (j), Λ ± f (j) (0)] = f (j) k=1 P 1 [T ± = k, Λ ± f (j) (0)] = f (j) k=1 (r ± 1 ) k-1 q ± 1 = q ± 1 1 -r ± 1 1 -(r ± 1 ) f (j) = q ± 1 1 -r ± 1 (1 + O((r ± 1 ) f (j) )).
(2.9)

A Taylor expansion of q ± 1 /(1 -r ± 1 ) according to the powers of 1/(j ∓ f (j)) together with the fact that (r ± 1 ) f (j) = o(1/(j ∓ f (j)) 3 ) provides that:

P 1 [T ± f (j), Λ ± f (j) (0)] = 2d r(j ∓ f (j)) 1 - 1 + 2d r j ∓ f (j) + O 1 (j ∓ f (j)) 2 .
(2.10)

We now consider the second term in the right-hand side of (2.8). On the event Λ ± f (j) (1), X ± first stays a time k 1 at 1, then jumps to 2, where it remains k 2 unit of times; it next goes to 1, and, after a time k 3 , jumps to 0. Further, since T ± f (j), we have

k 1 + 1 + k 2 + 1 + k 3 + 1 f (j). Denoting by k 1 = k 1 + k 3 the
time spent in position 1, we thus have:

P 1 [T ± f (j), Λ ± f (j) (1)] = k 1 +k 2 +k 3 f (j)-3 (r ± 1 ) k 1 p ± 1 (r ± 2 ) k 2 q ± 2 (r ± 1 ) k 3 q ± 1 = p ± 1 q ± 1 q ± 2 k 1 +k 2 f (j)-3 (r ± 1 ) k 1 (r ± 2 ) k 2 = p ± 1 q ± 1 q ± 2 (1 -r ± 1 )(1 -r ± 2 ) (1 + O((r ± 1 ∨ r ± 2 ) f (j)-2 )). (2.11)
As for the first term, using the fact that (r

± 1 ∨ r ± 2 ) f (j)-2 = o(1/(j ∓ f (j))
3 ) and a Taylor expansion according to the powers of 1/(j ∓ f (j)) gives that:

P 1 [T ± f (j), Λ ± f (j) (1)] = 4d 2 r(r + d)(j ∓ f (j)) 2 1 + O 1 (j ∓ f (j))
.

(2.12)

Finally, let us consider the third term

P 1 [T ± f (j), ∪ p 2 Λ ± f (j) (p)]. On ∪ p 2 Λ ± f (j) (p)
, the two first jumps to the right are either from 1 to 2 and 2 to 3, or twice from 1 to 2. Thus, extinction means that there is at least 3 jumps from 3 to 2, 2 to 1 and 1 to 0 or two jumps from 2 to 1 and 1 to 0. Since q ± i is an increasing function of i, we deduce:

P 1 [T ± f (j), ∪ p 2 Λ ± f (j) (p)] q ± 3 3 = d r + d 3 3 3 + j ∓ f (j) 3 = O 1 (j ∓ f (j)) 3 . (2.13)
From (2.8), (2.9), (2.11) and (2.13), we obtain (2.7).

Second step:

Study of P (1,j) [f (j) < τ 0 < ∞]. P (1,j) [f (j) < τ 0 < ∞] = k, 1 P (1,j) [f (j) < τ 0 < ∞|(X f (j) , Y f (j) ) = (k, )]P (1,j) [(X f (j) , Y f (j) ) = (k, )] = k, 1 p k, P (1,j) [(X f (j) , Y f (j) ) = (k, )], (2.14)
by using the strong Markov property. Introduce now a function g : N → N such that f (j) + g(j) < j for any j 1. We can split (2.14) into k, g(j)

p k, P (1,j) [(X f (j) , Y f (j) ) = (k, )] + g(j) k 1 and/or g(j) 1 p k, P (1,j) [(X f (j) , Y f (j) ) = (k, )]. ( 2
p k, P (1,j) [(X f (j) , Y f (j) ) = (k, )] 2 d r g(j)
.

In particular, if we choose g such that as j → ∞, g(j) → ∞ fast enough, then clearly the term above is negligible compared to (2.7). For the second sum in (2.15), g(j) k 1 and/or g(j)

1

p k, P (1,j) [(X f (j) , Y f (j) ) = (k, )] = g(j) k 1 p k, P (1,j) [(X f (j) , Y f (j) ) = (k, )],
since by assumption j -f (j) > g(j) so that Y cannot reach values g(j) in f (j) steps. Then we have

g(j) k 1 p k, P (1,j) [(X f (j) , Y f (j) ) = (k, )] g(j) k 1 P (1,j) [(X f (j) , Y f (j) ) = (k, )] P (1,j) [0 X f (j) g(j)].
(2.16)

To obtain an upper bound for (2.16) we are going to use, again, a one-dimensional random walk. Introduce X, a random walk on N which is killed at 0, homogeneous on N * with jumps

P k [ X 1 = k -1] = q, P k [ X 1 = k + 1] = p, P k [ X 1 = k] = r, q + p + r = 1,
where

q = d(1 + f (j)) (r + d)(1 + j -2f (j)) , p = r 2(r + d) .
This walk is again parameterized by j. By construction of X, we have

P (1,j) [0 X f (j) g(j)] P 1 [ X f (j) g(j)].
Denoting by m and σ 2 the mean and the variance of ( X 2 -X 1 ), respectively (they could easily be computed), we can write

P 1 [ X f (j) g(j)] = P 0 X f (j) -mf (j) σ f (j) g(j) -1 -mf (j) σ f (j) .
By a suitable choice of the functions f and g, for instance f (j) = j with ∈ (0, 1) and g(j) = j 3/4 , the central limit theorem gives that the latter is negligible compared to P (1,j) [τ 0 f (j)]. For this last term, using (2.7), (2.6) and letting tend to 0 provides (2.3). The proof is concluded.

Let us make some remarks on possible extensions of Proposition 2.3.

Remark 1. 1. The proof of Proposition 2.3 can easily be extended to the asymptotic of p i,j as j → ∞, for any fixed value of i. In particular, we have the following asymptotic behavior:

p i,j = 2d r i i! j i + O 1 j i+1 .
(2.17)

2. It is possible to generalize (2.8) by

P 1 [T ± f (j)] = k-1 p=0 P 1 [T ± f (j), Λ f (j) p ] + P 1 [T ± f (j), ∪ p k Λ f (j) p ]. ( 2

.18)

We can show as in the proof of Proposition 2. (2.13)) and that the probabilities

3 that P 1 [T ± f (j), ∪ p k Λ f (j) p ] (q ± k+1 ) k+1 = o(1/j k+1 ) (see
P 1 [T ± f (j), Λ f (j)
p ] admit Taylor expansions in powers of 1/(j ∓ f (j)) where the development for the p th probability has a main term in 1/(j ∓ f (j)) p+1 . This allows us to push the developments in (2.7) to higher orders.

For instance, the next term in (2.3) can be obtained by a long computation. First, we generalize (2.11) by writing that

P 1 [T ± f (j), Λ f (j) 2 ] = q ± 1 (p ± 1 ) 2 (q ± 2 ) 2 + p ± 1 p ± 2 q ± 3 q ± 2 k 1 + k 2 + k 3 f (j)-5 (r ± 1 ) k 1 (r ± 2 ) k 2 (r ± 3 ) k 3 (2.19)
where k 1 , k 2 and k 3 are the times spent by the random walk in the states 1, 2 and 3. Then, pushing further the Taylor expansion leads to:

p 1,j = 2d r 1 j ∓ f (j) - 2d(r 2 + dr + 2d 2 ) r 2 (r + d) 1 (j ∓ f (j)) 2 + d 2 r 1 + 2d r 2 - 24d r 2 + d r 2 (r + d) + 5r 2 d 2 2(r + d) 2 r 2 + d 3 1 (j ∓ f (j)) 3 + O 1 (j ∓ f (j)) 4 .

Green's function 3.1 A functional equation for the Green's function

In this section, we consider the Green's function P (x, y) defined in (1.4) associated with a solution of (1.7) in the same spirit as what can be found in Feller [START_REF] Feller | An introduction to probability theory and its applications[END_REF]Ch. XVII]. We show that it satisfies a non-classical linear PDE that can be solved (see Proposition 3.5).

Proposition 3.1. (i) The function P (x, y) satisfies formally:

AP (x, y) = h(x, y, P ), (3.1)
where:

AP (x, y) =Q(x, y) ∂P ∂x (x, y) + Q(y, x) ∂P ∂y (x, y) + R(x, y)P (x, y), (3.2)
Q(x, y) =(r + d)x - r 2 - r 2 x y -d x 2 , R(x, y) = r 2x + r 2y -dx -dy,
and where:

h(x, y, P ) = - r 2 x ∂ 2 P ∂x∂y (x, 0) + y ∂ 2 P ∂y∂x (0, y) + d xy 1 1 -x + 1 1 -y . (3.3)
(ii) For given (p i,1 ) i 1 , we have a unique classical solution to

(3.1)-(3.3) on ]0, 1[×]0, 1[.
The function h in (3.3) only depends on a boundary condition (∂ 2 P/∂x∂y at the boundaries x = 0 or y = 0, i.e. the p i,1 's for i ∈ N * ), which is non-classical, while the operator A is of first order and hence associated with some transport equations.

Proof of Proposition 3.1. Let us first establish (i). Using the Markov property at time t = 1:

p i,j = r 2(r + d) (p i+1,j + p i,j+1 ) + dj (r + d)(i + j) p i,j-1 + di (r + d)(i + j) p i-1,j ,
then multiplying by x i y j , and summing over i, j ∈ N * leads to:

(r + d) i,j 1 (i + j)p i,j x i y j = r 2 i,j 1 (i + j)(p i+1,j + p i,j+1 )x i y j + d i,j 1 jp i,j-1 x i y j + d i,j 1 ip i-1,j x i y j . (3.4)
The l.h.s. of (3.4) equals:

(r + d) i,j 1 (i + j)p i,j x i y j = (r + d) x ∂ ∂x + y ∂ ∂y P (x, y). (3.5)
For the r.h.s. of (3.4):

r 2 j 1 i 2 p i,j (i -1 + j)x i-1 y j + r 2 i 1 j 2 p i,j (i + j -1)x i y j-1 + d i 1 j 0 (j + 1)p i,j x i y j+1 + d i 0 j 1 (i + 1)p i,j x i+1 y j = r 2 j 1 i 2 p i,j (i -1)x i-1 y j + r 2 j 1 i 2 p i,j jx i-1 y j + r 2 i 1 j 2 p i,j ix i y j-1 + r 2 i 1 j 2 p i,j (j -1)x i y j-1 + d i 1 j 0 (j + 1)p i,j x i y j+1 + d i 0 j 1 (i + 1)p i,j x i+1 y j . (3.6)
For the first term in the r.h.s. of (3.6): r 2

j 1 i 2 p i,j (i -1)x i-1 y j = r 2 j 1 i 1 p i,j (i -1)x i-1 y j = r 2 j 1 i 1 p i,j ix i-1 y j - r 2 j 1 i 1 p i,j x i-1 y j = r 2
∂ ∂x P (x, y) - r 2 P (x, y) x . (3.7)
Similar computation holds for the 4th term of the r.h.s. of (3.6). For the 2nd term:

r 2 j 1 i 2 p i,j jx i-1 y j = r 2 y x j 1 i 2 p i,j jx i y j-1 = r 2 y x j 1 i 1 p i,j jx i y j-1 - j 1 p 1,j jxy j-1 = r 2 y x ∂ ∂y P (x, y) - r 2 j 1 p 1,j jy j . (3.8)
We handle the 3rd term of the r.h.s. of (3.6) similarly. Now for the 5th term:

d i 1 j 0
p i,j (j + 1)x i y j+1 =d y 2 ∂ ∂y P (x, y) + yP (x, y) + i 1

x i y .

(3.9)

Similar computation holds for the last term of (3.6). From (3.5), (3.6), (3.7), (3.8) and (3.9) we deduce that:

(r + d) x ∂ ∂x + y ∂ ∂y P (x, y) = r 2 ∂ ∂x P (x, y) - r 2 P (x, y) x + r 2
∂ ∂y P (x, y) - r 2 P (x, y) y + r 2 y x ∂ ∂y P (x, y) - r 2 j 1 p 1,j jy j + r 2 x y ∂ ∂x P (x, y) - r 2 i 1 p i,1 ix i +d y 2 ∂ ∂y P (x, y) + yP (x, y) + i 1 x i y + d x 2 ∂ ∂x P (x, y) + xP (x, y) + j 1
xy j , and finally i 1 p i,1 ix i = x ∂ 2 P ∂x∂y (x, 0). For point (ii), local existence and uniqueness stem from classical theorems [START_REF] Zachmanoglou | Introduction to partial differential equations with applications[END_REF]. Note that we construct an explicit, albeit complicated solution (3.5) using the method of characteristics. This concludes the proof.

Remark 2. In the case of homogeneous random walks, the operator AP (x, y) has the product form R(x, y)P (x, y), see [START_REF] Kurkova | Random walks in Z 2 + with non-zero drift absorbed at the axes[END_REF]. In some sense, this means that the inhomogeneity leads to partial derivatives in the functional equation.

Remark 3. This technique allowing to compute the solution of a discrete, linear problem thanks to generating series is also called the Z-transform method.

Characteristic curves

In Sections 3.2 and 3.3 we establish, by the methods of characteristic equations, an explicit formula for the solutions to (3.1), which proves the existence of the solution. Since A is a first-order differential operator, we have a transport-like PDE. We introduce the following characteristic ordinary differential equations (ODEs). Let (x s , y s ) s∈R + be the solution to the system:

     ẋs = dx ds (s) = Q(x s , y s ), ẏs = dy ds (s) = Q(y s , x s ), (3.10)
where Q has been defined in (3.2). The dynamical system (3.10) and its solutions will turn out to be decisive in the sequel-e.g. in Proposition 3.4, where we will use these characteristic equations in order to express the solutions to the fundamental functional equation (3.1).

Proposition 3.2. For any initial condition (x 0 , y 0 ) ∈ R 2 such that x 0 , y 0 = 0, there exists a unique solution to (3.10), defined for s ∈ R by: with:

λ = 2d x 0 y 0 -d(x 0 + y 0 ) (x 0 -y 0 )(r -d) , µ = -2d x 0 y 0 + r(x 0 + y 0 ) (x 0 -y 0 )(r -d) .
(3.12) y s = 0. Using the expression (3.2) of Q, the system (3.10) becomes:

ẋs =(r + d)x s - r 2 (1 + z s ) -dx 2 s , (3.13
x s = żs d(1 -z s ) , (3.15)
and differentiating (3.14) with respect to the time s gives:

ẋs = dz s (1 -z s ) + d ż2 s d 2 (1 -z s ) 2 .
Plugging this expression and (3.15) into (3.13) provides:

dz s (1 -z s ) + d ż2 s d 2 (1 -z s ) 2 = (r + d) żs d(1 -z s ) - r 2 (1 + z s ) -d ż2 s d 2 (1 -z s ) 2 ,
and, therefore,

zs (1 -z s ) + 2 ż2 s -(r + d)(1 -z s ) żs + rd 2 (1 + z s )(1 -z s ) 2 = 0. (3.16)
The system (3.15)-(3.16) can be solved explicitly. Let us define u s = (1 + z s )/(1 -z s ), so that

z s = u s -1 u s + 1 , ( 3
.17) from which we obtain:

żs = 2 us (1 + u s ) 2 , zs = 2ü s (1 + u s ) -4 u2 s (1 + u s ) 3 . (3.18)
Using these expressions in (3.16) provides:

2ü s (1 + u s ) -4 u2 s (1 + u s ) 3 2 u s + 1 + 2 4 u2 s (1 + u s ) 4 -(r + d) 2 u s + 1 2 us (1 + u s ) 2 + rd 2 2u s u s + 1 4 (u s + 1) 2 = 0 ⇔ 4ü s (1 + u s ) -4(r + d) us (1 + u s ) + 4rd u s (u s + 1) = 0 ⇔ üs -(r + d) us + rd u s = 0. (3.19)
Hence u satisfies a second-order ordinary differential equation, that solves in: The integration constants λ and µ can be expressed in terms of the initial conditions x 0 and z 0 = x 0 /y 0 :

u s = λe rs + µe d s , λ, µ ∈ R. ( 3
λ = 2d x 0 -d(1 + z 0 ) (1 -z 0 )(r -d) , µ = r -2dx 0 + z 0 r (1 -z 0 )(r -d) . (3.22)
This yields the announced result with y s = x s /z s .

Remark 4. Explicit expressions for x s and y s as functions of time and parameterized by r and d can be obtained. Using the expressions (3.12) of λ and µ in (3.21) and the relations between x s , y s and z s provides:

x s = - x 0 + y 0 -2x 0 y 0 x 0 + y 0 -2(d/r)x 0 y 0 exp(rs) + exp(ds) - d r x 0 + y 0 -2x 0 y 0 x 0 + y 0 -2(d/r)x 0 y 0 exp(rs) + exp(ds) + (1 -d/r)(y 0 -x 0 ) x 0 + y 0 -2(d/r)x 0 y 0 , (3.23)
y s = - x 0 + y 0 -2x 0 y 0 x 0 + y 0 -2(d/r)x 0 y 0 exp(rs) + exp(ds) - d r x 0 + y 0 -2x 0 y 0 x 0 + y 0 -2(d/r)x 0 y 0 exp(rs) + exp(ds) - (1 -d/r)(y 0 -x 0 ) x 0 + y 0 -2(d/r)x 0 y 0 . (3.24)
For the sequel, it is useful to notice that the constant

- λ r µ d = x 0 + y 0 -2x 0 y 0 x 0 + y 0 -2(d/r)x 0 y 0 (3.25)
which appears in (3.23) and (3.24) belongs to ]0, 1[ for all (x 0 , y 0 ) ∈]0, 1[ 2 ; this is due to the fact that d/r < 1.

Below, we list some properties of the solutions to the dynamical system (3.10). The trajectories of the solutions to (3.10) can be decomposed into four steps. In order to describe them, let us introduce

s 0 = 1 r -d log - µ d λ r = log x 0 + y 0 -2x 0 y 0 x 0 + y 0 -2(d/r)x 0 y 0 d -r (3.26)
and s ± as the only positive roots of the denominators of x s and y s in (3.11):

λe rs + µe d s ± 1. (3.27) Proposition 3.3. Let (x 0 , y 0 ) ∈]0, 1[ 2 .
(i) The stationary solutions to (3.10) are the saddle point (1, 1) and the attractive point (r/d, r/d).

(ii) s 0 ∈]0, ∞[ and s ± ∈]s 0 , ∞[.

(iii) inf{s + , s -} = s + (resp. s -) if and only if y 0 < x 0 (resp. y 0 > x 0 ).

(iv) lim s↑s + x s = -∞ and lim s↓s -y s = -∞.

The next points specify the behavior of the solution to (3.10).

(v) On [0, s 0 [, (x s , y s ) belongs to ]0, 1[ 2 and goes to (x s 0 , y s 0 ) = (0, 0) as s → s 0 .

(vi) On ]s 0 , inf{s + , s -}[, (x s , y s ) goes decreasingly to (x inf{s + ,s -} , y inf{s + ,s -} )-by "decreasingly" we mean that both coordinates decrease.

(vii) On ] inf{s + , s -}, sup{s + , s -}[, (x s , y s ) goes to (x sup{s + ,s -} , y sup{s + ,s -} ) decreasingly.

(viii) On ] sup{s + , s -}, ∞[, (x s , y s ) goes decreasingly to (r/d, r/d).

Proof of Item (i). To find the stationary solutions to (3.10), let us solve ẋ = 0 and ẏ = 0. With (3.10), we get Q(x, y) = 0 and Q(y, x) = 0. This directly implies that (x, y) = (1, 1) or (x, y) = (r/d, r/d), see Table 2. Let us now study the stability of these two equilibria. At (1, 1), the Jacobian of (3.10) is:

Jac(1, 1) = r 2 J -dI,
where J is the 2 × 2 matrix full of ones and where I is the 2 × 2 identity matrix. The eigenvalues of Jac(1, 1) are -d < 0 and r -d > 0, associated with the eigenvectors (1, -1) and (1, 1), respectively. By classical linearization methods (e.g. [START_REF] Verhulst | Nonlinear Differential Equations and Dynamical Systems[END_REF]Chap. 3]), we deduce that the point (1, 1) is a saddle point. At (r/d, r/d), the Jacobian of (3.10) is:

Jac(r/d, r/d) = d 2 J -rI.
The eigenvalues are -r < 0 and d -r < 0, associated with the eigenvectors (1, -1) and (1, 1), respectively. The point (r/d, r/d) is therefore attractive.

Proof of Item (ii) to (iv). Now we prove the different facts dealing with s 0 , s + and s -. First, (3.25) and the fact that r > d immediately imply that s 0 ∈]0, ∞[. Next, we show that (3.27) has on [0, ∞[ only one root, which belongs to ]s 0 , ∞[. For this, we shall start with proving that (3.27) is positive on [0, s 0 ]. Then, we shall show that (3.27) is decreasing in ]s 0 , ∞[ and goes to -∞ as s → ∞.

In order to prove the first point above, it is enough to show that (3.27) is positive at s = 0 and increasing on [0, s 0 [. (3.27) is positive at s = 0 simply because

λ + µ ± 1 = x 0 + y 0 y 0 -x 0 ± 1 = 2y 0 y 0 -x 0 or 2x 0 y 0 -x 0 = (1 -d/r)[x 0 + y 0 ± (y 0 -x 0 )] x 0 + y 0 -2(d/r)x 0 y 0 > 0.
To check that (3.27) is increasing on [0, s 0 [, we note that the derivative of (3.27) is positive on [0, s 0 [actually by construction of s 0 . Now we prove the second point. From (3.25) and since r > d, we obtain that (3.27) goes to -∞ as s → ∞. Also, by definition of s 0 , the derivative of (3.27) is negative on ]s 0 , ∞[, (3.27) is therefore decreasing on ]s 0 , ∞[. 14 x y q q q q q q q q q Figure 4: Solutions corresponding to several initial conditions. We see that depending on the initial condition, the solutions converge of diverge to infinity.

The fact that inf{s + , s -} equals s + (resp. s -) if and only if y 0 < x 0 (resp. y 0 > x 0 ) follows directly from (3.27).

Finally, since the numerators of x s and y s are negative on ]s 0 , ∞[, hence in particular at s ± , it is immediate that lim s↑s + x s = -∞ and lim s↓s -y s = -∞.

Proof of Item (v) to (viii). Let us first consider Item (v). By definition of s 0 , the numerators of x s and y s in (3.23) and (3.24) vanish for the first time at s 0 . Moreover, since s ± > s 0 , both denominators are non-zero at s 0 and x s 0 = y s 0 = 0. In particular, on [0, s 0 [, we have (x s , y s ) ∈ R 2 + . In fact, (x s , y s ) ∈]0, 1[ 2 . Indeed, on the segment {1}×]0, 1[, ẋs < 0 whereas on ]0, 1[×{1}, ẏs < 0: it is therefore not possible to go through these segments.

We turn to the proof of Item (vi). Thanks to (3.23) and (3.24), just after the time s 0 , (x s , y s ) belongs to the negative quadrant R 2 -. But for any (x, y) ∈ R 2 -, Q(x, y) 0 and Q(y, x) 0, see Table 2, in such a way that both x s and y s are decreasing as soon as they stay in this quarter plane, in other words for s ∈]s 0 , inf{s + , s -}[. At time inf{s + , s -}, one (or even the two if s + = s -, i.e. if x 0 = y 0 ) of x s and y s becomes infinite. In the sequel, let us assume that inf{s + , s -} = s + ; a similar reasoning would hold for the symmetrical case inf{s + , s -} = s -.

Let us show Item (vii). Just after s + , (x s , y s )

∈ (R + × R -) ∩ {(x, y) ∈ R 2 : Q(x, y) < 0, Q(y, x) < 0}.
The latter set is simply connected and bounded by the curve {(x, y) ∈ R 2 : Q(x, y) = 0}, see Table 2. Using classical arguments (see e.g. [START_REF] Verhulst | Nonlinear Differential Equations and Dynamical Systems[END_REF]), we obtain that it is not possible to go through this limiting curve on which ẋs = 0; this is why for any s ∈]s + , s -[, (x s , y s ) remains inside of this set.

We conclude with the proof of Item (viii). Just after the time s -, (x s , y s ) ∈ R 2 + ∩ {(x, y) ∈ R 2 : Q(x, y) < 0, Q(y, x) < 0}. For the same reasons as above, (x s , y s ) cannot leave this set and actually converges to (r/d, r/d).

Use of the characteristic curves to simplify the functional equation

Let us assume the existence of a solution P (x, y) to (3.1), and let us define g s = P (x s , y s ). Then: which looks like a first-order ODE for g, except that h depends on the boundary condition of P .

We first freeze the dependence on the solution in h, i.e. we solve the ODE (3.28) as if the term in the right-hand side were a known function. Using the solutions to the characteristic equations, we shall prove the following result: Proposition 3.4. Let h(x, y) be an analytical function on [0, 1[ 2 . Let (x 0 , y 0 ) ∈ R 2 . The solution to the ODE ġs + R(x s , y s )g s = h(x s , y s ), g 0 = P (x 0 , y 0 ),

where (x s , y s ) s 0 are the solutions (3.11) to the characteristic curve starting at (x 0 , y 0 ), is given by

g h s =P (x 0 , y 0 ) exp - s 0 R(x u , y u )du + s 0 h(x u , y u )e -s u R(xα,yα)dα du =:F (s, x 0 , y 0 , h). ( 3

.30)

Proof. Equation (3.29) is an inhomogeneous first-order ODE. The solution to the associated homogeneous equation is:

g s = R(x 0 , y 0 ) exp - s 0 R(x u , y u )du .
The announced result is deduced from the variation of constant method.

A solution P to (3.1) hence satisfies the following functional equation for all s, x 0 and y 0 :

P (x s , y s ) =P (x 0 , y 0 )e -s 0 R(xu,yu)du + s 0 h(x u , y u , P )e -s u R(xα,yα)dα du, (3.31)
with the function h defined in (3.3). Plugging the definitions (1.4) and (3.3) in (3.31), we obtain:

i,j 1 p i,j x i s y j s = P (x 0 , y 0 )e -s 0 R(xα,yα)dα - r 2 i 1 p i,1 s 0 i(x u ) i e -s u R(xα,yα)dα du - r 2 j 1 p 1,j s 0 j(y u ) j e -s u R(xα,yα)dα du +d s 0 x u y u 1 1 -x u + 1 1 -y u e -s u R(xα,yα)dα du. ( 3

.32)

Notice that the r.h.s. of (3.32) depends only on the p i,1 's and p 1,j 's, while the l.h.s. depends on all p i,j 's.

Computing the extinction probabilities by Monte-Carlo methods yields good results if we are interested in a given initial condition (i, j). We then have a complexity of order M × T . However, biologists may be interested in investigating the extinction probabilities when the initial condition (i, j) varies, and the method become computationally expensive.

Deterministic algorithm

For numerical approximations, we restrict ourselves to the computation of (p i,j ) i,j∈{1,...,N } for a positive (large) integer N . In this case, (1.7) can be approximated by the solution to a linear system.

Let us define p N = (p 1,1 , . . . , p 1,N , p 2,1 , . . . , p 2,N , . . . , p N,1 , . . . , p N,N ) T and T N is a N 2 × N 2 -matrix with five non-zero diagonals:

T N =          A 1 D 0 . . . 0 B 2,1 A 2 D . . . . . . 0 B 3,2 . . . . . . 0 . . . . . . . . . . . . D 0 . . . 0 B N,N -1 A N          where D = r 2(r+d) Id N , A i (i ∈ {1, . . . , N }) and B i,i-1 (i ∈ {2, . . . N }) are the N × N -matrices A i =           -1 r 2(r+d) 0 . . . 0 d r+d 2 i+2

Results

We start with r = 3 and d = 2. For the Monte-Carlo simulation, we use M = 200 and T = 5000. For the deterministic method, we use N = 50, so that (i, j) ∈ {1, . . . , 50} 2 . Estimators of the extinction probabilities p

i,j and p

i,j obtained respectively from the methods of Sections 4.1 and 4.2 are plotted in Figure 5. The results given by both methods are very similar, as shown by the statistics of Table 1.

In Table 1, we compute the square difference between the two predictions ( p i,j . For the latter, we consider only the couples (i, j) where p (1) i,j and p (2) i,j do not vanish (else, the fraction is either not defined or either 1 whatever the value of p (1) i,j ).

Mean

St.dev Min Max Square error 3.24 10 -5 2.47 10 -4 1.68 10 -36 4.63 10 -3 Absolute error 1.34 10 -3 5.53 10 -3 1.30 10 -18 6.81 10 -2 Relative error 4.49 10 -2 1.71 10 -1 9.80 10 -3 9.71 10 -1 Table 1: Square, absolute and relative differences between the predictions of the stochastic method when r = 3 and d = 2 (Section 4.1) and of the deterministic methods (Section 4.2). Recall that with M = 200, the width of the confidence interval (4.1) is 6.92 10 -2 .

To carry further the comparison of the stochastic and deterministic method, and to observe the influence of N on the quality of the approximation, we compute the relative quadratic error

Figure 1 :

 1 Figure 1: (a) Transition rates for the continuous-time pure-jump Markov process (X t , Y t ) t∈R + . (b) Transition probabilities of the embedded random walk, that we denote by (X t , Y t ) t∈N (here and throughout, N is the set {0, 1, 2, . . .}), with an abuse of notation.

  we have almost sure extinction of the population.

  x s = λre rs + µd e ds d(λe rs + µe d s + 1) , y s = λre rs + µd e ds d(λe rs + µe d s -1) , (3.11)

Figure 2 :Figure 3 :

 23 Figure 2: Differential system (3.10): Vector field (Q(x, y), Q(y, x)). The plain lines correspond to the sets {Q(x, y) = 0} or {Q(y, x) = 0}.Proof. First of all, since Q(x, y) and Q(y, x) are locally Lipschitz-continuous for x, y = 0, there is local existence and uniqueness. Let us introduce the new variable z s = x s /y s , which is defined as long as

.

  Let us also define the vector bN = (b 1N , . . . , b N N ) T ∈ R N ×N such that: b 1N = -

, 2 ( 1 r 2 ( 1 r 2 (

 21212 for i ∈ {2, . . . , N -1}, b N N = -r+d) p N +1,r+d) p N +1,2. . . r 2(r+d) p N +1,N -r+d) p N +1,N + p N,N +1 p i,N +1 and p N +1,N are approximations of p i,N +1 given by Proposition 2.3. With these notations, (1.7) rewrites as T N p N = b N .

  | and the relative difference | p

Figure 6 :

 6 Figure 6: Estimation of the extinction probabilities p i,j 's when r = 2.002 and d = 2: (a) with the Monte-Carlo method of Section 4.1. (b) with the deterministic method of Section 4.2.

Figure 7 :

 7 Figure 7: Estimation of the extinction probabilities p i,j 's when r = 2.002 and d = 2: Decrease with N of the log of the relative quadratic errors (4.2) between the deterministic and Monte-Carlo methods (M = 200).

  Let the probabilities (p i,1 ) i 1 be given. Then the probabilities (p i,j ) i,j 1 are completely determined.

	2 Existence of a solution
	2.1 Dirichlet problem
	We first establish that the extinction probabilities p i,j 's (1.3) solve the Dirichlet problem (1.6)-(1.7).
	Proposition 2.1. (i) The extinction probabilities (p i,j ) i,j 1 are solutions to the Dirichlet problem
	(1.7) with boundary condition (1.6). Uniqueness of the solution may not hold, but the extinction
	probabilities (p i,j ) i,j∈N define the smallest positive solution to this problem.
	(ii) Proof. We begin with Point (i). Equation (1.7) is obtained by using the strong Markov property at
	the time of the first event. Let us denote by K the transition kernel of the discrete-time Markov chain
	(X t , Y t ) t∈N
	This provides a new formulation of the
	solution of (1.7), that is however uneasy to work with numerically. Hence, in Section 4, we propose
	two different approaches leading to numerical approximations of the solution of the Dirichlet problem
	(1.6)-(1.7), that are based on stochastic and deterministic approaches.
	In conclusion, we provide here several approaches to handle the extinction probabilities of the in-
	homogeneous random walk (X t , Y t ) t∈R + of our problem. Estimates from [4] are recalled and the
	recurrence equation (1.7) is solved numerically and theoretically, pending further investigation of the
	PDE formulation.

* ; we have:

Acknowledgements

The authors warmly thank the referees for useful comments and suggestions.

Proposition 3.5. Let s 0 > 0 be defined in (3.26). We have: Before proving Proposition 3.5, let us show that the different quantities that appear in its statement are well defined-indeed, this is a priori not clear: as α → s 0 , x α → 0 and y α → 0, in such a way that R(x α , y α ) → ∞, see (3.2). Lemma 3.6. Let i, j ∈ N. Then lim u→s 0 (x u ) i (y u ) j e u 0 R(xα,yα)dα is finite if and only if i + j 1-and equals zero if and only if i + j 2.

Proof. First, since the only zero of any function of the form α exp(as) + β exp(bs) with αβ < 0 and a = b has order one, the following function has a simple zero at s 0 :

Thanks to this and since s + , s -> s 0 , both x s and y s have a zero of order 1 at s 0 . Moreover, with λ and µ defined in (3.12), we obtain that:

It is indeed easy to check that the derivative of the logarithm of (3.35) is equal to R(x u , y u ), for which we have an explicit expression, see (3.2), (3.23) and (3.24). From (3.35), we see that e u 0 R(xα,yα)dα has three poles, namely at s 0 , s + , s -. The zero at s 0 has order one by using again the considerations on the zeros of (3.34). In particular, Lemma 3.6 follows immediately.

Proof of Proposition 3.5. Start by multiplying (3.32) by e s 0 R(xα,yα)dα and then let s → s 0 . Since P (x, y) = xy i,j 1 p i,j x i-1 y j-1 , see (1.4), and since lim s→s 0 x s y s e s 0 R(xα,yα)dα = 0, see Lemma 3.6, we obtain that lim

which concludes the proof of Proposition 3.5.

Remark 5. When (x 0 , y 0 ) ∈ (0, 1) 2 , we also have (x u , y u ) ∈ (0, 1) 2 for all u ∈ (0, s 0 ). Thus it is possible to plug approximations of the p 1,j 's and p i,1 's into (3.33) thanks to the terms (x u ) i and (y u ) j . Using that p 1,i = 2d/(ri) + o(1/i) when i → +∞, it is possible to find I 0 sufficiently large so that:

Thus:

The latter expression shows that numerically, one can restrict to the computation of a finite number of probabilities p i,1 , for i I 0 .

Numerical results

In this section, we present two different ways of approximating the extinction probabilities p i,j .

Probabilistic algorithm

A first possibility, if we are interested in a given initial condition (i, j), is to approximate p i,j by Monte-Carlo simulations. For T > 0 large, we simulate M paths (X t , Y t ) t∈{1,...,T } started at (i, j), for ∈ {1, . . . , M }, independent and distributed as the process (X t , Y t ) t∈{1,...,T } . The extinction probability is estimated by:

The estimator p M,T is the proportion of paths that have gone extinct before time T .

Proposition 4.1. Let (i, j) be the initial condition. The estimator p M,T has the following properties:

(i) It is a convergent and unbiased estimator of

, and hence we have the following asymptotic 95% confidence interval for p i,j :

Proof. These results are straightforward consequences of the law of large numbers and central limit theorem, given that 1l {∃t T,X t Y t =0} are independent Bernoulli random variables with parameter P i,j [τ 0 T ].

i,j is the deterministic approximation for N ∈ {10, . . . , 50} and p

i,j is either given by the deterministic approximation with N = 50, or by the stochastic approximation p

i,j , the decrease in the quadratic errors stops around N = 18 around 0.0891. This corresponds roughly to the stochastic error of the law of large numbers (4.1) which depends only on M . In the second case, when p

i,j is the deterministic approximation with N = 50, the relative quadratic errors decrease exponentially fast in exp(-0.6842 N ) (R 2 = 99.92%).

In a second experiment, we choose r = 2.002 and d = 2. This case is more interesting in population ecology, since small populations are of interest when they are fragile and endangered species. For the Monte-Carlo simulation, we use M = 200 and T = 5000. For the deterministic method, we use N = 100, so that (i, j) ∈ {1, . . . , 100} 2 . The estimated extinction probabilities p (1) i,j and p (2) i,j are plotted in Figure 6, and statistics are computed in Table 2. Again, results from both methods are similar. This is confirmed by computing the relative quadratic errors, with the p (2) i,j 's obtained from the deterministic method and the p

i,j 's from the Monte-Carlo method. The decrease of this error is exponential with N in exp(-0.0619 N ) (R 2 = 98.98%) as shown in Figure 6(c). It can be noticed that in this case, the per! formances of the Monte-Carlo method match better the one of the deterministic algorithm. This is due to the fact that Monte-Carlo methods fail to produce good estimates of small probabilities (see [START_REF] Clémençon | On computer-intensive simulation and estimation methods for rare event analysis in epidemic models[END_REF] and references therein). When the probabilities p (3) i,j 's are given by the deterministic method with N = 50 in (4.2), we have as in the previous case (r = 3) an exponential decrease of the relative quadratic error in exp(-0.1092 N ) (R 2 = 95.07%).

Mean

St.dev Min Max Square error 6.27 10 -3 6.25 10 -3 8.42 10 -10 5.31 10 -2 Absolute error 6.89 10 -2 3.90 10 -2 2.90 10 -5 2.31 10 -1 Relative error 2.22 10 -1 2.17 10 -1 1.02 10 - 4 1

Table 2: Square, absolute and relative differences between the predictions of the stochastic method when r = 2.002 and d = 2 (Section 4.1) and of the deterministic methods (Section 4.2). As in Table 1, since M = 200, the width of the confidence interval (4.1) is 6.92 10 -2 .


Emmanuel Soubiès

Pierre Weiss
email: pierre.armand.weiss@gmail.com

Xavier Descombes
email: xavier.descombes@inria.fr

Emmanuel Soubies
email: esoubies@gmail.com

A 3D SEGMENTATION ALGORITHM FOR ELLIPSOIDAL SHAPES. Application to nuclei extraction

Keywords: Nuclei segmentation, 2D and 3D images, graph-cuts, marked point processes, ellipses and ellipsoids, multiple object detection, multiple birth and cut, bio-imaging

A 3D segmentation algorithm for ellipsoidal shapes.

Application to nuclei extraction.

INTRODUCTION

Cell or nuclei segmentation in 2D and 3D is a major challenge in bio-medical imaging. New microscopes provide images at higher resolutions, deeper into biological tissues, leading to an increasing need for automatic cell delineation. This task may be easy in certain imaging modalities where images are well resolved and contrasted, but it remains mostly unresolved in emerging fluorescent microscopes dedicated to live imaging such as confocal, bi-photon, or selective plane illumination microscopes. These modalities suffer from multiple degradations such as light attenuation in the sample, heavy noise and spatially varying blur that make the segmentation task hard even for human experts Our aim in this work is to propose a segmentation algorithm robust to such situations. Since images are heavily deteriorated, standard methods aiming at finding contours based on a sole regularity assumption such as active contours or Mumford-Shah derivatives fail for the segmentation. This observation led us to introduce strong shape priors: cells are modelled as ellipses or ellipsoids that should fit the image contents. Unfortunately, adding geometrical constraints makes the optimization problems highly non convex and appeal for the development of new global opti- mization methods.

Following recent works [START_REF] Descombes | Object extraction using a stochastic birth-and-death dynamics in continuum[END_REF][START_REF] Descombes | Stochastic geometry for image analysis[END_REF][START_REF] Gamal Eldin | Multiple birth and cut algorithm for multiple object detection[END_REF], we use randomized algorithms that allow to escape from local minima. These algorithms are based on marked point processes. The Marked Point Process (MPP) approach [START_REF] Baddeley | Stochastic geometry models in high-level vision[END_REF][START_REF] Dong | Detection of rolling leukocytes by marked point processes[END_REF] consists in estimating a configuration of geometric objects (in our case ellipses or ellipsoids) whose number, location and shape are un-known. It has proved to be very efficient in numerous image analysis applications as it allows the combination of radiometric information with strong geometrics constraints on the objects but also at a global scale. Defined by a density against the Poisson process measure, its main advantage is to consider a random number of objects and can be considered as an extension of the Markov Random Field approach. A review of this approach and its applications can be found in [START_REF] Descombes | Stochastic geometry for image analysis[END_REF].

The objects are defined on a state space χ = I × M by their location and their marks (i.e. geometric attributes). The associated marked point process X is a random variable whose realisations are random configurations of objects. Considering a Gibbs process, the modeling consists of an energy construction. Similarly to the Bayesian framework, this energy can be written as the sum of a data term and a prior. In this paper we consider a pairwise interactions prior that forbids intersections between objects. Once the model defined, the solution is obtained by minimizing the energy. This energy being highly non-convex requires stochastic dynamics, such as MCMC methods, to be minimized. The Reversible Jump MCMC embeded in a simulated annealing framework is a natural candidate for this task [START_REF] Green | Reversible jump markov chain monte carlo computation and bayesian model determination[END_REF]. However, in case of simple constraints such as non overlap, the recently proposed multiple birth and death algorithm is preferable [START_REF] Descombes | Object extraction using a stochastic birth-and-death dynamics in continuum[END_REF]. To avoid the fastidious calibration of annealing parameters, we propose to revisit the combination of the multipe births principle with the graph cut paradigm proposed by [START_REF] Gamal Eldin | Multiple birth and cut algorithm for multiple object detection[END_REF].

The paper is organized as follows. We formalize the segmentation problem as a minimization problem in section 2. Section 3 begins by a global algorithm description and is followed by a precise description of each algorithm step. We finish by presenting numerical results in section 4.

PROBLEM STATEMENT

Figure [1] contains typical examples of images encountered in biology. It is readily seen from these images that most nuclei contours can be well approximated by ellipses or ellipsoids, at least at a coarse scale. Moreover these nuclei cannot overlap due to obvious physical considerations. We thus formulate our segmentation problem as that of finding a set of non overlapping ellipsoids that fit the image contents. We formalize this statement in the latter.

Let C n , n ∈ N denote the set of configurations containing n objects that do not overlap. An element

x ∈ C n is a set of n non overlapping objects. Since the number of nuclei in the configuration is unknown, we aim both at finding this number n * and the best configuration x ∈ C n * with respect to a certain data fi- delity term f (x). Our optimization problem can thus be formulated as follows. Let

g(n) = min x∈C n f (x)
denote the minimum value of f in the set C n . We wish to find both n * = arg min n∈N g(n)

and

x * = arg min

x∈C n * f (x).
By convention, we assume that C 0 = / 0 and that min x∈C 0 f (x) = 0. The data term f should thus be negative for configurations that are likely to represent the nuclei parameters and positive otherwise. We detail how the ellipses are parametrized and the construction of such a function in the following paragraphs.

Object modelling In 2 dimensions, ellipses are parameterized using 5 parameters (see Figure 2):

• (x, y) ∈ Ω: center coordinates which should belong to the image domain Ω.

• θ ∈ [0, 2π[: angle with the horizontal direction.

• 0 < λ -< b < a < λ + : describe the ellipses minor and major axes size. λ -and λ + are user defined parameters that describe the nuclei maximal size and ellipticy. In 3 dimensions, nuclei are parameterized using 9 parameters:

• (x, y, z) ∈ Ω: center coordinates.

• φ, θ, ψ ∈ [0, 2π[ 3 : Euler angles to define the ellipsoids orientations.

• 0 < λ -< c < b < a < λ + : axes lengths.

Overall, it can be seen that objects belong to a state space χ defined as a parallelepiped:

χ = Ω × [0, 2π[×[λ -, λ + ] 2 (1)
in 2D and

χ = Ω × [0, 2π[ 3 ×[λ -, λ + ] 3 (2) in 3D.
In this paper, the objects are denoted ω and their boundary is denoted ∂ω.

Data term Let u : Ω → R denote a grayscale image. In order to define the data term f (x), we associate an elementary energy U d (ω, u) to each element ω ∈ x and set:

f (x) = ∑ ω∈x U d (ω, u).
(3)

The function U d (ω, u) ∈ [-1, 1] should be negative if the object ω is well positioned on the image and positive otherwise.

In fluorescence microscopy, nuclei are usually characterized by bright region surrounded by a dark background since they are stained or genetically modified in order to express a fluorescent protein. Unfortunately, their radiometry is not constant due to local bleaching or light attenuation in the deepest layers. We thus need to construct an energy that is contrast invariant, meaning that local modifications of the radiometry shall not affect the energy. Such an energy can be constructed easily by considering the normal to the image level lines ∇u |∇u| where ∇u denotes the usual gradient in R d and |∇u| denotes the gradient magnitude in the standard Euclidean norm. This tool is well known to be contrast invariant. Let us define an energy U for a given object ω as:

U(ω) = 1 |∂ω| ∂ω ∇u(x) |∇u(x)| 2 + ε 2 , n(x) dx (4)
where •, • denotes the standard scalar product, |∂ω| denotes the length of the object boundary, n(x) denotes the outward normal to ω at location x ∈ ∂ω and ε is a regularization parameter that discard faint transitions. The behavior of this energy is illustrated on Fig. 3. Overall, it does what is expected, but as can be seen on the illustration b) and d) in Fig. 3, badly located ellipses might have a negative energy and be kept in the final configuration. It is thus necessary to modify U in order to promote well located objects only. A simple way to do so consists in setting:

U d (ω, u) = ψ(U(ω), s)
where s ∈] -1, 0] is an acceptance threshold for the objects. and

ψ(α, s) = min( 1 s + 1 α - s s + 1
, 1). This function is illustrated on Figure 4 Other data terms based on the contrast between the object interior and the background as presented by [START_REF] Gamal Eldin | Multiple birth and cut algorithm for multiple object detection[END_REF] (in dimension 2) could also be used but present two drawbacks: first they require to compute an integral over the interior of the domain while the proposed approach consist in computing a boundary integral which is faster. Second, such measures might be inaccurate in the case of very dense media, where the background can be difficult to extract. Finally our measure is contrast invariant, which is central for the targeted applications.

MULTIPLE BIRTH AND CUT ALGORITHM (MBC)

The Multiple Birth and Cut algorithm (MBC) has been proposed by [START_REF] Gamal Eldin | Multiple birth and cut algorithm for multiple object detection[END_REF] for counting flamingos in a colony. In this section, we describe the different steps of the MBC algorithm (Algorithm 1).

The main idea consists in generating two random configurations of non-overlapping objects x and x ′ (birth step) and then keep the subset of objects in x ∪ x ′ that minimizes f (cut step). This process is iterated and decreases f at each iteration. The cut step can be performed efficiently using a Graph Cut algorithm [START_REF] Boykov | Fast approximate energy minimization via graph cuts[END_REF][START_REF] Kolmogorov | What energy functions can be minimized via graph cuts[END_REF]. We describe this algorithm more formally below:

Interestingly, this algorithm contains only one parameter N (the number of objects generated in a con- Generate of a new configuration x ′ using Algorithm 2.

5:

x

[n+1] ←-Cut(x n ∪ x ′ ) 6:
n ←n + 1 7: end while figuration). We observed experimentally that this parameter might affect slightly the speed of convergence but not the segmentation accuracy. This algorithm is thus much easier to tune than more standard RJM-CMC based dynamics.

Birth step

A new configuration x ′ of non-overlapping objects is generated. Note that only objects which are in the same configuration have to respect the nonoverlapping constraint, but two objects in different configurations can intersect as can be seen on Figure 5. The birth step is detailed in Algorithm 2. The fourth step of this algorithm can be efficiently implemented using a lookup table and the fast intersection algorithm proposed in the latter.

Cut step

This step consists in selecting the best configuration of non-overlapping objects in (x [n] ∪ x ′ ). To perform this optimization, a weighted graph is constructed. The nodes of this graph are the objects ω of the two configurations x [n] and x ′ . This graph also possesses Algorithm 2 Birth step Require: N, n max .

1: Set k = 0, n = 0, x ′ = / 0. 2: while k < N and n < n max do 3:

Construct an object ω ′ by generating a random vector uniformly in χ.

4:

If ω ′ intersects an object in x ′ , set n = n + 1 and go back to 3.

5:

Otherwise set x ′ = x ′ ∪ {ω ′ }, k = k + 1, n = 0 and go back to 3. 6: end while two special nodes, the source 's' and the sink 't'. The weights should belong to [0, 1] ∪ {+∞} and are defined using the data term U d (ω, u) by:

W (ω) = (1 -U d (ω, u))/2.

Graph construction

Each object of the configuration (x [n] ∪ x ′ ) is linked to the source and the sink. The difference between the objects ω i ∈ x [n] and the objects ω j ∈ x ′ is that the objects ω i ∈ x [n] are linked to the source with a weight equal to the data energy W (ω) and to the sink with a weight equal to 1 -W (ω), while it is the reverse for the objects ω j ∈ x ′ .

The weights associated to edges linking two objects are non zero only when two objects intersect. If ω 1 ∈ x [n] (current configuration) intersects with ω 2 ∈ x ′ (new configuration), the link from ω 1 to ω 2 is set to ∞ and the link from ω 2 to ω 1 is set to zero1 . This ensures that the cut step generates an admissible configuration (with no overlapping objects). Figure 6 summarises the graph construction of the configurations on Figure 5. The nuclei to detect are represented by black ellipses.

Cut

Once the graph is constructed, we perform a cut that consists in partitioning the vertices into two disjoint subsets. One contains the source and the other the sink. The cut realized is the one with minimal cost (the one minimizing the sum of the weights of the removed edges).

After the cut step, if ω i ∈ x [n] is in the sub-graph containing the source, we keep it, otherwise we remove it. On the contrary the objects ω j ∈ x ′ are Figure 6: Graph corresponding to the figure 5 kept only if they belong to the sub-graph that contains the sink. This difference of interpretation between the two configurations combined with the different weights to the source and the sink, ensure that in case where an object of x [n] and an object of x ′ intersect, only one can be kept.

The cut step is implemented using the graph-cut code developed by Yuri Boykov and Vladimir Kolmogorov in [START_REF] Boykov | Fast approximate energy minimization via graph cuts[END_REF][START_REF] Boykov | An experimental comparison of min-cut/max-flow algorithms for energy minimization in vision[END_REF][START_REF] Kolmogorov | What energy functions can be minimized via graph cuts[END_REF].

A FAST DETERMINATION OF ELLIPSES INTERSECTION

One of the proposed algorithm bottleneck is the fast determination of whether two ellipsoids intersect or not. In this section, we present a fast algorithm to answer that question and prove theoretically that only a few arithmetic operations suffice to provide the answer with a low error rate.

Let ω be an ellipse or an ellipsoid. It can be defined using a quadratic function

Q ω as ω = {x ∈ R d , Q ω (x) = 1}.
The quadratic function Q can be defined by:

Q ω (x) = A(x -c), (x -c) (5)
where c denotes the object center and A is positive definite matrix defined by:

A = P -1 DP = P T DP.
where P is a rotation matrix and D is a positive diagonal matrix. In 2D, P is defined by:

P = cos(θ) sin(θ) -sin(θ) cos(θ)
and

D = 1 a 2 0 0 1 b 2 .
In 3D, the notation become cumbersome and we leave them to the reader.

Let ω 1 and ω 2 be two ellipses or ellipsoids. In order to know whether they intersect or not, we can find the minimal level set of Q ω 2 which intersects the boundary of ω 1 . If this level set is associated to a value greater than 1, the ellipses are separated, otherwise they overlap. This idea can be formulated as the following minimization problem:

min x∈R d ,Q ω 1 (x)≤1 Q ω 2 (x) (6)
This problem consists of minimizing a quadratic function over convex set. Projected descent methods can thus be used. Unfortunately, there exists no closed form solution to the problem of projection of a point on an ellipse. We thus need to simplify the constraint set:

min

Q ω 1 (x)≤1 Q ω 2 (x) = min A 1 (x-c 1 ),(x-c 1 ) ≤1 A 2 (x -c 2 ), (x -c 2 ) = min √ A 1 (x-c 1 ), √ A 1 (x-c 1 ) ≤1 A 2 (x -c 2 ), (x -c 2 ) min y- √ A 1 c 1 2 2 ≤1 A 2 (A -1 2 1 y -c 2 ), (A -1 2 1 y -c 2 ) .
where

y = √ A 1 x. In this reformulation, the constraint set Y = {y ∈ R d , y- √ A 1 c 1 2 2 ≤ 1} is a simple l 2 -ball and the function F(y) = A 2 (A -1 2 1 y -c 2 ), (A -1 2 1 y - c 2
) is a strongly convex differentiable function. We can thus use a projected gradient descent that writes:

Algorithm 3 Detection of overlapping ellipsoids

Require: Q ω 1 , Q ω 2 , ε > 0. 1: Set k = 0, y 0 = c 1 +c 2 2 . 2: Set µ = b 2 1 a 2 2 , L = a 2 1 b 2 2 . 3: Set τ = 2 µ+L . 4: while y k+1 -y k ≥ ε do 5: y k+ 1 2 = y k -τ∇F(y k ). 6: y k+1 = Π Y y k+ 1 2 . 7: k = k + 1.
8: end while 9: If F(y k ) >= 1 return 0 (the ellipsoids do not intersect with high probability). 10: If F(y k ) < 1 return 1 (the ellipsoids intersect).

Let y * denote the solution of the above problem. The previous algorithm comes with the following guarantees: Theorem 1. After k iterations, y k satisfies:

F(y k ) -F(y * ) ≤ µ 2 y 0 -y * 2 2 Q F -1 Q F + 1 2k y k -y * 2 2 ≤ y 0 -y * 2 2 Q F -1 Q F + 1 2k where Q F = a 2 1 b 2 2 a 2 2 b 2 1 ≤ λ 4 + λ 4 - . Proof The Hessian of F is H F (y) = 2A -1 2 1 A 2 A -1 2 1 .
Since A 1 and A 2 are products of orthogonal and diagonal matrices (A = P T DP), the eigenvalues of H F (y) can be easily bounded:

λ min [H F (y)] ≥ b 2 1 a 2 2 λ max [H F (y)] ≤ a 2 1 b 2 2 The function F is thus µ-strongly convex with µ ≥ b 2 1 a 2 2 and its gradient is L-Lipschitz with L ≤ a 2 1 b 2 2
. Using standard convergence theorems in convex analysis [START_REF] Bertsekas | Nonlinear programming[END_REF], we obtain the announced result.

The conditioning number Q F depends solely on the ratio between the major axis and the minor axis sizes and not on the dimension d. This algorithm will thus be as efficient in 3D as in 2D. For two circles the ratio Q F is equal to a b = 1 and the algorithm provides the exact answer after one iteration. For elliptic ratios of 2, Q f = 16 and in the worst case, after 18 iterations, the algorithm returns a point y k that is 100 times closer to the intersection than y 0 . We also tested an accelerated algorithm by [START_REF] Nesterov | Introductory lectures on convex optimization: A basic course[END_REF], where the convergence rate is of order

√ Q F -1 √ Q F +1 2k
but it did not improve the computing times.

In our problems the ratio between a and b is always less than 2 and the algorithm usually converges in just a few iterations (2 to 10 depending on the problem).

Acceleration by local perturbations

When the objects variability is important, the state space size increases and affects the convergence speed of the MBC algorithm. This problem is particulary important in 3D since ellipsoids are defined by 9 parameters instead of 5 for the 2 dimensional case.

In order to improve the convergence speed, [START_REF] Gamal-Eldin | A fast multiple birth and cut algorithm using belief propagation[END_REF] proposed to insert a selection phase in the birth step. This selection phase consists in generating a dense configuration of objects at similar locations and to keep the best ones using Belief Propagation in order to form the new configuration.

In this paper, we propose another heuristic in order to increase the convergence speed. We propose to alternate between two different kinds of birth steps. The first one is that proposed in algorithm 2. The second one consists in perturbating locally the current configuration. This principle mimics the proposition kernels used in RJMCMC algorithms [START_REF] Perrin | A marked point process model for tree crown extraction in plantations[END_REF]. The idea behind this modification is that after a while, most objects are close to their real location and that local perturbations may allow much faster convergence than fully randomized generation. This algorithm is described in details in Algorithm 4. Construct an object ω ′ by local perturbation of ω k ∈ x[n].

3:

If ω ′ intersects an object in x ′ , set k = k + 1 and go back to 2.

4:

Otherwise set x ′ = x ′ ∪ {ω ′ }, k = k + 1 and go back to 2. 5: end while

Local perturbations

A given object ω in x [n] is described by a set of parameters λ ∈ χ (see equations 1 and 2). We generate the new object ω ′ by setting its parameters λ ′ = λ + z where z is the realization of a random vector Z dis-tributed uniformly in χ ε where :

χ ε = [-δ xy , δ xy ] 2 × [-δ ab , δ ab ] 2 × [0, 2π[ in 2D and χ ε = [-δ xyz , δ xyz ] 3 × [-δ abc , δ abc ] 3 × [0, 2π[ 3 in 3D.
The value of the different δ describes the perturbation extent. We observed that small values accelerates the convergence speed.

Comparison of the convergence speed

We have tested this method in order to compare the speed of convergence of the MBC algorithm and the MBC algorithm with local perturbation. Figure 7 presents the energy evolution with respect to time for both MBC and MBC with local perturbations (denoted MBC with LP) on the same image (the 3D nuclei of Drosophila embryo). The segmentation result is presented on Figure 14 (image size 700 × 350 × 100). These results show that the MBC with LP algorithm strongly improve the MBC algorithm.

RESULTS

In this section, we present some practical results in 2D and 3D. Figure 9 shows the segmentation result on a Drosophila embryo obtained using SPIM imaging. This is a rather easy case, since nuclei shapes vary little. The images are impaired by various defects: blur, stripes and attenuation. Despite this relatively poor image quality, the segmentation results are almost perfect. The computing time is 5 minutes using a c++ implementation. The image size is 700 × 350. Figure 10 presents a more difficult case, where the image is highly deteriorated. Nuclei cannot be identified in the image center. Moreover, nuclei variability is important meaning that the state space size χ is large. Some nuclei are in mitosis (see e.g. top-left). In spite of these difficulties, the MBC algorithm provides acceptable results. They would allow to make statistics on the cell location and orientation, which is a major problem in biology. The computing times for this example is 30 minutes. Nuclei segmentation is a major open problem with a large number of other applications. In Figure 11, we attempt to detect the spermatozoon heads. The foreseen application is tracking in order to understand their collective motion. Figure 12 presents a multicellular spheroid surrounded by a circle of pillars. A typical problem consists in determining whether the pillars have moved afetr a while, due to the tumor dynamic. They thus need to detect the pillars and evaluate their motion. Figure 11 and Figure 12 show that the segmentation results are extremely precise. The spheroid segmentation presented in Figure 16 is less precise than the previous ones due to an important cell variability and to the fact that the images are extremly blurry in the Z direction. For that case, image restoration algorithms and the design of new energies robust to strong perturbations seem important.

Figure 1 :

 1 Figure 1: Example of a SPIM image (Multicellular tumor spheroid).

Figure 2 :

 2 Figure 2: Parameters of the ellipse.

Figure 3 Figure 4 :

 34 Figure 3

Figure 5 :

 5 Figure 5: Two configuration on an image (the black ellipses are the object to detect)

Figure 7 :

 7 Figure 7: Comparison of the MBC and MBC with LP algorithms

Figure 8 :

 8 Figure 8: 2D nuclei of Drosophila embryo

Figure 10 :

 10 Figure 10: 2D segmentation of a multicellular tumor spheroid (Fig 1)

Figure 11 :

 11 Figure 11: Segmentation of a spermatozoon colony (5 minutes). Image size: 2000 x 1024.

Figure 12 :

 12 Figure 12: Nano pillars detection (less than 1 minute). Image size: 840 x 800.

Figure 13 :

 13 Figure 13: 3D Drosophila embryo nuclei

Figure 15 :

 15 Figure 15: 3D multicellular tumor spheroid

  Algorithm 1 Multiple Birth & Cut algorithm

	Require: N
	1: Generate a configuration x [0] with Algorithm 2
	2: n ←-0 3: while (Not converged) do
	4:

  Algorithm 4 MBC algorithm with local perturbations

	Require: N
	1: Generate a configuration x [0] using Algorithm 2.
	2: n ←-0 3: while (Not converged) do
	4:	Generate a uniformly distributed random num-
	5: 6:	ber r ∈ [0, 1]. if r < p then Generate a new configuration x ′
		using Algorithm 2.
	7:	else
	8:	Generate a new configuration x ′
		using Algorithm 5.
	9:	end if
	10: 11: 12: end while x [n+1] ←-Cut(x n ∪ x ′ ) n ←-n + 1
	Algorithm 5 Birth step with local perturbation
	Require: x [n] .
	1: while k < size(x [n] ) do
	2:

When two objects intersect the link affected by a weight of ∞ is always the link from the object of the current configuration to the object of the new configuration.

ACKNOWLEDGEMENTS

This work was partially funded by the Mission pour l'interdisciplinarité from CNRS, Région Midi Pyrénées, PEPII CASPA3D and ANR SPHIM3D.

The authors wish to thank F. Malgouyres and J. Fehrenbach for interesting discussions. They also thank Valérie Lobjois, Charlotte Emery, Jacques Figure 16: 3D segmentation of the multicellular tumor spheroid (Fig 15)

Thomazeau, Paul Escande and Bernard Ducommun for their help in this project. They thank all the ITAV staff for their warm welcome in a biology laboratory.