Rodolphe Garbit
email: rodolphe.garbit@univ-angers.fr

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

ON THE EXIT TIME FROM A CONE FOR RANDOM WALKS WITH DRIFT

Keywords: Mathematics Subject Classification. 60G40, 60G50, 05A16 Random walk, Cones, Exit time, Laplace transform

We compute the exponential decay of the probability that a given multidimensional random walk stays in a convex cone up to time n, as n goes to infinity. We show that the latter equals the minimum, on the dual cone, of the Laplace transform of the random walk increments. As an example, our results find applications in the counting of walks in orthants, a classical domain in enumerative combinatorics.

Introduction and main results

1.1. General context. For general random processes in R d , d 1 (including in particular Brownian motion and random walks), it is at once important and natural to study the first exit times τ K from certain domains K. Precisely, for discrete-time random processes (S n ) n 0 , τ K is defined by [START_REF] Bañuelos | Brownian motion in cones[END_REF] τ K := inf{n 1 : S n / ∈ K}.

Indeed, these random times carry much valuable information on the process. As an example, the fruitful theory of random walks fluctuations (see, e.g., Spitzer [START_REF] Spitzer | Principles of random walk[END_REF]) is based on the analysis of the τ K for compact domains K.

In a recent past (1990 to present), the case of cones K has arisen a great interest in the mathematical community, due to interactions with many areas: First, certain random walks in conical domains can be treated with representation theory [START_REF] Biane | Quantum random walk on the dual of SU(n)[END_REF][START_REF] Biane | Minuscule weights and random walks on lattices[END_REF] (in that case, the cones are Weyl chambers related to Lie algebras). Further, the exit times τ K are crucial to construct conditioned random walks in cones, which appear in the theory of quantum random walks [START_REF] Biane | Quantum random walk on the dual of SU(n)[END_REF][START_REF] Biane | Minuscule weights and random walks on lattices[END_REF], random matrices [START_REF] Dyson | A Brownian-motion model for the eigenvalues of a random matrix[END_REF], non-colliding random walks [START_REF] Denisov | Conditional limit theorems for ordered random walks[END_REF][START_REF] Eichelsbacher | Ordered random walks[END_REF], etc. In another direction, the probability (2)

P x [τ K > n]
admits a direct combinatorial interpretation in terms of the number of walks starting from x and staying in the cone K up to time n. These counting numbers are particularly important in enumerative combinatorics [START_REF] Bousquet-Mélou | Walks with small steps in the quarter plane[END_REF][START_REF] Fayolle | Some exact asymptotics in the counting of walks in the quarter plane[END_REF][START_REF] Johnson | Towards a combinatorial understanding of lattice path asymptotics[END_REF], and are the topic of many recent studies.

For processes with no drift, the exit times τ K from cones are now well studied in the literature. The case of Brownian motion was solved by DeBlassie [START_REF] Deblassie | Exit times from cones in R n of Brownian motion[END_REF] (see also Bañuelos and Smits [START_REF] Bañuelos | Brownian motion in cones[END_REF]): He showed that the probability (2) satisfies a certain partial differential equation (the heat equation), and he solved it in terms of hypergeometric functions. Concerning discrete-time random processes, in the one-dimensional case, the asymptotic behavior of the non-exit probability (2) is well known, as well as that of (3) P x [S n = y, τ K > n]

(called a local limit theorem), thanks to the theory of fluctuations of random walks [START_REF] Spitzer | Principles of random walk[END_REF]. In higher dimension, some sporadic cases have first been analyzed: We may cite [START_REF] Doumerc | Exit problems associated with finite reflection groups[END_REF], for which there exists a strong underlying algebraic structure (certain reflexion groups are finite), or the case of Weyl chambers, which has been considered in [START_REF] Denisov | Conditional limit theorems for ordered random walks[END_REF][START_REF] Eichelsbacher | Ordered random walks[END_REF]. For more general cones, but essentially for random walks with increments having a finite support, Varopoulos [START_REF] Varopoulos | Potential theory in conical domains[END_REF] gave lower and upper bounds for the probability [START_REF] Biane | Quantum random walk on the dual of SU(n)[END_REF]. The first author of the present article showed in [START_REF] Garbit | Temps de sortie d'un cône pour une marche aléatoire centrée[END_REF] that for general random walks, the probability (2) does not decay exponentially fast. More recently, Denisov and Wachtel [START_REF] Denisov | Random walks in cones[END_REF] provided the exact asymptotics for both [START_REF] Biane | Quantum random walk on the dual of SU(n)[END_REF] and [START_REF] Biane | Minuscule weights and random walks on lattices[END_REF].

For processes with drift, much less is known. Concerning Brownian motion, one of the first significant results is due to Biane, Bougerol and O'Connell [START_REF] Biane | Littelmann paths and Brownian paths[END_REF], who derived the asymptotics of the non-exit probability [START_REF] Biane | Quantum random walk on the dual of SU(n)[END_REF] in the case of Weyl chambers of type A, when the drift is inside of the cone. Later on, by using different techniques, Pucha la and Rolski [START_REF] Pucha La | The exact asymptotic of the collision time tail distribution for independent Brownian particles with different drifts[END_REF] obtained (also in the context of Weyl chambers) the asymptotics of (2) without any hypothesis on the drift. In [START_REF] Garbit | On the exit time from a cone for Brownian motion with drift[END_REF] we gave, for Brownian motion with a given arbitrary drift, the asymptotics of (2) for a large class of cones.

As for random walks (S n ) n 0 with increments having a common distribution µ, the exponential decay of (3) is known: It equals the global minimum on R d of the Laplace transform of µ:

(4) L µ (x) := E µ [e x,S n+1 -Sn ] = R d e x,y µ(dy).
This was first proved by Iglehart [START_REF] Iglehart | Random walks with negative drift conditioned to stay positive[END_REF] for one-dimensional random walks. For more general walks, this was shown and used by many authors (see, e.g., [START_REF] Denisov | Random walks in cones[END_REF][START_REF] Garbit | Temps de sortie d'un cône pour une marche aléatoire centrée[END_REF]). Regarding now the asymptotic behavior of the probability (2), the case d = 1 is known, see [START_REF] Doney | On the asymptotic behaviour of first passage times for transient random walk[END_REF].

It is the aim of this paper to give, for a very broad class of random walks and cones, in any dimension, the exponential decay of the non-exit probability [START_REF] Biane | Quantum random walk on the dual of SU(n)[END_REF]. We shall also relate its value to the Laplace transform, by proving that it equals the minimum of this function on the dual cone; we give the exact statement (Theorem 1) in Subsection 1.4. Our main motivation comes from the possible applications to lattice path enumeration. Indeed, our results provide the first unified treatment of the question of determining the growth constant for the number of lattice paths confined to the positive orthant. They also solve a conjecture on these numbers stated in [START_REF] Johnson | Towards a combinatorial understanding of lattice path asymptotics[END_REF]. However, we would like to emphasize that our results are much more general (see Section 1.3).

Simultaneously and independently of us, Duraj [START_REF] Duraj | Random walks in cones: The case of nonzero drift[END_REF] obtained in some particular case the exact asymptotics of the non-exit probability (2) for lattice random walks. In the following section we introduce our main ideas and tools, and we discuss the difference between our results and his. n ) n 0 be the canonical random walk on R d . Given any probability measure µ on R d and x ∈ R d , we denote by P x µ the probability measure under which (S n ) n 0 is a random walk started at x whose independent increments (S n+1 -S n ) n 0 have common distribution µ.

The standard idea to handle the case of random walks with non-zero drift is to carry out an exponential change of measure. More precisely, if z is a point in R d such that L µ (z) is finite, then we can consider the new probability measure

µ z (dy) = e z,y L µ (z) µ(dy).
It is theoretically possible to compare the behavior of the random walk under P µ with its behavior under P µz thanks to Cramér's formula (see Lemma 4). For example, for the local probabilities, this formula gives

P x µ [S n = y, τ K > n] = L µ (z) n e z,x-y P x µz [S n = y, τ K > n].
Since the asymptotic behavior of those probabilities are now well known when the random walk has no drift (see [START_REF] Denisov | Random walks in cones[END_REF]), the general problem can be solved if one can find a point z such that the distribution µ z is centered. It is also well known that this condition is fulfilled if and only if z is a critical point for L µ (under the assumption that L µ be finite in a neighborhood of z). By convexity of L µ , this means that one has to find a local, hence

global minimum point z = x 0 in R d .
This approach is used by Duraj in [START_REF] Duraj | Random walks in cones: The case of nonzero drift[END_REF] to analyze the non-exit probability [START_REF] Biane | Quantum random walk on the dual of SU(n)[END_REF]. Indeed, for a lattice random walk, one can sum the contribution of each y to eventually obtain (below µ 0 is an abbreviation for µ x 0 )

P x µ [τ K > n] = L µ (x 0 ) n y∈K∩Z d e x 0 ,x-y P x µ 0 [S n = y, τ K > n].
But then, one needs to impose an additional condition on the position of the global minimum point x 0 with respect to K so as to ensure that the infinite sum of asymptotics will be convergent as well. This technical assumption on x 0 done in [START_REF] Duraj | Random walks in cones: The case of nonzero drift[END_REF] happens to have a very natural interpretation in the light of our analysis. Indeed, for the non-exit probability, Cramér's formula (applied with any z) gives ( 5)

P x µ [τ K > n] = L µ (z) n e z,x E x µ 0 [e -z,Sn , τ K > n],
and one sees that the main difficulty will arise because of the exponential term inside the expectation.

Let K * denote the dual cone associated with K, that is, the closed convex cone defined by ( 6)

K * := {z ∈ R d : x, z 0, ∀x ∈ K},
where x, z denotes the standard inner product. If z belongs to K * , it immediately follows from (5) that

P x µ [τ K > n] L µ (z) n e z,x .
Hence, the infimum ρ of the Laplace transform on K * is always an upper bound of the exponential rate, i.e., lim sup n→∞

P x µ [τ K > n] 1/n ρ := inf K * L µ .
Our main result shows that, in fact, when the infimum ρ is a minimum, it is also a lower bound of the above quantity. Thus ρ is the value of the exponential decreasing rate of the non-exit probability. It is now easily seen that assumptions 1 and 5 in [START_REF] Duraj | Random walks in cones: The case of nonzero drift[END_REF] on the global minimum x 0 imply that it belongs to the interior (K * ) o of the dual cone (and in this case, clearly, x 0 = x * ). In [START_REF] Duraj | Random walks in cones: The case of nonzero drift[END_REF] the author then obtains the precise asymptotics of the non-exit probability [START_REF] Biane | Quantum random walk on the dual of SU(n)[END_REF] in this specific case.

The general philosophy of our work is different (and in a sense complementary). We shall only focus on the exponential rate

ρ x := lim sup n→∞ P x µ [τ K > n] 1/n
, and we answer completely the question of determining its value under fairly broad assumptions, regardless the position of the global minimum point x 0 .

1.3. Cones and random walks considered. In this work, we consider a closed convex cone K with non-empty interior. Recall that we denote by K * its dual cone, which turns out to be particularly relevant for our problem. It is the closed convex cone defined in [START_REF] Bousquet-Mélou | Walks with small steps in the quarter plane[END_REF]. We also set

K δ := K + δv,
where δ ∈ R and v is some fixed vector in K o , the interior of K.

Throughout this paper, we shall make the assumption that µ is truly d-dimensional in the following sense:

(H1) The support of the probability measure µ is not included in any linear hyperplane.

For a square-integrable probability measure µ with mean m and variance-covariance matrix Γ, it is well known that the minimal (with respect to inclusion) affine subspace A such that µ(A) = 1 is m + (ker Γ) ⊥ . Hence, the condition in (H1) holds if and only if m + (ker Γ) ⊥ is not included in any hyperplane (or equivalently, if and only if Γ is non-degenerate or dim(ker Γ) = 1 and m / ∈ (ker Γ) ⊥ ). Notice that in the case where m = 0, the assumption (H1) is equivalent to ker Γ = {0}, i.e., Γ is non-degenerate.

As pointed out in the preceding section, our analysis of the decreasing rate of the nonexit probability requires the existence of a minimum point for the Laplace transform L µ on the dual cone. Thus we shall impose the following technical condition: (H2) There exists a point x * ∈ K * and an open neighborhood V of x * in R d such that L µ (x) is finite for all x ∈ V , and x * is a minimum point of L µ restricted to K * ∩ V .

It is worth noting that we do not assume the existence of moments of µ. Hypothesis (H2) implies the existence of these moments only in the case where x * = 0.

In view of applications, we will prove in Subsection 2.3 that for random walks with all exponential moments (i.e., L µ (x) is finite for all x ∈ R d ), the condition (H2) is equivalent to the more geometric-flavoured condition: (H2') The support of µ is not included in any half-space u -:= {x ∈ R d : u, x 0} with u ∈ K * \ {0}.

1.4. Main results. We are now in position to state our main result: Theorem 1. Suppose µ satisfies (H1) and (H2). Then,

lim n→∞ P x µ [τ K > n] 1/n = L µ (x * ),
for all x ∈ K δ , for some constant δ 0.

For a large class of random walks and cones, Theorem 1 gives the universal recipe to compute the exponential decay of the non-exit probability. Notice that the latter is independent of the starting point x. This is not the case when (H2) is not satisfied and we shall illustrate this phenomenon in Section 4, with the walks in the quarter-plane having transition probabilities as in Figure 1.

Let us point out that, in general, there is no explicit link between the position of the drift m of the random walk (if it exists), the position of x * and the value L µ (x * ) of the decreasing rate, except in the case where m belongs to the cone K. As shown in the next lemma, the fact that m ∈ K is a necessary and sufficient condition for having L µ (x * ) = 1 (i.e., a non-exponential decay of the non-exit probability).

Lemma 2. Assume (H1) and (H2). Then L µ (x * ) = 1 if and only if x * = 0. In addition, if the drift m = R d yµ(dy) exists (i.e., if µ admits a moment of order 1), then m belongs to K if and only if x * = 0. Theorem 1 in itself does not provide any explicit value for the constant δ, but such a value can be found a posteriori thanks to the following: Proposition 3. The statement in Theorem 1 holds for any δ 0 for which there exists n 0 1 such that

P 0 µ [τ K -δ > n 0 , S n 0 ∈ K o ] > 0.
Proofs of Theorem 1, Lemma 2 and Proposition 3 are postponed to Subsection 2.2.

In order to illustrate Theorem 1, it is interesting to compare its content with the corresponding result known for Brownian motion with drift, in the light of the recent paper [START_REF] Garbit | On the exit time from a cone for Brownian motion with drift[END_REF]. It is proved there that, for Brownian motion (B t ) t 0 with drift a ∈ R d , the non-exit probability admits the asymptotics (in the continuous case, the exit time from K is defined by τ

K := inf{t > 0 : B t / ∈ K}) (7) P x [τ K > t] = κh(x)t -α e -γt (1 + o(1)), t → ∞,
where γ := d(a, K) 2 /2. Therefore

lim t→∞ P x [τ K > t] 1/t = e -d(a,K) 2 /2 .
Let us compare with the value given by Theorem 1 for the random walk (B n ) n 0 . Its distribution µ is Gaussian with mean a and identity variance-covariance matrix; therefore

L µ (x) = e |x| 2 /2+ x,a . ✲ ✻ | ■ ✠ q p q x Figure 1
. Random walks considered in Section 4 (p + 2q = 1, p, q > 0), for different starting points x

The minimum on K * of |x| 2 /2 + x, a is obviously the minimum on the polar cone

K ♯ = -K * of |x| 2 /2 -x, a = |x -a| 2 /2 -|a| 2 /2.
It is reached at x = p ⊥ K ♯ (a), the orthogonal projection of a on K ♯ , and an easy computation shows that the minimum value is

|p ⊥ K ♯ (a) -a| 2 /2 -|a| 2 /2 = -|a -p ⊥ K (a)| 2 /2 = -d(a, K) 2 /2
, where we have used Moreau's decomposition theorem which asserts that, for any convex cone K, a is the orthogonal sum of p ⊥ K (a) and p ⊥ K ♯ (a). We thus have min

K * L µ = e -d(a,K) 2 /2
, which means that the exponential decreasing rate is the same for Brownian motion (B t ) t 0 and for the "sampled" Brownian motion (B n ) n 0 , as one could expect. 1.5. Plan of the paper. The rest of our article is organized as follows: In Section 2 we prove Theorem 1. In Section 3 we present an important consequence of Theorem 1 in the counting of walks in orthants (a topical domain in enumerative combinatorics), see Corollaries 8 and 9. In Section 4 we consider the walks of Figure 1, for which we prove that contrary to the walks satisfying hypothesis (H2), the exponential decay depends on the starting point x. Finally, in Section 5 we prove the non-exponential decay of the non-exit probability for random walks with drift in the cone (a refinement of a theorem of [START_REF] Garbit | Temps de sortie d'un cône pour une marche aléatoire centrée[END_REF]), which is needed for proving our main result.

Proof of the main results

This section is organized as follows: In Subsection 2.1, we review some elementary properties of the Laplace transform and present Cramér's formula. Then, we prove Theorem 1, Lemma 2 and Proposition 3 in Subsection 2.2. In Subsection 2.3 we provide a geometric interpretation of our main assumption on the random walk distribution. It is clearly a convex function. If L µ is finite in a neighborhood of the origin, say B(0, r), then it is well known that L µ is (infinitely) differentiable in B(0, r), and that its partial derivatives are given by

∂L µ (x) ∂x i = R d y i e x,y µ(dy), ∀i ∈ 1, d .
Therefore, the expectation of µ is equal to the gradient of L µ at the origin:

E[µ] = ∇L µ (0). Notice that µ is centered if and only if 0 is a critical point of L µ .
Suppose now that z is a point where L µ is finite, and let µ z denote the probability measure defined by [START_REF] Denisov | Conditional limit theorems for ordered random walks[END_REF] µ z (dy) := e z,y L µ (z) µ(dy).

The Laplace transform of µ z is related to that of µ by the formula

L µz (x) = L µ (z + x) L µ (z) .
If in addition L µ is finite in some ball B(z, r), then L µz is finite in B(0, r). By consequence, L µz is differentiable in B(0, r) and

E[µ z ] = ∇L µz (0) = ∇L µ (z) L µ (z) .
The distribution of the random walk under P µz is linked to the initial distribution by the following: Lemma 4 (Cramér's formula). For any measurable and positive function F : R n → R, we have

E x µ [F (S 1 , . . . , S n )] = L µ (z) n e z,x E x µz [e -z,Sn F (S 1 , . . . , S n )].
Proof. It follows directly from the definition (8) of µ z that

µ ⊗n (dy 1 , . . . , dy n ) = L µ (z) n e -z, n i=1 y i µ ⊗n z (dy 1 , . . . , dy n ).
The conclusion is then straightforward.

Applied to the function F (s 1 , . . . , s n ) = Π n i=1 ½ K (s i ), Cramér's formula reads

P x µ [τ K > n] = L µ (z) n e z,x E x µz [e -z,Sn , τ K > n],
and this implies that for all x ∈ R d , (9) lim sup n→∞

P x µ [τ K > n] 1/n inf K * L µ ,
as already observed at the end of Subsection 1.2.

2.2. Proofs of Theorem 1, Lemma 2 and Proposition 3. In order to obtain the lower bound for the non-exit probability, we shall use Cramér's formula at z = x * . The following lemma gives some useful information on the position of the drift of the random walk under the measure changed at x * .

Lemma 5. Suppose µ satisfies (H1) and (H2). Then, the gradient ∇L µ (x * ) belongs to K and is orthogonal to x * .

Proof. We first notice that under (H2), the Laplace transform is finite, hence differentiable, in some neighborhood of x * . It is well known that the equality (K * ) * = K holds for any closed convex cone, see [START_REF] Rockafellar | Convex analysis[END_REF]Theorem 14.1]. Hence ∇L µ (x * ) belongs to K if and only if

∇L µ (x * ), y 0, ∀y ∈ K * .
So, let y ∈ K * . Since K * is a convex cone and x * ∈ K * , x * + ty also belongs to K * for all t 0. Hence, thanks to (H2), the function

t ∈ [0, ∞) → f y (t) := L µ (x * + ty)
is differentiable in some neighborhood of t = 0 and reaches its minimum at t = 0. This implies ∇L µ (x * ), y = f ′ y (0) 0. Now, if we take y = x * , we have a stronger result since x * + tx * = (1 + t)x * belongs to K * for all t -1: the function f x * is differentiable in some open neighborhood of t = 0 and has a local minimum point on [-1, ∞) at t = 0, hence

∇L µ (x * ), x * = f ′ x * (0) = 0.
The proof is completed.

We are now in position to conclude the proof of Theorem 1.

Proof of Theorem 1. As already observed, L µ (x * ) is an upper bound for the exponential decreasing rate lim sup

n→∞ P x µ [τ K > n] 1/n ,
and it remains to prove that it is also the right lower bound. By performing the Cramér transformation at x * , we get

P x µ [τ K > n] = ρ n e x * ,x E x µ * [e -x * ,Sn , τ K > n] ρ n e x * ,x e -α √ n P x µ * [| x * , S n | α √ n, τ K > n],
where ρ = L µ (x * ), µ * (dy) = ρ -1 e x * ,y µ(dy), and α is any positive number. Notice that µ * is truly d-dimensional, because µ has this property and both measures have the same support. Note also that by assumption (H2), L µ * (x) = ρ -1 L µ (x * + x) is finite in some neighborhood of x = 0, and therefore µ * has all moments. Since the new drift

E[µ * ] = m * = ρ -1 ∇L µ (x * )
belongs to K and is orthogonal to x * (by Lemma 5), it follows from Proposition 12 in Section 5 that there exist α > 0 and δ 0 such that

lim n→∞ P x µ * [| x * , S n | α √ n, τ K > n] 1/n = 1, ∀x ∈ K δ .
Hence, we reach the conclusion that lim inf n→∞

P x µ [τ K > n] 1/n ρ
for all x ∈ K δ , and the theorem is proved.

We now give the proof of Lemma 2, which provides a necessary and sufficient condition for having L µ (x * ) = 1.

Proof of Lemma 2. Assume that L µ (x * ) = 1 and at the same time x * = 0. It is well known that L µ is then finite on [0, x * ], thus strictly convex on that segment (see Subsection 2.3). Since L µ (0) = L µ (x * ) = 1 it follows that L µ (x) < 1 for all x ∈ (0, x * ). But this open interval is a subset of K * , hence this contradicts the hypothesis (H2) asserting that x * is a local minimum point on K * . Conversely, that x * = 0 implies L µ (x * ) = 1 is trivial.

We now turn to the second part of the lemma. First, we know from Lemma 5 that ∇L µ (x * ) belongs to K. So, if x * = 0, then m = ∇L µ (0) belongs to K (here we do not need to assume the existence of m: it exists because (H2) at x * = 0 ensures that L µ is infinitely differentiable in some neighborhood of 0, and therefore µ has all moments). Conversely, assume that m exists and belongs to K and suppose that x * = 0. Consider the function g(t) = L µ (tx * ), which is finite on [0, 1]. Under the assumption |y|µ(dy) < ∞, it follows by standard arguments that g(t) has a right derivative at t = 0 given by g ′ (0+) = x * , m . Since m belongs to K and x * to K * , this derivative is non-negative. So, g(t) must be increasing since it is strictly convex. Thus L µ (tx * ) = g(t) < g(1) = L µ (x * ) for all t ∈ [0, 1), and x * cannot be a local minimum.

To conclude this section, we explain how to find a δ 0 for which the statement of Theorem 1 holds.

Proof of Proposition 3. Recall that v ∈ K o is fixed and K δ = K + δv. We assume that there exist δ 0 and k 1 such that

P 0 µ [τ K -δ > k, S k ∈ K o ] > 0.
Therefore, we can find ǫ > 0 such that

P 0 µ [τ K -δ > k, S k ∈ K ǫ ] = γ > 0,
and since K is a convex cone, it satisfies the relation K + K ⊂ K, thus

P x µ [τ K > k, S k -x ∈ K ǫ ] γ,
for all x ∈ K δ (by inclusion of events). From this, we shall deduce by induction that (10)

p ℓ := P x µ [τ K > ℓk, S ℓk -x ∈ K ℓǫ ] γ ℓ ,
Let S d-1 denote the unit sphere of R d . Standard arguments involving the convexity of L µ and the compactness of C ∩S d-1 show that L µ (x) goes to infinity as |x| → ∞ uniformly on C if and only if [START_REF] Duraj | Random walks in cones: The case of nonzero drift[END_REF] lim

t→∞ L µ (tu) = ∞, ∀u ∈ C ∩ S d-1 .
Hence, the condition ( 12) is sufficient for the existence of a global minimum on C. Indeed, if it is satisfied, then there exists R > 0 such that L µ (x) 1 for all x ∈ C with |x| R. By continuity, L µ reaches a minimum on B(0, R) ∩ C which is less than or equal to 1 (since L µ (0) = 1), thus it is a global minimum on C.

The next lemma gives some interesting information on the behavior at infinity of L µ . Recall that u -denotes the half-space {y ∈ R d : u, y 0}.

Lemma 6. Suppose that µ has all C-exponential moments. For every u ∈ C ∩ S d-1 , the following dichotomy holds:

(1) If µ(u -) < 1, then lim t→∞ L µ (x + tu) = ∞, ∀x ∈ C. (2) If µ(u -) = 1, then lim t→∞ L µ (x + tu) = u ⊥ e x,y µ(dy), ∀x ∈ C.
Proof. If µ(u -) < 1, then we can find ǫ > 0 such that the set {y ∈ R d : u, y ǫ} has positive measure, and the inequality L µ (x + tu) {y∈R d : u,y ǫ} e x,y e tǫ µ(dy) ce tǫ , proves the first assertion of Lemma 6. Suppose now on the contrary that µ(u -) = 1. We then may write

L µ (x + tu) = u ⊥ e x,y µ(dy) + {y∈R d : u,y <0}
e x+tu,y µ(dy).

The second integral on the right-hand side of the above equation goes to zero as t goes to infinity by the dominated convergence theorem, thus proving the second assertion of Lemma 6.

Lemma 7. Suppose that µ satisfies (H1) and has all C-exponential moments. Then the Laplace transform L µ has a global minimum on the closed convex cone C if and only if there does not exist any u = 0 in C such that µ(u -) = 1.

Proof. If µ(u -) < 1 for all u = 0 in C, then lim t→∞ L µ (tu) = ∞ by Lemma 6, and the function L µ has a global minimum on C as explained earlier. Now, suppose on the contrary that µ(u -) = 1 for some u = 0 in C. Then, by Lemma 6 again, the limit h(x) := lim t→∞ L µ (x + tu) exists and is finite for all x. Since any convex cone is a semi-group, x + tu ∈ C for all x ∈ C and t 0, and consequently

h(x) inf C L µ , ∀x ∈ C.
But our assumption that µ satisfies (H1) implies that L µ is strictly convex, so that L µ (x) > h(x) for all x ∈ C (for else the strictly convex function t → L µ (x + tu) on [0, ∞) would not have a finite limit). We thus reach the conclusion that L µ (x) > inf C L µ for all x ∈ C, thereby proving that L µ has no global minimum on C.

For random walks with all exponential moments, the equivalence between conditions (H2) and (H2') under (H1) is now an easy consequence of Lemma 7.

Application to lattice path enumeration

In this section we present an application of our main result (Theorem 1) in enumerative combinatorics: Given a finite set S of allowed steps, a now classical problem is to study S-walks in the orthant

Q := (R + ) d = {x ∈ R d : x i 0, ∀i ∈ 1, d },
that is walks confined to Q, starting at a fixed point x (often the origin) and using steps in S only. Denote by f S (x, y; n) the number of such walks that end at y and use exactly n steps. Many properties of the counting numbers f S (x, y; n) have been recently analyzed (the seminal work in this area is [START_REF] Bousquet-Mélou | Walks with small steps in the quarter plane[END_REF]). First, exact properties of them were derived, via the study of their generating function (exact expression and algebraic nature). Such properties are now well established for the case of small steps walks in the quarter-plane, meaning that the step set S is included in {0, ±1} 2 . More qualitative properties of the f S (x, y; n) were also investigated, such as the asymptotic behavior, as n → ∞, of the number of excursions f S (x, y; n) for fixed y, or that of the total number of walks, [START_REF] Dyson | A Brownian-motion model for the eigenvalues of a random matrix[END_REF] f S (x; n) := y∈Q f S (x, y; n).

Concerning the excursions, several small steps cases have been treated by Bousquet-Mélou and Mishna [START_REF] Bousquet-Mélou | Walks with small steps in the quarter plane[END_REF] and by Fayolle and Raschel [START_REF] Fayolle | Some exact asymptotics in the counting of walks in the quarter plane[END_REF]. Later on, Denisov and Wachtel [START_REF] Denisov | Random walks in cones[END_REF] obtained the very precise asymptotics of the excursions, for a quite large class of step sets and cones. As for the total number of walks ( 13), only very particular cases are solved, see again [START_REF] Bousquet-Mélou | Walks with small steps in the quarter plane[END_REF][START_REF] Fayolle | Some exact asymptotics in the counting of walks in the quarter plane[END_REF]. In a most recent work [START_REF] Johnson | Towards a combinatorial understanding of lattice path asymptotics[END_REF], Johnson, Mishna and Yeats obtained an upper bound for the exponential growth constant, namely,

lim sup n→∞ f S (x; n) 1/n ,
and proved by comparison with results of [START_REF] Fayolle | Some exact asymptotics in the counting of walks in the quarter plane[END_REF] that these bounds are tight for all small steps models in the quarter-plane. In the present article, we find the exponential growth constant of the total number of walks [START_REF] Dyson | A Brownian-motion model for the eigenvalues of a random matrix[END_REF] in any dimension for any model such that:

(H1") The step set S is not included in a linear hyperplane; (H2") The step set S is not included in a half-space u -, with u ∈ Q \ {0}.

Our results provide the first unified treatment of this problem of determining the growth constant for the number of lattice paths confined to the positive orthant. In the sequel we shall say that a step set S is proper if it satisfies to (H1") and (H2"). Note in particular that the well-known 79 models of walks in the quarter-plane studied in [START_REF] Bousquet-Mélou | Walks with small steps in the quarter plane[END_REF][START_REF] Fayolle | Some exact asymptotics in the counting of walks in the quarter plane[END_REF] (including the so-called 5 singular walks) satisfy both hypotheses above.

Corollary 8. Let S be any proper step set. The Laplace transform of S,

L S (x) := s∈S e x,s
, reaches a global minimum on Q at a unique point x 0 , and there exists δ 0 such that for any starting point

x ∈ Q δ , lim n→∞ f S (x; n) 1/n = L S (x 0 ).
Suppose that in addition:

(H3") The step set allows a path staying in Q from the origin to some point in the interior of Q.

Then it follows from Proposition 3 that the result in Corollary 8 holds with δ = 0, i.e., it is valid for all x ∈ Q. Note that this assumption is not restrictive from a combinatorial point of view, since if (H3") is not satisfied, the counting problem is obvious.

Proof of Corollary 8. Consider a random walk (S n ) n 0 starting from x such that µ is the uniform law on S. Let then τ Q denote the first exit time from Q. The enumeration problem is related to probabilities in a simple way:

(14) P x µ [τ Q > n] = f S (x; n) |S| n .
Further, it is immediate from our definitions that L S (x) = |S|L µ (x). Corollary 8 then follows from Theorem 1 and from the fact that Q * = Q.

As a consequence, we obtain the following result, which was conjectured in [START_REF] Johnson | Towards a combinatorial understanding of lattice path asymptotics[END_REF]:

Corollary 9. Let S ⊂ Z d be a proper step set (hypotheses (H1") and (H2")), which additionally satisfies (H3"), and let K S be the growth constant for the total number of walks [START_REF] Dyson | A Brownian-motion model for the eigenvalues of a random matrix[END_REF]. Let P be the set of hyperplanes through the origin in R d which do not meet the interior of the first orthant. Given p ∈ P, let K S (p) be the growth constant of the walks on S which are restricted to the side of p which includes the first orthant. Then K S = min p∈P K S (p).

Proof. Let us first notice that P can be described as the set of hyperplanes u ⊥ such that u ∈ Q ∩ S d-1 , and that the side of p = u ⊥ which includes the first orthant is then the half-space u + = {x ∈ R d : x, u 0}. By Theorem 1, the exponential rate for the random walk associated to the step set S and confined to u + is the minimum of L S /|S| on the dual cone (u + ) * = {tu : t 0}. Therefore, the growth constant K S (p) equals min t 0 L S (tu), and the equality min

x∈Q L S (x) = min u∈Q∩S d-1 min t 0 L S (tu) ✲ ✻ • • • • • • • • • • • • • • • • | ■ ✠ 2N 2N + 2 2N 2N + 2 0 q q • p Figure 2.
Random walks considered in the proof of Proposition 10 on the lines {(i, j)

∈ Q : i + j = 2N } for N 0 immediately translates into K S = min p∈P K S (p).
The proof of Corollary 9 is completed.

An example of half-space walks

In this section we illustrate the following phenomenon: If the support of the random walk is included in a certain half-space (chosen so as contradicting (H2)), such a universal result as Theorem 1 does not hold; in particular, the exponential decay of the non-exit probability may depend on the starting point.

Let (S n ) n 0 be the random walk on Q starting at x and with transition probabilities to (1, -1), (-1, 1) and (-1, -1), with respective probabilities q, q and p, where p + 2q = 1 and p, q > 0, see Figures 1 and2. Let τ Q be the exit time (1) of this random walk from the quarter-plane Q. Finally, define for fixed N the segment

D 2N = {(i, j) ∈ Q : i + j = 2N }.
Proposition 10. For any N 1 and any x ∈ D 2N , we have [START_REF] Fayolle | Some exact asymptotics in the counting of walks in the quarter plane[END_REF] lim

n→∞ P x [τ Q > n] 1/n = 2q cos π 2N + 2 .
We shall need the following result on the simple symmetric random walk on Z (the proof can be easily derived from the identities in [22, page 243]): Lemma 11 ([22]). For the simple symmetric random walk ( S n ) n 0 on Z (with jumps to the left and to the right with equal probabilities 1/2), we have, for any x ∈ 0, 2N ,

lim n→∞ P x [ S 1 , . . . , S n ∈ 0, 2N ] 1/n = cos π 2N + 2 .
Proof of Proposition 10. We shall prove Proposition 10 by induction over N 1. For N = 1, we have three choices for x (see Figure 2). We write the proof when x = (1, 1), since the arguments for other values of x are quite similar. For this choice of x, the origin (0, 0) can be reached only at odd times n, and in that event, the random walk gets out of Q at time n + 1. From this simple remark we deduce that (below, we note (X k ) k 1 the increments of the random walk (S n ) n 0 )

P x [τ Q > 2n] = P x [τ Q > 2n, X k = (-1, -1), ∀k ∈ 1, n ] = P x [τ Q > 2n| X k = (-1, -1), ∀k ∈ 1, n ](2q) n .
Further, the random walk conditioned on never making the jump (-1, -1) is a simple symmetric random walk on the segment D 2 . Therefore,

P x [τ Q > 2n] = P 1 [ S 1 , . . . , S n ∈ 0, 2 ](2q) 2n .
Using Lemma 11, we conclude that for x = (1, 1),

lim n→∞ P x [τ Q > 2n] 1/(2n) = 2q cos π 4 .
The fact that P x [τ > n] is decreasing in n implies that the above equation holds with 2n + 1 instead of 2n. This achieves the proof of Proposition 10 for N = 1.

Let us now assume that equation ( 15) holds for a fixed value of N 1. For

x ∈ D 2N +2 , introduce H := inf{n > 0 : S n ∈ D 2N }
the hitting time of the set D 2N , see Figure 2. We can write ( 16)

P x [τ Q > n] = P x [τ Q > n, H > n] + n k=1 P x [τ Q > n, H = k].
The first term in the right-hand side of ( 16) can be written as

P x [τ Q > n, H > n] = P x [τ Q > n|H > n](2q) n ,
where (for the same reasons as for the case N = 1)

lim n→∞ P x [τ Q > n|H > n] 1/n = cos π 2N + 4 .
As for the second term in the right-hand side of ( 16),

P x [τ Q > n, H = k] = E x [τ Q > k, H = k, P S k [τ Q > n -k]] C • P x [τ Q > k, H = k]P x 0 [τ Q > n -k] C • P x [τ Q > k -1, H > k -1]P x 0 [τ Q > n -k] := Ca k b n-k .
The first equality above comes from the strong Markov property. The first inequality follows from the fact that for any fixed x 0 ∈ D 2N , there exists a constant C > 0 such that, for any n 0 and any y

∈ D 2N , P y [τ Q > n] CP x 0 [τ Q > n].
The second inequality is obvious, and the last line has to be read as a definition.

Using on the one hand the same reasoning as for the case N = 1, and on the other hand the induction hypothesis, we obtain To summarize, with the help of ( 16) we have written

P x [τ Q > n] = A n + B n , where lim n→∞ A 1/n n = 2q cos π 2N + 4 , lim sup n→∞ B 1/n n 2q cos π 2N + 4 .
The formula (15) therefore holds for N + 1, and Proposition 10 is proved.

The case of random walks with drift in the cone

In this section we refine a result announced in [START_REF] Garbit | Temps de sortie d'un cône pour une marche aléatoire centrée[END_REF] concerning the non-exponential decay of the non-exit probability (2) from a cone for a centered, square-integrable and non-degenerate multidimensional random walk. We prove that the result in [START_REF] Garbit | Temps de sortie d'un cône pour une marche aléatoire centrée[END_REF] still holds if the random walk has a drift in the cone and if the hypothesis that its variance-covariance matrix is non-degenerate is weakened to (H1). This result is one of the main ingredients of the proof of Theorem 1. We would like to notice that the proof of (a weakened form of) Proposition 12 below was only sketched in [START_REF] Garbit | Temps de sortie d'un cône pour une marche aléatoire centrée[END_REF], so that the extended proof we shall give here is new.

As before, we only assume that K is a closed convex cone with non-empty interior K o . We fix some v ∈ K o and define K δ = K + δv. In this setting, we shall prove the following: Proposition 12. Assume that the distribution µ of the random walk increments is squareintegrable and truly d-dimensional (H1). Suppose in addition that the drift m = E 0 [S 1 ] belongs to the cone K and that v is a vector orthogonal to m. Then there exists α > 0 and δ 0 such that, for all x ∈ K δ ,

lim n→∞ P x [τ K > n, | v, S n | α √ n] 1/n = 1.
If m = 0, as we have already explained, (H1) is equivalent to the fact that the variancecovariance matrix of the increments distribution is non-degenerate. Hence, we exactly recover [START_REF] Garbit | Temps de sortie d'un cône pour une marche aléatoire centrée[END_REF]Theorem]. However, if m = 0, Proposition 12 can not be derived from [START_REF] Garbit | Temps de sortie d'un cône pour une marche aléatoire centrée[END_REF]Theorem]. Indeed, under the hypothesis of Theorem 12, it is clear that

P x [τ K > n, | v, S n | α √ n] P x [τ K ( S) > n, | v, S n | α √ n],
where ( S n = S nnm) n 0 is the centered random walk associated with (S n ) n 0 . But the variance-covariance matrix of S 1 is equal to that of S 1 and might be degenerate, so that [START_REF] Garbit | Temps de sortie d'un cône pour une marche aléatoire centrée[END_REF]Theorem] would not apply to the walk ( S n ) n 0 . This is for example the case when (S n ) n 0 is the uniform two-dimensional random walk with step set S = {(0, 1), (1, 0)}. In order to prove Theorem 12, we will need a series of lemmas. We begin with some geometric considerations. Lemma 13. Let m be a point in K and V be a linear subspace of R d . If (m+V )∩K o = ∅, then m + V is included in a hyperplane.

Proof. Clearly, dim V < d, for else m + V = R d would intersect K o . Since m + V is included in the linear subspace generated by m and V , whose dimension is less or equal to dim V + 1, the affine space m + V is always included in a hyperplane if dim V < d -1. Thus, it remains to consider the case where V is a hyperplane, that is,

V = u ⊥ = {x ∈ R d : u, x = 0}
for some u = 0. Assume m / ∈ V . Possibly changing u to -u, we can assume in addition that m, u > 0. Then,

λ>0 (λm + V ) = λ>0 (λu + V ) = {x ∈ R d : u, x > 0} =: u + * .
Further, by homogeneity of V and K o , we have (λm + V )∩ K o = ∅ for every λ > 0. Hence, K o does not intersect u + * , and is therefore included in u -= {x ∈ R d : u, x 0}. This is a contradiction since m ∈ K = (K o ) (this equality holds for any convex set with non-empty interior) and m, u > 0. Thus m belongs to V and m + V = V is a hyperplane. Lemma 14. Let Y ∈ R d be a random vector with Gaussian distribution N (m, Γ). If m belongs to K and if m + (ker Γ) ⊥ is not included in a hyperplane, then

P[Y ∈ K o ] > 0.
Proof. It is well known that the Gaussian distribution N (m, Γ) admits a positive density with respect to Lebesgue measure on the affine space m + (ker Γ) ⊥ . Thus, it suffices to show that (m + (ker Γ) ⊥ ) ∩ K o is a non-empty open set in m + (ker Γ) ⊥ . But this follows from Lemma 13 since m + (ker Γ) ⊥ is not contained in a hyperplane. Lemma 15. Under the hypotheses of Theorem 12, there exist k 1 and δ 0 such that

P[τ K -δ > k, S k ∈ K o ] > 0.
Proof. For n 1 define Z n := (S nnm)/ √ n. The homogeneity and convexity properties of K ensure that

{S n ∈ K o } = {Z n ∈ K o - √ nm} ⊃ {Z n ∈ K o -m}
Let Γ denote the variance-covariance matrix of µ. We notice that m + (ker Γ) ⊥ is not included in a hyperplane since µ satisfies (H1). The central limit theorem asserts that (Z n ) n 1 converges in distribution to a random vector Y with Gaussian distribution N (0, Γ). Hence, applying the Portmanteau theorem [5, Theorem 2.1], we obtain the inequality lim inf

n→∞ P[S n ∈ K o ] P[Y ∈ K o -m],
where the right-hand side is positive according to Lemma 14. Now, to conclude it suffices to fix k so that P[S k ∈ K o ] = 2ǫ > 0, and then choose δ so large that P

[τ K -δ > k] 1 -ǫ (this is possible since K -δ ↑ R d as δ ↑ ∞).
Proof of Theorem 12. The proof follows the same kind of arguments as in [START_REF] Garbit | Temps de sortie d'un cône pour une marche aléatoire centrée[END_REF]. We shall first use Lemma 15 in order to push the random walk inside the cone at a distance √ n from the boundary, and then apply the functional central limit theorem.

Recall that v ∈ K o is fixed and that K δ = K + δv. We know by Lemma 15 that there exist δ 0 and k 1 such that

P[τ K -δ > k, S k ∈ K o ] > 0.
Therefore, we can find a closed ball B := B(z, ǫ) ⊂ K o , with center at z ∈ K o and radius ǫ > 0, such that

P[τ K -δ > k, S k ∈ B] = γ > 0. Since K is a convex cone, it satisfies the relation K + K ⊂ K, thus P x [τ K > k, S k -x ∈ B] γ,
for all x ∈ K δ (by inclusion of events). From this, we shall deduce by induction that [START_REF] Garbit | On the exit time from a cone for Brownian motion with drift[END_REF] p

ℓ := P x [τ K > ℓk, S ℓk -x ∈ ℓB] γ ℓ ,
for all ℓ 1 and x ∈ K δ . Indeed, by the Markov property of the random walk, which we shall prove now. Since ℓk ≪ n does not play any significant role in the last probability, we will neglect it in order to simplify notations. Also, for any ǫ ′ > ǫ, we have

p ℓ+1 E x [τ K > ℓk, S ℓk -x ∈ ℓB, P S ℓk [τ K > k, S k -x ∈ (ℓ + 1)B]]
x + ℓB = ℓ x ℓ + B(z, ǫ) ⊂ ℓB(z, ǫ ′ ) for all large enough n. Since an ǫ ′ > ǫ can be found so that B(z, ǫ ′ ) ⊂ K o , we may replace x + ℓB by ℓB without loss of generality. Finally, since ℓ √ n, we may replace ℓB by √ nB. With these simplifications, it remains to consider

q n := inf {y∈K:y∈ √ nB} P y [τ K > n, | v, S n | α √ n].
By mapping y to y/ √ n, we may write q n = inf y∈B q n (y), where q n (y) := P 0 [τ K (y

√ n + S) > n, | v, y √ n + S n | α √ n].
Let S = ( S n = S n -nm) n 0 denote the centered random walk associated with S = (S n ) n 0 . By inclusion of events we get the lower bound

q n (y) P 0 [τ K (y √ n + S) > n, | v, y √ n + S n | α √ n],
where we used the fact that m ∈ K and v, m = 0. Finally, let us denote by Z n = (Z n (t)) t the random process with continuous paths that coincides with S k / √ n for t = k/n and which is linearly interpolated elsewhere. By definition of Z n and convexity of K, the last inequality immediately rewrites q n (y) P 0 [τ K (y + Z n (t)) > 1, | v, y + Z n (1) | α].

The functional central limit theorem ensures that Z n converges in distribution to a Brownian motion (b(t)) t with variance-covariance matrix Γ. Suppose that the sequence (y n ) n 0 converges to some y ∈ B. Then (y n + Z n ) n 0 converges in distribution to y + b, and it follows from the Portmanteau theorem that Therefore, the probability in ( 18) is bounded from below by the probability

P max t∈[0,1]
|b(t)| < ηǫ that the Brownian motion (b(t)) t stays near the origin for all t ∈ [0, 1], and this event happens with positive probability, regardless Γ be positive definite or not.

To summarize, we have proved that lim inf n→∞ q n (y n ) > 0, for any sequence (y n ) n 0 ∈ B that converges to some y. Thus, by standard compactness arguments, we reach the conclusion that lim inf n→∞ q n = lim inf n→∞ inf y∈B q n (y) > 0, and Proposition 12 is proved.

1. 2 .

 2 Preliminary discussion. Let (S n ) n 0 = (S

2. 1 .

 1 Cramér's formula. The Laplace transform of a probability distribution µ is the function L µ defined for x ∈ R d by L µ (x) := R d e x,y µ(dy).

  p ℓ • inf {y∈K:y-x∈ℓB} P y [τ K > k, S kx ∈ (ℓ + 1)B].Noticing that yx ∈ ℓB and S ky ∈ B yields S kx ∈ (ℓ + 1)B, we obtainp ℓ+1 p ℓ • inf {y∈K:y-x∈ℓB} P y [τ K > k, S ky ∈ B].But x ∈ K δ and yx ∈ ℓB ⊂ K imply y ∈ K δ . Hence p ℓ+1 p ℓ • γ and (17) is proved. Now, for x ∈ K δ , definep n := P x [τ K > n, | v, S n | α √ n],where α > 0 will be fixed latter. Write ℓ = ⌊ √ n⌋ for the lower integer part of √ n. Using the Markov property at time ℓk and the estimate in (17) leads top n P x [τ K > n, S ℓkx ∈ ℓB, | v, S n | α √ n] E x [τ K > ℓk, S ℓkx ∈ ℓB, P S ℓk [τ K > nℓk, | v, S n-ℓk | α √ n]] γ ℓ • inf {y∈K:y-x∈ℓB} P y [τ K > nℓk, | v, S n-ℓk | α √ n].Therefore, Proposition 12 will follow from the fact that lim inf n→∞ inf {y∈K:y-x∈ℓB}P y [τ K > nℓk, | v, S n-ℓk | α √ n] > 0,

  y n ) P 0 [τ K o (y + b(t)) > 1, | v, y + b(1) | < α].Now, it is time to choose α. To do this, first recall that B = B(z, ǫ) ⊂ K o . Choose η > ǫ so that B(z, η) ⊂ K o and set (notice that we could have done this at the very beginning of the proof) α = |v|(|z| + η). If |b(t)| < ηǫ for all t ∈ [0, 1], then y + b(t) ∈ B(z, η) ⊂ K o for all t ∈ [0, 1]. Furthermore, | v, y + b(1) | |v|(|y| + |b(1)|) < |v|(|z| + ǫ + ηǫ) = α.

Acknowledgments

We thank Marni Mishna for motivating discussions. Many thanks also to Marc Peigné for his valuable comments and encouragement. We thank two anonymous referees for useful comments and suggestions.

for all ℓ 1 and x ∈ K δ . Indeed, by the Markov property of the random walk,

where the infimum is taken over all y ∈ K such that yx ∈ K ℓǫ . Noting that yx ∈ K ℓǫ and S ky ∈ K ǫ imply S kx ∈ K (ℓ+1)ǫ , we obtain

But x ∈ K δ and yx ∈ K ℓǫ ⊂ K imply y ∈ K δ . Hence p ℓ+1 p ℓ • γ and ( 10) is proved. Now Theorem 1 asserts the existence of some δ 0 0 such that

for all y ∈ K δ 0 , and we want to prove that the result also holds for x ∈ K δ . To do this, we shall simply use [START_REF] Doney | On the asymptotic behaviour of first passage times for transient random walk[END_REF] in order to push the walk from K δ to K δ 0 . More precisely, choose ℓ 1 such that δ + ℓǫ δ 0 . Then for all x ∈ K δ the inclusion x + K ℓǫ ⊂ K δ 0 holds, and thanks to [START_REF] Doney | On the asymptotic behaviour of first passage times for transient random walk[END_REF],

By the Markov property, for all n m, we have

where the last inequality follows by inclusion of events. This implies immediately lim inf

But since the inequality lim sup n→∞

holds for all x (see ( 9)), Proposition 3 is proved.

Geometric interpretation of condition (H2

). The aim of this subsection is to give a geometric interpretation of condition (H2) under some additional condition on the exponential moments. Throughout this section, we fix a closed convex cone C and assume that µ has all C-exponential moments, that is,

The strict convexity of the exponential function ensures that (11) L µ (ax 1 + bx 2 ) aL µ (x 1 ) + bL µ (x 2 ), for all x 1 = x 2 ∈ C and a, b > 0 with a + b = 1, and that equality occurs if and only if

where (x 1x 2 ) ⊥ denotes the hyperplane orthogonal to x 1x 2 . Thus, if (H1) is satisfied, the equality in [START_REF] Doumerc | Exit problems associated with finite reflection groups[END_REF] never occurs and L µ is strictly convex on C.


Olivier Commowick
email: olivier.commowick@inria.fr

Adil Maarouf

Jean-Christophe Ferré

Jean-Philippe Ranjeva

Gilles Edan

Christian Barillot

Diffusion MRI Abnormalities Detection with Orientation Distribution Functions: a Multiple Sclerosis Longitudinal Study

Keywords: Diffusion MRI, patient to controls comparison, Orientation distribution functions

We propose a new algorithm for the voxelwise analysis of orientation distribution functions between one image and a group of reference images. It relies on a generic framework for the comparison of diffusion probabilities on the sphere, sampled from the underlying models. We demonstrate that this method, combined to dimensionality reduction through a principal component analysis, allows for more robust detection of lesions on simulated data when compared to classical tensor-based analysis. We then demonstrate the efficiency of this pipeline on the longitudinal comparison of multiple sclerosis patients at an early stage of the disease: right after their first clinically isolated syndrome (CIS) and three months later. We demonstrate the predictive value of ODF-based scores for the early detection of lesions that will appear or heal.

Introduction

Diffusion MRI has been widely used in recent studies for the detection of abnormalities between brains of patients suffering from different diseases and those of healthy controls. To do so, several methodologies have been proposed that may be split into two categories. First, several groups [START_REF] Whitcher | Statistical group comparison of diffusion tensors via multivariate hypothesis testing[END_REF]; [START_REF] Lepore | Generalized tensor-based morphometry of HIV/AIDS using multivariate statistics on deformation tensors[END_REF]) have presented studies comparing a population of patient brains with respect to a population of healthy controls. Such approaches allow to see, at the population level, how a disease affects the brain white matter structure and are therefore of great interest to understand globally a disease. The second category focuses more on the study of a specific patient with respect to a population of controls [START_REF] Commowick | Detection of DTI white matter abnormalities in multiple sclerosis patients[END_REF]), allowing to detect for that specific person where the brain may be affected (i.e. presence of lesions or other diffuse pathology-related change).

All these methods have been widely used to perform studies on multiple sclerosis (MS) with the goal of better understanding the disease and possibly find early biomarkers of the pathology evolution that would allow an adaptation of the treatment to the patient. Among those studies, [START_REF] Filippi | Diffusion tensor magnetic resonance imaging in multiple sclerosis[END_REF] studied, in different regions of interest, differences of diffusion derived scalar coefficients (fractional anisotropy -FA, apparent diffusion coefficient -ADC) between controls and MS patients. [START_REF] Mani | Morphological changes in the corpus callosum: A study using a joint Riemannian feature spaces[END_REF] on the other hand performed studies on the shape of specific white matter fiber tracts, corpus callosum, and found evidence of shape changes related to multiple sclerosis. Longitudinal studies of MS in diffusion tensor imaging have also been performed, trying to find correlation between clinical scores such as EDSS (Expanded Disability Status Scale) and DTI derived scalar values. [START_REF] Giorgio | Relationships of brain white matter microstructure with clinical and MR measures in relapsingremitting multiple sclerosis[END_REF] or [START_REF] Liu | Whole brain white matter changes revealed by multiple diffusion metrics in multiple sclerosis: A TBSS study[END_REF] performed tract-based spatial statistics and found some correlation in several white matter regions, such as the splenium of the corpus callosum, between decreasing FA and EDSS in relapsing remitting patients. [START_REF] Rocca | Weekly diffusion-weighted imaging of normal-appearing white matter in MS[END_REF] conducted an interesting study on the predictive value of diffusion MRI in a relatively small number of patients with scans every week, showing an increase in mean diffusivity in normal appearing white matter areas subsequently transforming into lesions. [START_REF] Werring | The pathogenesis of lesions and normal-appearing white matter changes in multiple sclerosis: a serial diffusion MRI study[END_REF] studied the evolution of ADC in relationship with the apparition of lesions with contrast in Gadolinium-enhanced MR images, concluding that structural damage in lesions can cause damage in connected areas of normal appearing white matter. Finally, [START_REF] Grigis | Longitudinal change detection in diffusion MRI using multivariate statistical testing on tensors[END_REF] presented a longitudinal study of change detection in diffusion MRI, based on DWI signal bootstrap, demonstrating the ability to detect lesions or diffuse damage as a change of diffusion properties. Only few studies have looked at diffusion imaging in multiple sclerosis from the beginning of the disease (i.e. from the first clinically isolated syndrome -CIS). Such studies could bring important insights into the early development of the disease, leading to a better patient management in daily practice. Almost all previously mentioned MS studies have considered the comparison of scalar values derived from DTI, such as ADC or FA, and are thereby losing part of the information contained in the tensor. [START_REF] Whitcher | Statistical group comparison of diffusion tensors via multivariate hypothesis testing[END_REF]; [START_REF] Lepore | Generalized tensor-based morphometry of HIV/AIDS using multivariate statistics on deformation tensors[END_REF] showed in studies on other pathologies that considering the whole tensor can be more sensitive to detect differences. Moreover, it has been demonstrated that the tensor model is not sufficient to capture regions of crossing fibers in the brain, as it allows to code only for one principal direction of diffusion when there may be several. Recently, several papers proposed higher order diffusion models able to capture several diffusion directions [START_REF] Assaf | Composite hindered and restricted model of diffusion (CHARMED) MR imaging of the human brain[END_REF]; [START_REF] Descoteaux | Regularized, Fast, and Robust Analytical Q-Ball Imaging[END_REF]; [START_REF] Scherrer | Parametric Representation of Multiple White Matter Fascicles from Cube and Sphere Diffusion MRI[END_REF]; [START_REF] Stamm | A new multi-fiber model for low angular resolution diffusion MRI[END_REF]; [START_REF] Zhang | NODDI: practical in vivo neurite orientation dispersion and density imaging of the human brain[END_REF]). Utilizing the full information coming from those diffusion models could potentially improve further the detection power of abnormalities in patients. Following this idea, recent works have studied frameworks based on multiple fascicle models. [START_REF] Jbabdi | Crossing fibres in tract-based spatial statistics[END_REF] presented an extension of the popular tractbased statistics framework to better handle fiber crossing regions, based on a ball-and-stick model. [START_REF] Taquet | A Mathematical Framework for the Registration and Analysis of Multi-Fascicle Models for Population Studies of the Brain Microstructure[END_REF] presented a framework to perform group comparison between patients and controls relying on the multi-tensor model. Both approaches compare information along tracts of interest at the group level to detect any abnormalities related to a disease. However, as group comparison frameworks, these papers are not directly applicable to a single patient to help predict disease evolution. Moreover, they consider tract-based statistics while in prospective studies we may be also interested in detecting differences at the voxel level over the whole brain.

We propose a new method to detect abnormalities by comparing the orientation distribution functions (ODF) at the voxel level between a patient and a population of controls, in a similar fashion to [START_REF] Commowick | Detection of DTI white matter abnormalities in multiple sclerosis patients[END_REF]. It relies on a statistical comparison, combined to a principal component analysis (PCA) for robustness, of the diffusion probability values sampled on the sphere from virtually any model that has a probability density function defined on the sphere. We applied this method first to simulated data with lesions to evaluate the detection power of ODF-based scores with respect to DTI-based scores. Then, we applied our pipeline to the longitudinal detection of differences between a database of CIS patients with scans at two time points, to evaluate the ability of ODF-based scores to highlight early differences that may lead to the appearance or disappearance of lesions at the second time point.

We present in Section 2 our new methodology for the comparison of higher order models on the sphere, the databases on which our method was applied and our overall pipeline for longitudinal comparison. Then, we present in Section 3 our experiments, first based on simulated data with lesions where the ground truth of differences is known. This first section demonstrates the evaluation of our new framework compared to existing algorithms. Then, we focus on the results of the longitudinal study of MS patients. Finally, we discuss our results and conclude in Section 4.

Material and Methods

Patient to Group Diffusion MRI Comparison

The central part of our framework, illustrated later on with experimental data in Section 2.3, relies on the voxelwise comparison of a diffusion model to a group of diffusion models, e.g. the comparison of the DTI value of a patient vs the DTI values of a set of controls. In the following, we assume that all images (patient and controls) are aligned into a common space where each voxel of each image describes the exact same spatial position. Based on this assumption, we present a change detection method that applies both to diffusion tensor images and higher order models such as orientation distribution functions (ODF). We detail in the following the general comparison framework and its application to tensor and ODF images.

General Abnormalities Detection Framework

In [START_REF] Commowick | Detection of DTI white matter abnormalities in multiple sclerosis patients[END_REF], we introduced a methodology for the comparison of diffusion tensors between a patient and a group of controls. Here, we generalize and extend it to virtually any model. We wish to compare, at a given voxel, the patient model M P and models from the control subjects database M j , j = 1 . . . N . From each of these diffusion models, let us assume that it can be represented as a K-dimensional vector of features V (V P or V j ). We could then apply any comparison framework to those vectors. However, artifactual detections may appear due to noise in the images or registration errors. This is particularly true for complex diffusion models whose estimation from diffusion MRI may be more sensitive to these errors. The comparison is therefore preceded by a Principal Component Analysis (PCA) as follows:

• Perform a PCA from the control database to retain the h most meaningful eigenvectors

• Project V P and the V j onto a reduced h-dimensional space utilizing the eigenvector basis

• Compute a comparison score and p-value

To perform the PCA analysis, the covariance matrix Σ C is derived from the controls database. The eigenvalues and eigenvectors of Σ C are then computed and we retain only the h largest eigenvalues and eigenvectors. These vectors are in turn utilized to project the input data (patient and controls) on an h-dimensional space (h K). After projection, we therefore seek to compare a vector V P,h and a set of vectors V j,h from the controls database. The controls population V = {V j,h , j = 1 . . . N } is assumed to follow a multivariate Normal distribution N ( V , Σ V ), where V and Σ V denote respectively the average and covariance matrix of the population V. Under this assumption, we compute the difference statistic between V P,h and N ( V , Σ V ) as a Mahalanobis distance:

d 2 (V P,h ) = V P,h -V T Σ -1 V V P,h -V (1)
d 2 varies between 0 and infinity, getting larger as the patient vector gets less likely to belong to the multivariate Normal distribution of the controls. A pvalue can then be computed from this distance by recalling that the statistic

T = N (N -h) h(N 2 -1)
d 2 follows a Fisher distribution with parameters h and N -h: T ∼ F (h, N -h). The difference test p-value is therefore written as:

p(V P,h ) = 1 -F h,N -h d 2 (V P,h ) (2)
where F h,N -h is the cumulative distribution function of a Fisher distribution with parameters h and N -h.

DTI Abnormalities Detection

In the particular case of DTI, we wish to compare tensor values. Tensors lie in the space of 3 × 3 positive definite matrices and some care should be taken so that computations account for this specific space. We utilize the log-Euclidean framework proposed by Arsigny et al. (2006b) to perform the comparison. Arsigny et al. have indeed defined a Lie group for tensors (relying on their matrix principal logarithm), and Euclidean operations can therefore be computed on their logarithms. For DTI, we therefore chose to compare the vector representations of the logarithms of tensor values Vec [log(D)]where D is a tensor, log is the matrix logarithm and Vec the vector operator as proposed by Arsigny et al. (2006b). PCA and vector comparisons are then computed on these values, and in this case h may vary in between 0 and 6.

ODF Abnormalities Detection

Higher order models may allow us to capture more differences as they are better able to capture regions where fibers are crossing. In the following experiments, we will study Orientation Distribution Functions (ODF) represented on a modified Spherical Harmonics (SH) basis as proposed by [START_REF] Descoteaux | Regularized, Fast, and Robust Analytical Q-Ball Imaging[END_REF]:

Y j (θ, φ) =    √ 2 2 (-1) m Y -m l (θ, φ) + Y m l (θ, φ) , if m < 0 Y 0 l , if m = 0 i √ 2 2 (-1) m Y -m l (θ, φ) -Y m l (θ, φ) , if m > 0 (3)
where the couple (θ, φ) denotes a point on the sphere, l = 0, 2, 4, . . . , L is the order of the SH basis, m = -l, . . . , 0, . . . , l is an index, j := j(l, m) = (l 2 + l + 1) + m. Y j is the j-th element of the SH basis, and Y m l is the complex SH value. An ODF M expressed in this basis is therefore a vector of

L 2 + L + 1 parameters c j : M (θ, φ) = j c j Y j (θ, φ).
From each of the diffusion models, we sample the probability distribution function (pdf) of diffusion in K directions regularly placed over the unit sphere thereby constructing for each model a K-dimensional vector V (V P or V j ). The comparison of the patient with respect to the control database is then performed using Eqs (1) and ( 2), with the projected data from PCA. Interestingly, following such an approach allows us to be independent of the diffusion model used: if a probability density function is available for the model, then the presented statistical test may be used.

Diffusion MRI Databases 2.2.1. Simulated Data

To quantitatively evaluate the detection performance of our "one patient vs control group" comparison framework, we have created two simulated databases where lesions are known. The first database consists of crossing fibers with and without lesions (illustrated in Fig. 1). A noise-free DWI image of two fibers crossing at 90 degrees was first constructed using a multitensor model, the relative weights of each fiber being equal in the crossing region. The simulated diffusion weighted image consists of one B0 and 81 directions with a b-value of 1000 s.mm -2 (image size of 128×128). Then, a simulated lesion was added inside the crossing region, so that each compartment ADC (respectively FA) was increased (respectively decreased). Then, Rician noise was added repeatedly to these two images so that a database of 50 simulated controls and 100 simulated cases was created. Different noise levels were tested to study its influence on detection results: low level (3% of the image peak value), mid level (5% of the image peak value) and high level of noise (7% of the image peak value).

The second dataset (illustrated in Fig. 2) was created from a real control subject acquisition, for whom a diffusion sequence with 30 directions (image size: 128×128×60, b-value: 1000 s.mm -2 , 2×2×2 mm 3 voxel resolution) was acquired. We first added a lesion in the control DWI volume, located in a known region of crossing fibers: fibers going through the corpus callosum and corticospinal tracts. The diffusion properties of the lesion were then computed by comparing values of DWI signals (for each gradient direction) in known patient lesions and corresponding DWI signals in controls at the same positions (after registration). For each DWI image, we therefore computed a ratio indicating the relative change between a lesion signal and a normal signal. This intensity ratio was then applied to the selected control subject in the lesion region delineated. Rician noise was also added repeatedly to the two images so that a database of 50 simulated cases and 50 simulated controls were obtained. As for the first dataset, three different levels of noise were also tested (3%, 5% and 7% of the image peak value). While the first database aims at proving the interest of higher order models for regions of crossing fibers when HARDI sequences are acquired, the second database will be used to evaluate the proposed approach in a realistic frame where lower angular resolution images are acquired.

Multiple Sclerosis CIS Patients

We have used a database of 15 patients acquired from their first Clinically Isolated Syndrome (CIS), suggestive of multiple sclerosis. Among those 15 patients, 53.3 % were women and the average age was 29.8 ± 8.5 years old. MRI data was acquired using the same protocol at two time points: right after the first CIS event, hereafter called M 0 , and 3 months later, hereafter called M 3 . At each of these dates, T1w, T2w, PDw and FLAIR images were acquired on a Siemens Verio 3T scanner, as well as diffusion weighted MRI on a single shell. The T1w sequence was a 3D-MPRAGE with a 1 mm 3 isotropic voxel resolution, on a 256×256×160 matrix. T2w and PDw were acquired using a dual Turbo Spin Echo sequence with a 1×1×3 mm 3 voxel resolution, on a 192×256×44 matrix. The diffusion weighted MRI acquisition was a 30 directions acquisition on a single shell of b-value 1000 s.mm -2 , with a 2×2×2 mm 3 voxel resolution, on a 128×128×60 matrix.

A trained radiologist then delineated all lesions at all time points for those patients, based on T2w and FLAIR images which are conventionally used in clinical routine for lesion delineation.

Healthy Controls

In addition to the patients, a database of 46 matched control subjects was acquired. These subjects were selected so that their age and gender matched those from the patients (60.9 % of the controls were women and the average age was 29.3 ± 6.9 years old). The same protocol was performed on the same scanner for those subjects but at only one time point since no evolution was expected. T1w, T2w, PDw and FLAIR images were acquired on a Siemens Verio 3T scanner. The T1w sequence was a 3D-MPRAGE with a 1 mm 3 isotropic voxel resolution, on a 256×256×160 matrix. T2w and PDw were acquired using a dual Turbo Spin Echo sequence with a 1×1×3 mm 3 voxel resolution, on a 192×256×44 matrix. A diffusion weighted MRI acquisition was added with 30 directions on a single shell of b-value 1000 s.mm -2 , with a 2×2×2 mm 3 voxel resolution, on a 128×128×60 matrix.

Multiple Sclerosis Processing Pipeline

The processing pipeline we developed is divided into two parts. First, the input images were pre-processed as illustrated in Fig. 3. All gradient volumes of each DWI image was first registered onto the corresponding B0 image to correct for motion during the acquisition. To do so, a global rigid transform was computed using block-matching registration as proposed by [START_REF] Ourselin | Block matching: A general framework to improve robustness of rigid registration of medical images[END_REF]. Then, EPI distortion correction was applied by non-linearly registering each gradient image to the B0, utilizing the FFD algorithm from [START_REF] Rueckert | Nonrigid registration using free-form deformations: Application to breast MR images[END_REF] with few control points (5 in each direction). Then, tensors were estimated as well as ODF using [START_REF] Descoteaux | Regularized, Fast, and Robust Analytical Q-Ball Imaging[END_REF] estimation method. To conclude this pre-processing pipeline, T1w images and diffusion images were registered together up to a global rigid transformation to get an anatomical reference for our results. Once this preprocessing was performed for all input images (controls and patients at M 0 and M 3 ), we then utilized the pipeline illustrated in Fig. 4 to compare, voxel by voxel, each patient to the database of controls. First, a DTI atlas was computed from the controls images, utilizing an updated method derived from [START_REF] Guimond | Average brain models: A convergence study[END_REF] to handle transformations in the log-Euclidean framework on diffeomorphisms proposed by Arsigny et al. (2006a). Non linear registration for atlas construction was performed based on full tensors utilizing a block-matching strategy introduced by [START_REF] Suarez | Automated delineation of white matter fiber tracts with a multiple region-ofinterest approach[END_REF]. This registration proceeds by iteratively computing matches between blocks in the tensor images to estimate dense field corrections. We used the following parameters for the registration: three pyramid levels, block size 5 × 5 × 5 voxels, block neighborhood of 2 voxels in each direction. The obtained transformations were then applied on T1w images, and on ODF images, using the local rigid rotations extracted from the transformation fields [START_REF] Blanco | Evaluation of the rotation matrices in the basis of real spherical harmonics[END_REF]; [START_REF] Ritchie | Fast computation, rotation and comparison of low resolution spherical harmonic molecular surfaces[END_REF]; Geng et al. ( 2009)). The first time point (M 0 ) of each patient was then registered onto the average atlas image utilizing the same registration methods with the same parameters. Then, for subsequent time points, only a supplementary global rigid registration to M 0 is necessary to bring the DTI and ODF image onto the atlas, as we can assume that the global brain shape of the patient did not change within three months. Once M 0 and M 3 time points are aligned onto the atlas, we applied the framework proposed above to compare the ODFs extracted from the diffusion weighted images and report the results in the next section.

Experiments and Results

Diffusion Models Comparison on Simulated Data

Based on the simulated datasets presented in Section 2.2, we have estimated tensors and ODFs in a spherical harmonics basis of order 4. We have then performed an abnormality detection for each of the simulated cases independently, each time comparing the detection results with the known ground truth. Such an approach was applied on the two databases and results are reported in the following sub-sections.

Results on Simulated Crossing Fibers with Lesions

Fig. 5 presents a representative example of the FDR-corrected detection results at p < 0.05 for the mid-level of noise, with our approach (when keeping four eigenvectors from PCA analysis) based either on tensors or on ODFs. In addition, we tested the validity of the proposed statistical tests on this dataset: we calculated the the ratio of false positive detections in the background before FDR correction. The numbers are close to 5% (respectively 5.39% and 5.27%), confirming the validity of the proposed tests. From Fig. 5, it appears that both methods are well capable to detect the true lesion. However, tensor-based detection seems more prone to false positives than ODF-based detection.

In addition to visual inspection and individual tests, we have then computed Dice scores between the detections obtained by both algorithms and the known ground truth of lesion. We report in Fig. 6 the mean and standard deviation of these scores over the 100 simulated cases and for the three different levels of noise. As the ODF-based method is dependent on the number of eigenvectors kept from PCA, we are also reporting the results depending on this number for ODF-based detection. This figure illustrates a clear difference between the ODF and tensorbased method, with an advantage to ODF for all noise levels. This demonstrates that, due to the inability of DTI to estimate water diffusion in crossing fibers regions, DTI-based comparison may result in degraded detection scores. On the contrary, ODF, being able to model multiple fiber directions in a voxel, is more able to capture differences between a patient and controls. ODF-based detection achieves very good results for low and mid-levels of noise, with only a small decrease for 14 and 15 eigenvalues, which comes from the fact that adding more eigenvectors is only bringing more noise at that stage. For a high level of noise, all detection results get lower scores but ODF detection performs still better than tensor-based comparison and gets high scores (above 0.8), although it is more impacted by noise when reaching a high number of kept eigenvectors, compared to other levels of noise.

Results on Simulated Images of Patients

Fig. 7 presents FDR-corrected detection results at p < 0.05 obtained on the second database of simulated data with a mid-level of noise (on real images) with our proposed approach based on either tensors or ODFs (when keeping three or six eigenvectors from PCA analysis). Even with lower angular resolution, both methods are well capable to detect the true lesion. Tensor-based detection is more prone to false positives than ODF-based detection. In addition, it is clear that the number of eigenvalues kept is important: keeping only three eigenvalues leads to a clear under-estimation of the abnormality both for tensor and ODF-based detection. We have then computed Dice scores between the detections obtained for each algorithm and the known ground truth of lesion. We report in Fig. 8 the mean and standard deviation of these scores over the 50 simulated cases for each noise level. This figure illustrates again the advantage of ODF-based detection, which outperforms tensor-based detection in almost all cases. The ODF-based detection best result is reached for 6 eigenvalues (Dice score respectively of 0.98, 0.97 and 0.86 for low, mid and high levels of noise, compared to tensor-based scores of 0.82, 0.73 and 0.56). Interestingly, this number of eigenvalues where the maximum is reached differs from previous experiments, indicating that either the change to a real image or to a lower angular resolution is affecting the optimal number of eigenvalues kept after PCA. Overall, although scores are lower with respect to the previous experiment, this figure confirms that ODF performs better than tensor-based comparison even for low angular resolution images.

Multiple Sclerosis Patients Analysis

We have then applied our processing pipeline to our database of multiple sclerosis patients, running it on each patient and each time point to obtain difference scores as expressed in Eq. 1. In addition, we have derived p-values from these scores and detect abnormalities as those voxels having a p-value p < 0.05 after FDR correction. We present in Figs. 9 and 10 results obtained for two representative patients out of the 15 patients studied. Appearing lesions at M 3 are highlighted in these figures with arrows. The ODF scores at M 0 show a higher score when a lesion will appear at M 3 (see arrows).

Both figures clearly show detections of abnormalities inside the regions of lesions in diffusion weighted MRI. In addition, it may be seen on Fig. 9 that some regions of the brain that appear as normal on both time points on T2w images are showing abnormalities in diffusion weighted imaging (right part of Fig. 9.c and 9.f). Comparing the results at time points M 0 and M 3 on regions of appearing lesions at M 3 (arrows on images b and e on Figs. 9 and 10), it seems that for both patients the ODF derived score is higher at M 0 when compared to the value in the normal appearing white matter. These score differences do not seem to be sufficient to have a definite large abnormality detection at M 0 which may be explained by the fact that the lesion has not yet developed. However, the scores seem to be a good indicator of lesions appearing three months later.

To further evaluate ODF scores as a marker of future lesion apparition, we have performed a quantitative study based on the ODF scores obtained Although the ODF scores at M 0 do not always highlight a significant difference from controls, these scores are higher when a lesion will appear at M 3 (see arrows). Legend: vertical axis: Mahalanobis distance.

for each patient. Thanks to the manual delineations at both M 0 and M 3 time points, we have computed the average for all patients of ODF-based scores in four regions: normal appearing white matter at both time points, regions of appearing, staying and disappearing lesions at M 3 . We report the Box and Whiskers plot of the ODF scores for all 15 patients in Fig. 11. These quantitative comparisons bring up several important results. First, normal appearing white matter (NAWM) does not significantly change of score between M 0 and M 3 (see NAWM in Fig. 11, paired t-test, p = 0.472). As expected, appearing and disappearing lesions have significantly different scores in between M 0 and M 3 (see appearing and disappearing in Fig. 11, paired t-tests, p = 0.017 for appearing lesions and p = 0.009 for disappearing lesions). Interestingly, lesions that were already present at M 0 and that stay at M 3 also have significantly different scores (see Staying in Fig. 11, paired t-test, p = 0.021), which highlights the fact that, although these lesions stay, their internal structure may also change in time.

Another important result is that the NAWM score at M 0 is significantly lower than the score at M 0 for those tissues where lesions will appear at M 3 (NAWM vs appearing lesions taken at M 0 , one-way ANOVA, p < 0.001). These results confirm the visual results on ODF score maps, presented in Figs. 9 and 10, that showed different scores in regions of future appearing lesions. In addition, staying lesions have a significantly higher score at M 0 than lesions that will disappear at M 3 (staying vs disappearing taken at M 0 , one-way ANOVA, p = 0.0015). Both these results clearly show that ODFbased scores are able to show in advance some regions where lesions will appear or disappear.

Conclusion

We have presented a new method for the individual comparison of a patient diffusion MRI with respect to a population of healthy subjects. It is based on a generic framework to compare high dimensional models of diffusion such as Orientation Distribution Functions expressed in a modified spherical harmonics basis. We have demonstrated how such models allow for a better detection of differences with respect to simpler diffusion models as the classical diffusion tensor on simulated datasets.

We have applied this framework for the detection of white matter structure abnormalities at a very early stage of multiple sclerosis. This study performed on two time points right after the first CIS event, separated by three months, showed significant differences inside the lesions but also significant differences in ODF-related scores in regions where lesions were set to appear or disappear at M 3 . This is a very significant result illustrating the potential prediction power of diffusion imaging for multiple sclerosis evolution, at the individual level, and could therefore lead to adapted treatments in the future based on single time points evaluation of the lesions.

Future works will further study this predictive power, namely by integrating evaluation of false positive detections among the results. Furthermore, other modalities may be of great interest in addition to diffusion imaging to characterize white matter focal and diffuse lesions, such as magnetization transfer imaging or relaxometry. We will further investigate how these modalities are providing complementary results and how to integrate them into one framework for the detection of abnormalities both at the patient and population level. In addition, these modalities could provide a further insight in white matter diffusion abnormalities, thereby helping to define a white matter degeneracy ground truth for the validation of our detection technique. Although not developed for MS lesion segmentation, our algorithm could also be coupled to other conventional modalities to help traditional MS lesion segmentation algorithms.

At a more technical level, we will also investigate more time points, and specifically designed ODF metrics such as the ones proposed by [START_REF] Goh | A nonparametric Riemannian framework for processing high angular resolution diffusion images and its applications to ODF-based morphometry[END_REF]. The influence of registration and normal variability among controls on the detection power could also be further studied, for example by adding typical transformations to each simulated database image, or by varying the simulated DWI intensities. Such an extensive comparison would bring further insight into our approach. Our comparison framework is also generic and could be applied to multi-compartment models such as multi-tensor [START_REF] Scherrer | Parametric Representation of Multiple White Matter Fascicles from Cube and Sphere Diffusion MRI[END_REF]) or DDI [START_REF] Stamm | A new multi-fiber model for low angular resolution diffusion MRI[END_REF]), which may be more adapted to low angular acquisitions and do not consider only the directional component of the diffusion.

Figure 1 :

 1 Figure 1: Illustration of Simulated Data from the First Database. Tensor maps derived from simulated diffusion weighted images. Legend: Tensor map without lesions (a,b) and with a simulated lesion in the crossing region (c,d). (b,d) are close-ups on (a,c).

Figure 2 :

 2 Figure 2: Illustration of Simulated Data from the Second Database. Tensor maps derived from simulated diffusion weighted images overlaid on top of average diffusion coefficient images. Legend: Tensor map without lesions (a,b) and with a simulated lesion in the crossing region (c,d). (b,d) are close-ups on the lesion region (a,c).

Figure 3 :

 3 Figure 3: Preprocessing Pipeline Illustration.

Figure 4 :

 4 Figure 4: Processing Pipeline Illustration.

Figure 5 :

 5 Figure 5: Lesion Difference Evaluation on Simulated Fiber Crossing Data. Zscore maps (top row) and thresholded p-value maps (bottom row) illustrating the differences of detection between tensor and ODF-based abnormality detection. Legend: difference score maps for tensor-based (a) and ODF-based (b) methods, (c): color range of values, (d): abnormality ground truth, detection maps (p < 0.05 after FDR correction) for tensor-based (e) and ODF-based (f) methods.

Figure 6 :

 6 Figure 6: Lesion Detection Dice Scores on Simulated Fiber Crossing Data. Dice overlap scores for lesion detection on simulated crossing fibers with known lesion ground truth, for tensor and ODF-based detection , on three different levels of noise. Results for tensor-based detection above six eigenvalues are set to the results with six eigenvalues for comparison purposes.

Figure 7 :

 7 Figure 7: Abnormality Detection on Simulated Patient Data. Abnormalities detected at a p < 0.05 level after FDR correction with the evaluated methods. Patient ADC image (a), patient ADC image overlayed with (b): true simulated lesion, tensor-based detections keeping the 3 (c) or 6 (d) eigenvectors in the PCA step, ODF-based detections keeping the 3 (e) or 6 (f) eigenvectors.

Figure 8 :

 8 Figure 8: Lesion Detection Dice Scores on Simulated Patient Data. Dice overlap scores for lesion detection on a simulated patient with known ground truth, on three different levels of noise. Results for tensor-based detection above six eigenvalues are set to the results with six eigenvalues for comparison purposes.

Figure 9 :

 9 Figure 9: ODF-Based Difference Scores on one MS Patient with Appearing Lesions. Illustration of the ODF scores at the two time points for one patient with appearing lesions. First row: at M 0 , second row: at M 3 . (a,d) are axial views of T2-w images overlaid with the manual lesion delineation from the expert. (b,e) correspond to the same axial views of the ODF scores as computed from our pipeline. (c,f) correspond to the abnormalities with respect to normal anatomy detected (at p < 0.05 after correction).The ODF scores at M 0 show a higher score when a lesion will appear at M 3 (see arrows).

Figure 10 :

 10 Figure 10: ODF-Based Difference Scores on a Second MS Patient with Appearing Lesions. Illustration of the ODF scores at the two time points for a second patient with appearing lesions. First row: at M 0 , second row: at M 3 . (a,d) are coronal views of T2w images overlaid with the manual lesion delineation from the expert. (b,e) correspond to the same views of the ODF scores as computed from our pipeline. (c,f) correspond to the abnormalities with respect to normal anatomy detected (at p < 0.05 after correction).Although the ODF scores at M 0 do not always highlight a significant difference from controls, these scores are higher when a lesion will appear at M 3 (see arrows). Legend: vertical axis: Mahalanobis distance.

Figure 11 :

 11 Figure 11: ODF-Based Difference Scores on 15 MS Patients. Illustration of the average difference scores computed from ODFs over the 15 MS patients, respectively on NAWM, appearing, staying and disappearing lesions. Scores in lesions appearing at M 3 are significantly different at M 0 from normal appearing white matter.

Acknowledgments

The authors acknowledge Neurinfo MRI research facility for logistic support and MRI data acquisition. This work was supported by ARSEP, France (Association pour la Recherche de la Sclérose en Plaques).