Mireille Bousquet-M Élou

THE EXPECTED NUMBER OF INVERSIONS AFTER n ADJACENT TRANSPOSITIONS

We give a new expression for the expected number of inversions in the product of n random adjacent transpositions in the symmetric group S m+1 . We then derive from this expression the asymptotic behaviour of this number when n ≡ nm scales with m in various ways. Our starting point is an equivalence, due to Eriksson et al., with a problem of weighted walks confined to a triangular area of the plane.

Introduction

Let S m+1 denote the group of permutations on the set {0, 1, 2, . . . , m}. Define the adjacent transposition s i to be the two-cycle (i, i + 1). We consider the adjacent transposition Markov chain π (0) , π (1) , π (2) , . . . on S m+1 . That is, π (0) is the identity, and π (n+1) = π (n) s i , where i ∈ {0, . . . , m -1} is chosen uniformly at random. For instance, denoting π = π 0 π 1 • • • π m if π(i) = π i for all i, we may have, when m = 3, π (0) = 0123 = id, π (1) = 0213 = π (0) s 1 , π (2) = 2013 = π (1) s 0 , π (3) = 2103 = π (2) s 1 .

Markov chains on finite groups have attracted a lot of interest in, at least, the past 30 years, at the interface of algebra and probability theory [START_REF] Aldous | Random walks on finite groups and rapidly mixing Markov chains[END_REF][START_REF] Diaconis | Group representations in probability and statistics[END_REF][START_REF] Diaconis | Comparison theorems for reversible Markov chains[END_REF][START_REF] Diaconis | Generating a random permutation with random transpositions[END_REF][START_REF] Saloff-Coste | Probability on groups: Random walks and invariant diffusions[END_REF][START_REF] Saloff-Coste | Random walks on finite groups[END_REF]. Such chains are also studied in computational biology, in connection with models of gene mutations [START_REF] Berestycki | A phase transition in the random transposition random walk[END_REF][START_REF] Durrett | Probability models for DNA sequence evolution[END_REF][START_REF] Eriksen | Estimating the expected reversal distance after a fixed number of reversals[END_REF][START_REF] Eriksson | Expected number of inversions after a sequence of random adjacent transpositions[END_REF][START_REF] Wang | Estimating true evolutionary distances between genomes[END_REF]. A central question is to estimate the mixing time of the chain, for large values of m. (The mixing time is the number of steps the chain takes to approach its equilibrium distribution.) The above chain is periodic of period 2, since π (2n) is always in the alternating group. Thus it has no equilibrium distribution. However, an elementary modification (which consists in choosing π(n+1) = π(n) with probability 1/(m + 1) and otherwise multiplying by a transposition s i chosen uniformly) makes it aperiodic, and the equilibrium distribution is then uniform over S m+1 . The mixing time is known to be Θ(m 3 log m) (see Aldous [START_REF] Aldous | Random walks on finite groups and rapidly mixing Markov chains[END_REF], Diaconis and Saloff-Coste [START_REF] Diaconis | Comparison theorems for reversible Markov chains[END_REF] and Wilson [START_REF] Wilson | Mixing times of lozenge tiling and card shuffling Markov chains[END_REF]).

More recently, Berestycki and Durrett [START_REF] Berestycki | Limiting behavior of the distance of a random walk[END_REF] studied a continuous time version of this chain, denoted (X m t ) t≥0 , where multiplications by adjacent transpositions occur according to a Poisson process of rate 1. The connection with the chain π (n) described above is straightforward. The authors focussed on the inversion number D m t of X m t , as a function of m and t. They established the convergence in probability of this random variable, suitably normalized, when t ≡ t m scales with m in various ways. The limit is usually described in probabilistic terms, except in one simple regime, m ≪ t ≪ m 3 , where D m t / √ mt is shown to converge to 2/π. Of course, the interesting regimes are those occurring before stationarity.

This paper actually stems from a remark made in Berestycki & Durrett's paper. They quote a paper by Eriksen [START_REF] Eriksen | Expected number of inversions after a sequence of random adjacent transpositions -an exact expression[END_REF], which gives a rather formidable expression for I m,n , the expected number of inversions in π (n) : .

I m,n = n r=1 1 m r n r r s=1 r -1 s -1 (-4) r-s g s,m h s,m , (1)
The authors underline that "it is far from obvious how to extract useful asymptotic from this formula". This observation motivated our interest in the expectation I m,n . As the problem is of an algebraic nature, it should have much structure: would it be possible to find an alternative expression of I m,n , with a neater structure than [START_REF] Aldous | Random walks on finite groups and rapidly mixing Markov chains[END_REF], and to derive asymptotic results from this expression? This is precisely what we do in this paper. Our alternative formula for I m,n reads as follows.

Theorem 1. Let m ≥ 1. For k ∈ Z, denote

α k = (2k + 1)π 2m + 2 , c k = cos α k , s k = sin α k and x jk = 1 - 4 m (1 -c j c k ). ( 2
)
The expected number of inversions after n adjacent transpositions in S m+1 is

I m,n = m(m + 1) 4 - 1 8(m + 1) 2 m k,j=0 (c j + c k ) 2 s 2 j s 2 k x jk n .
Equivalently, the generating function

I m (t) = n≥0 I m,n t n is I m (t) = m(m + 1) 4(1 -t) - 1 8(m + 1) 2 m k,j=0 (c j + c k ) 2 s 2 j s 2 k 1 1 -tx jk . ( 3
)
This result is proved in Section 2, and asymptotic results are derived in Sections 3 to 5 for three main regimes: linear (n m = Θ(m)) and before, cubic (n m = Θ(m 3 )) and beyond, intermediate (m ≪ n m ≪ m 3 ). For the moment, let us give a few comments and variants on Theorem 1.

Limit behaviour. The chain π (n) has period 2, as π (n) is an even (resp. odd) permutation if n is even (resp. odd). But the sub-chains π (2n) and π (2n+1) are aperiodic on their respective state spaces, the alternating group A m+1 (the group of even permutations) for π (2n) and its complement S m+1 \ A m+1 for π (2n+1) . Moreover, each of these chains is irreducible and symmetric, and thus admits the uniform distribution as its equilibrium distribution. For m ≥ 3, the average number of inversions of an element of A m+1 is m(m + 1)/4, and the same is true of elements of S m+1 \ A m+1 . Thus when m ≥ 3 is fixed and n → ∞, we expect I m,n to tend to m(m + 1)/4. This can be seen on the above expression of I m,n , upon observing that, for m ≥ 3 and 0 ≤ j, k ≤ m with j + k = m, there holds x jk ∈ (-1, 1). (The condition j + k = m is equivalent to c j + c k = 0.) If m ≥ 8, the stronger property x jk ∈ (0, 1) holds, which shows that I m,n is an increasing function of n.

Eigenvalues of the transition matrix. Another consequence of the above theorem is that we have identified a quadratic number of eigenvalues of the transition matrix of the chain (see [START_REF] Eriksen | Expected number of inversions after a sequence of random adjacent transpositions -an exact expression[END_REF] for a detailed account of the connection between these eigenvalues and I m,n ).

Corollary 2. Let c k = cos (2k+1)π
2m+2 . The transition matrix of the adjacent transposition Markov chain on S m+1 admits -among others-the following eigenvalues:

x jk = 1 - 4 m (1 -c j c k ), for 0 ≤ j, k ≤ m and j + k = m.
This result is not really new: as pointed to us by David Wilson, one can derive from the proof of Lemma 9 in [START_REF] Wilson | Mixing times of lozenge tiling and card shuffling Markov chains[END_REF] that for 0 ≤ p, q ≤ m, with p = q,

1 - 2 m 2 -cos pπ m + 1 -cos qπ m + 1
is also an eigenvalue. This collection is larger than the one we have exhibited, as can be seen using

2c j c k = cos (k -j)π m + 1 + cos (k + j + 1)π m + 1 .
The transition matrix has still more eigenvalues, for instance -1, with eigenvector (ε(σ)) σ∈Sm+1 where ε denotes the signature. The paper [START_REF] Edelman | Codes, transforms and the spectrum of the symmetric group[END_REF], originally motivated by a coding problem, gives a linear number of eigenvalues of the form i/m, where i is an integer. There exists a description of all eigenvalues in terms of the characters of S m+1 (see [START_REF] Diaconis | Generating a random permutation with random transpositions[END_REF]Thm. 3]). This description is valid in a much more general framework, and it seems that the complete list of eigenvalues is not explicitly known.

Rationality of I m (t). That the series I m (t) is always rational should not be a surprise: this property is clear when considering the transition matrix of the chain. Here are the first few values: A related aperiodic chain. If we consider instead the aperiodic variant of the chain (obtained by choosing π(n+1) = π(n) with probability 1/(m + 1) and otherwise multiplying by a transposition s i chosen uniformly), it is easy to see that the expected number of inversions after n steps is now

I 1 (t) = t (1 -t)(1 + t) , I 2 (t) = t(2 + t) (1 -t)(2 -t)(1 + t) , I 3 (t) = 3t(27 + 9 t -7 t 2 -t 3 ) (1 -t)(9 + 6 t -t 2 )(9 -
Ĩm,n = n k=0 n k m k (m + 1) n I m,k ,
so that the associated generating function is also rational:

Ĩm (t) := n≥0 Ĩm,n t n = 1 1 -t/(m + 1) I m tm m + 1 -t .
More generally, if a transposition occurs with probability p, the generating function of the expected number of inversions is

1 1 -t(1 -p) I m tp 1 -t(1 -p) .
(Above, p = m/(m + 1), but in [START_REF] Wilson | Mixing times of lozenge tiling and card shuffling Markov chains[END_REF] for instance, p = 1/2.)

Alternative expressions. Theorem 1 gives a partial fraction expansion of I m (t). One of the advantages of (3) is that the coefficients involved in the sum over j and k are negative, which makes the list of poles of I m (t) clear (we have actually already used this to state Corollary 2).

A number of variations are possible, and we will use some of them in the asymptotic study of the numbers I m,n . For instance:

I m (t) = m(m + 1) 4(1 -t) - 1 8(m + 1) 2 m k,j=0 c j + c k (1 -c j )(1 -c k ) 1 1 -tx jk (4) = 1 8(m + 1) 2 m k,j=0 c j + c k (1 -c j )(1 -c k ) 1 1 -t - 1 1 -tx jk . (5)
Both formulas are proved in Section 2.

The expected number of inversions

We prove in this section Theorem 1 and its variants (4-5). Our starting point is a functional equation satisfied by a series related to I m (t). We then solve this equation using the bivariate kernel method, which has proved successful in the past few years in various enumerative problems related to lattice paths, permutations, and other combinatorial objects [START_REF] Bousquet-Mélou | Four classes of pattern-avoiding permutations under one roof: generating trees with two labels[END_REF][START_REF] Bousquet-Mélou | Walks in the quarter plane: Kreweras' algebraic model[END_REF][START_REF] Bousquet-Mélou | Walks with small steps in the quarter plane[END_REF][START_REF] Bousquet-Mélou | On partitions avoiding 3-crossings[END_REF][START_REF] Janse Van Rensburg | Partially directed paths in a wedge[END_REF][START_REF] Mishna | Two non-holonomic lattice walks in the quarter plane[END_REF][START_REF] Xin | Enumeration of bilaterally symmetric 3-noncrossing partitions[END_REF]. The equation under consideration here can be interpreted in terms of weighted lattice paths confined to a triangular portion of the square lattice. Another problem of walks in a triangle, also related to the adjacent transposition Markov chain, was studied by Wilson [START_REF] Wilson | Mixing times of lozenge tiling and card shuffling Markov chains[END_REF] via a diagonalization of the adjacency matrix. The enumeration of (unweighted) walks in a triangle was performed by Flajolet [START_REF] Ph | The evolution of two stacks in bounded space and random walks in a triangle[END_REF] using the reflection principle, in connection with a problem of storage raised by Knuth [24, Ex. 2.2.2.13].

A functional equation

For π ∈ S m+1 , let us write π = π 0 π 1 • • • π m if π(i) = π i for all i. For n ≥ 0 and 0 ≤ i ≤ j < m, let p (n) i,j = P(π (n) i > π (n) j+1
). Then the expected number of inversions in π (n) is

I m,n = 0≤i≤j<m p (n) i,j . (6)
Examining how the numbers p (•) i,j may change at the n th step gives a recurrence relation for these numbers, first obtained by Eriksson et al. [START_REF] Eriksson | Expected number of inversions after a sequence of random adjacent transpositions[END_REF]. As shown by the lemma below, it converts our problem into the study of weighted walks confined to a triangular region of the square lattice. Consider the subgraph G m of the square lattice Z × Z induced by the points (i, j), with 0 ≤ i ≤ j < m (Figure 1). We use the notation (i, j) ↔ (k, ℓ) to mean that the points (i, j) and (k, ℓ) are adjacent in this graph.

(m -1) 1 (m -1) 0 1 Figure 1. The graph G m .

Lemma 3 ([20]). The inversion probabilities p (n)

i,j are characterized by the following recursion:

p (0) i,j = 0 for 0 ≤ i ≤ j < m, and for n ≥ 0, p (n+1) i,j = p (n) i,j + 1 m (k,ℓ)↔(i,j) p (n) k,ℓ -p (n) i,j + δ i,j m 1 -2p (n) i,j
, where δ i,j = 1 if i = j and 0 otherwise.

As often, it is convenient to handle the numbers p

(n)
i,j via their generating function:

P (t; u, v) ≡ P (u, v) := n≥0 t n 0≤i≤j<m p (n) i,j u i v j .
Multiplying the above recursion by t n+1 , and then summing over n, gives the following functional equation for P (u, v).

Corollary 4. The series P (u, v) satisfies

1 -t + t m (4 -u -ū -v -v) P (u, v) = t m 1 -u m v m (1 -uv)(1 -t) -(ū -1)P ℓ (v) -(v -1)v m-1 P t (u) -(u + v)P d (uv) , (7)
where ū = 1/u, v = 1/v, and the series P ℓ , P t and P d describe the numbers p

i,j on the three borders (left, top, and diagonal) of the graph G m :

P ℓ (v) = n≥0 t n 0≤j<m p (n) 0,j v j , P t (u) = n≥0 t n 0≤i<m p (n) i,m-1 u i , P d (u) = n≥0 t n 0≤i<m p (n) i,i u i . (8)
In view of ( 6), the generating function we are interested in is

I m (t) = n≥0 I m,n t n = P (1, 1),
which, according to the functional equation [START_REF] Bousquet-Mélou | Four classes of pattern-avoiding permutations under one roof: generating trees with two labels[END_REF], may be rewritten

I m (t) = t (1 -t) 2 - 2tP d (1) m(1 -t) . (9)
In the next subsection, we solve [START_REF] Bousquet-Mélou | Four classes of pattern-avoiding permutations under one roof: generating trees with two labels[END_REF], at least to the point where we obtain a closed form expression of P d (1), and hence a closed form expression of I m (t), as announced in Theorem 1.

Solution of the functional equation

We first establish a symmetry property of the series P (u, v).

Lemma 5. The series P (u, v) satisfies

P (u, v) = u m-1 v m-1 P (v, ū) with ū = 1/u, v = 1/v.
In particular, the "diagonal" generating function P d (u) satisfies

P d (u) = u m-1 P d (ū). ( 10
)
Proof. This can be derived from the functional equation satisfied by P (u, v), but we prefer to give a combinatorial (or probabilistic) argument.

Let τ be the permutation of S m+1 that sends k to m -k for all k. Note that τ is an involution. Let Φ denote the conjugacy by τ : that is, Φ(σ) = τ στ . Of course, Φ(id) = id and Φ(σσ ′ ) = Φ(σ)Φ(σ ′ ). Also, Φ has a simple description in terms of the diagram of σ, in which σ(i) is plotted against i: the diagram of Φ(σ) is obtained by applying a rotation of 180 degrees to the diagram of σ. In particular, with s i = (i, i + 1), one has Φ(s i ) = s m-1-i for 0 ≤ i < m. These properties imply that the sequence of random permutations (Φ(π (0) ), Φ(π (1) ), Φ(π (2) ), . . .) follows the same law as the original Markov chain (π (0) , π (1) , π (2) , . . .). In particular, for 0

≤ i ≤ j < m, p (n) i,j = P(Φ(π (n) ) i > Φ(π (n) ) j+1 ) = P(π (n) m-j-1 > π (n) m-i ) = p (n) m-j-1,m-i-1 .
This is equivalent to the first statement of the lemma. The second follows by specialization.

The next ingredient in our solution is the "obstinate" kernel method of [START_REF] Bousquet-Mélou | Four classes of pattern-avoiding permutations under one roof: generating trees with two labels[END_REF][START_REF] Bousquet-Mélou | Walks in the quarter plane: Kreweras' algebraic model[END_REF][START_REF] Bousquet-Mélou | Walks with small steps in the quarter plane[END_REF][START_REF] Bousquet-Mélou | On partitions avoiding 3-crossings[END_REF][START_REF] Janse Van Rensburg | Partially directed paths in a wedge[END_REF][START_REF] Mishna | Two non-holonomic lattice walks in the quarter plane[END_REF][START_REF] Xin | Enumeration of bilaterally symmetric 3-noncrossing partitions[END_REF]. The kernel of the functional equation ( 7) is the coefficient of P (u, v), namely

K(u, v) = 1 -t + t m (4 -u -ū -v -v).
Let (U, V ) be a pair of Laurent series in t that cancel the kernel: K(U, V ) = 0. The series P (U, V ) is well-defined (because P (u, v) is a polynomial in u and v). Setting u = U and v = V in ( 7) cancels the left-hand side, and thus the right-hand side. That is, denoting as usually

Ū = 1/U and V = 1/V , V m-1 (1 -U m V m ) (1 -U V )( Ū -1)(V -1)(1 -t) - V m-1 P ℓ (V ) V -1 - P t (U ) Ū -1 - V m-1 (U + V ) ( Ū -1)(V -1) P d (U V ) = 0 provided U = 1 and V = 1.
Let us now exploit the symmetries of the kernel: obviously, K(u, v) is invariant by the transformations u → ū and v → v. Hence the pairs ( Ū , V ), ( Ū , V ) and (U, V ) also cancel K, and it follows that

V m-1 (1 -Ū m V m ) (1 -Ū V )(U -1)(V -1)(1 -t) - V m-1 P ℓ (V ) V -1 - P t ( Ū ) U -1 - V m-1 ( Ū + V ) (U -1)(V -1) P d ( Ū V ) = 0, V m-1 (1 -Ū m V m ) (1 -Ū V )(U -1)( V -1)(1 -t) - V m-1 P ℓ ( V ) V -1 - P t ( Ū ) U -1 - V m-1 ( Ū + V ) (U -1)( V -1) P d ( Ū V ) = 0, V m-1 (1 -U m V m ) (1 -U V )( Ū -1)( V -1)(1 -t) - V m-1 P ℓ ( V ) V -1 - P t (U ) Ū -1 - V m-1 (U + V ) ( Ū -1)( V -1) P d (U V ) = 0.
Let us form the alternating sum of the four previous equations: all occurrences of P ℓ and P t vanish, and, using [START_REF] Bousquet-Mélou | On partitions avoiding 3-crossings[END_REF], one obtains, after multiplying by

U m V m (1 -U )(1 -V ): (U V + 1) U m+1 + V m+1 P d (U V ) + V m-1 (U + V ) 1 + U m+1 V m+1 P d U V = U V 1 -t (1 -U m V m ) (U m + V m ) 1 -U V - (U m -V m ) (1 + U m V m ) V -U (11)
as soon as K(U, V ) = 0, U = 1 and V = 1. This equation involves only one unknown series, namely P d . So far, the series U and V are coupled by a single condition, K(U, V ) = 0. Let q be a complex root of q m+1 = -1, and let us add a second constraint on the pair (U, V ) by requiring that U = qV . That is, V must be a root of

K(qV, V ) = 1 -t + t m 4 -(1 + q)(V + q V ) = 0 ( 12
)
where q = 1/q and V = 1/V . We further assume q = -1 (otherwise K(qv, v) is independent of v). Then V = 1, U = qV = 1, and the first term of (11) vanishes. We thus obtain an explicit expression of P d (q), which we write below in terms of V and its conjugate root V ′ = q V . Lemma 6. Let q = -1 satisfy q m+1 = -1, and let V ≡ V (t) and V ′ ≡ V ′ (t) be the two roots of [START_REF] Diaconis | The cutoff phenomenon in finite Markov chains[END_REF]. Then V V ′ = q = 1/q and

P d (q) = q (1 -t)(V m+1 + V ′m+1 ) V m + V ′m 1 -q + 1 - q 1 + q V m -V ′m V -V ′ .
The series V and V ′ are algebraic (of degree 2) over C(t). But their symmetric functions are rational, and thus P d (q) is rational too, as expected. The following lemma gives an explicit rational expression of P d (q). Lemma 7. Let q = -1 satisfy q m+1 = -1. Assume q = q k := e i 2k+1 m+1 π . Then

P d (q k ) = 2it q k (1 -t)m(m + 1)s k m j=0 (c k + c j )(1 -c k c j ) 1 -tx jk with c j = cos (2j + 1)π 2m + 2 , s j = sin (2j + 1)π 2m + 2 and x jk = 1 - 4 m (1 -c j c k ).
Equivalently,

P d (q k ) = ic k 2q k (1 -t)s k - i 2q k (m + 1)s k m j=0 c k + c j 1 -tx jk .
Proof. We will establish a closed form expression for the coefficient of t n in (1 -t)P d (q), which is clearly equivalent to the first expression of P d (q) given above: for n ≥ 1,

a n := [t n ](1 -t)P d (q) = 2i qm(m + 1)s k m j=0 (c k + c j )(1 -c k c j ) x jk n-1 . (13)
In order to obtain this expression, we begin with applying Cauchy's formula to the expression of (1 -t)P d (q) given in Lemma 6. Let V be the root of ( 12) that vanishes at t = 0. (The other root V ′ = q V has a term in O(1/t) in its expansion.) Then

a n = 1 2iπ R(V ) dt t n+1
where the integral is taken along a small circle around the origin, in counterclockwise direction,

R(v) = q v m+1 + (qv) m+1 v m + (qv) m 1 -q + 1 - q 1 + q v m -(qv) m v -qv , and v = 1/v. Obviously, R(v) = R(qv). By (12), 1 t = 1 - 4 m + 1 + q m (V + q V ) so that dt t 2 = - 1 + q m (V -q V ) dV V .
Thus the integral expression of a n reads

a n = 1 2iπ S(v) dv v
where the integral is taken along a small circle around the origin, in counterclockwise direction, and

S(v) = - 1 + q m (v -qv) 1 - 4 m + 1 + q m (v + qv) n-1 R(v)
satisfies S(qv) = -S(v). Note that the only possible poles of S(v)/v are 0 and the (2m + 2)th roots of unity (because q m+1 = -1). Thus a n is simply the residue of S(v)/v at v = 0. Performing the change of variables v = q w, with w = 1/w, gives

a n = 1 2iπ S(w) dw w
where the integral is now taken along a large circle around the origin, in clockwise direction. This integral thus collects (up to a sign) all residues of S(v)/v. The residue formula thus gives:

2a n = 1 2iπ S(w) dw w + 1 2iπ S(v) dv v = - v 2m+2 =1 1 v Res v (S) = 1 2m + 2 v 2m+2 =1 P (v) (14)
where

P (w) = (1 -w 2m+2 )S(w) = 1 + q m w m+1 (w -q w) w m + (q w) m 1 -q + 1 - q 1 + q w m -(q w) m w -q w 1 - 4 m + 1 + q m (w + q w) n-1
.

Take v = e iα with α = ℓπ/(m + 1) and 0 ≤ ℓ < 2m + 2 and recall that q = e iθ with θ = (2k+1)π m+1 . Then 1 + q = 2e -iθ/2 cos(θ/2), 1 -q = -2ie iθ/2 sin(θ/2), 1q = 2ie -iθ/2 sin(θ/2), 1 + q = 2e iθ/2 cos(θ/2),

v m+1 = (-1) ℓ , v -qv = 2ie -iθ/2 sin(α + θ/2), v + qv = 2e -iθ/2 cos(α + θ/2), v m + (qv) m = 2i(-1) ℓ+1 e iθ/2 sin(α + θ/2), v m -(qv) m = 2(-1) ℓ e iθ/2 cos(α + θ/2).
Putting these identities together, one obtains

P (v) = 4i qm sin(θ/2) cos(θ/2) sin 2 (α + θ/2) + sin 2 (θ/2) cos(α + θ/2) × 1 - 4 m + 4 m cos(θ/2) cos(α + θ/2) n-1
, or, with the notation of the lemma,

P (v) = 4i qms k (c j + c k )(1 -c j c k ) x jk n-1
with j = k + ℓ. Returning to (14) now gives

a n = i qm(m + 1)s k m j=-m-1 (c k + c j )(1 -c k c j ) x jk n-1 ,
which is equivalent to [START_REF] Diaconis | Asymptotic analysis of a random walk on a hypercube with many dimensions[END_REF], upon noting that c j = c -j-1 . The first expression of P d (q) given in the lemma follows.

For the second expression, we simply perform a partial fraction expansion in the variable t, based on t

(1 -t)(1 -xt) = 1 1 -x 1 1 -t - 1 1 -xt , (15)
and use m j=0 c j = 0 (this identity follows for instance from c m-j = -c j ).

Recall that P d (u) is a polynomial in u of degree m -1. The above lemma gives its values at m, or even m + 1, distinct points. Thus P d (u) is completely determined by these values, and we can recover it by interpolation. We use a version of Lagrange's interpolation that is well-suited to symmetric polynomials.

Lemma 8. Let P (u) be a polynomial of degree m -1 with coefficients in some field containing

C. Assume P (u) is symmetric. That is, P (u) = u m-1 P (ū), with ū = 1/u. Let ℓ = ⌊ m-1
2 ⌋, and let q 0 , . . . , q ℓ be distinct elements of C such that q j q k = 1 for all k and j. Then

P (u) = (1 + u) χm,0 ℓ k=0 P (q k ) (1 + q k ) χm,0 j =k (u -q j )(u -1/q j ) (q k -q j )(q k -1/q j ) , ( 16
)
where χ m,0 = 1 if m is even, and 0 otherwise. When q k = e iθ k with θ k = (2k+1)π m+1 , this can be rewritten as

P (u) = - 2i m + 1 (u m+1 + 1) ℓ k=0 P (q k ) q k sin θ k u 2 -2u cos θ k + 1 .
In particular, with the notation (2),

P (1) = - 2i m + 1 ℓ k=0 P (q k ) q k c k s k .
Proof. For the first part, it suffices to observe that the expression given on the right-hand side of ( 16) has degree m -1, and takes the same values as P (u) at the 2ℓ + 2 distinct points q 0 , . . . , q ℓ , 1/q 0 , . . . , 1/q ℓ , and also at -1 if m is even. This gives a total of m + 1 correct values, which is more than enough to determine a polynomial of degree m -1.

For the second part, we observe that q 0 , . . . , q ℓ , 1/q 0 , . . . , 1/q ℓ , together with -1 if m is even, are the m + 1 roots of u m+1 + 1. Hence

(1 + u) χm,0 j =k (u -q j )(u -1/q j ) = u m+1 + 1 (u -q k )(u -1/q k ) ,
which, in the limit u → q k , gives

(1 + q k ) χm,0 j =k (q k -q j )(q k -1/q j ) = - m + 1 q k (q k -1/q k )
.

The second result follows. The third one is obtained by setting u = 1, and noticing that sin

θ k = 2c k s k , while 1 -cos θ k = 2s 2 k .
We can finally combine this interpolation formula with Lemma 7 to obtain an explicit expression of P d (u). As we are mostly interested in P d (1) (see ( 9)), we give only this value.

Proposition 9. Let ℓ = ⌊ m-1
2 ⌋, and adopt the notation [START_REF] Bayer | Trailing the dovetail shuffle to its lair[END_REF]. Then the series P d (u) defined by (8) satisfies:

P d (1) = m 2(1 -t) - 1 (m + 1) 2 ℓ k=0 m j=0 c k s 2 k c j + c k 1 -tx jk = m 2(1 -t) - 1 4(m + 1) 2 m k=0 m j=0 (c j + c k ) 2 (1 -c j c k ) s 2 j s 2 k 1 1 -tx jk .
Proof. We apply Lemma 8 to the second expression of P d (q k ) obtained in Lemma 7. This gives the first expression of P d (1) above, provided

ℓ k=0 c 2 k s 2 k = m(m + 1) 2 . ( 17
)
The latter identity is obtained by applying Lemma 8 to

P (u) = 1 + u + • • • + u m-1 = (1 - u m )/(1 -u).
We now seek a symmetric formula in j and k. Using

c m-j = -c j , c m-k = -c k , s m-k = s k and x j,k = x m-j,m-k , we write ℓ k=0 m j=0 c k s 2 k c j + c k 1 -tx jk = ℓ k=0 m j=0 c m-k s 2 m-k c m-j + c m-k 1 -tx m-j,m-k = m k=m-ℓ m j=0 c k s 2 k c j + c k 1 -tx jk = 1 2 m k=0 m j=0 c k s 2 k c j + c k 1 -tx jk = 1 4 m k=0 m j=0 c k s 2 k + c j s 2 j c j + c k 1 -tx jk
and this gives the second expression of P d (1).

Proof of Theorem 1. Let us now return to the generating function I m (t) whose coefficients give the expected number of inversions. It is related to P d (1) by ( 9). Theorem 1 is obtained by combining the second expression of Proposition 9, a partial fraction expansion in t (based on ( 15)), and finally the identity m j,k=0

(c j + c k ) 2 s 2 j s 2 k = 2m(m + 1) 3 . ( 18
)
To prove this identity, we write m j,k=0

(c j + c k ) 2 s 2 j s 2 k = m j,k=0 c 2 j + c 2 k s 2 j s 2 k (as c m-j = -c j and s m-j = s j ) = 2   m j=0 c 2 j s 2 j   m k=0 1 s 2 k = 2   m j=0 c 2 j s 2 j   m k=0 1 + c 2 k s 2 k
and complete the proof thanks to [START_REF] Edelman | Codes, transforms and the spectrum of the symmetric group[END_REF].

To obtain the expression (4) of I m (t), we write m j,k=0

(c j + c k ) 2 s 2 j s 2 k 1 1 -tx jk = 1 2 m j,k=0 c j + c k (1 -c j )(1 -c k ) - c j + c k (1 + c j )(1 + c k ) 1 1 -tx jk = m j,k=0 c j + c k (1 -c j )(1 -c k ) 1 1 -tx jk .
The latter identity follows from replacing j by m -j and k by m -k.

Let us now extract the coefficient of t 0 in (4). This gives

0 = m(m + 1) 4 - 1 8(m + 1) 2 m j,k=0 c j + c k (1 -c j )(1 -c k ) , ( 19
)
and the expression (5) of I m (t) follows.

Small times: linear and before

When m is fixed and n → ∞, the asymptotic behaviour of the numbers I m,n is easily derived from Theorem 1, as sketched just after the statement of this theorem. For m ≥ 3,

I m,n = m(m + 1) 4 + O(x 00 n )
where x 00 = 1 -4 m sin 2 π 2m+2 . In this section and the next two ones, we consider the case where n ≡ n m depends on m, and m → ∞. As in [START_REF] Berestycki | Limiting behavior of the distance of a random walk[END_REF], three main regimes appear: linear (n m = Θ(m))), cubic (n m = Θ(m 3 )) and intermediate (m ≪ n m ≪ m 3 ). This can be partly explained using the following simple bounds. In particular, if n ≡ n m and m → ∞,

I m,n = m(m + 1) 4 -Θ(m 2 x 00 n ).
Proof. These inequalities follow from the expression of I m,n given in Theorem 1. The upper bound is obtained by retaining only, in the sum over j and k, the term obtained for j = k = 0. The lower bound follows from |x jk | ≤ x 00 and [START_REF] Eriksen | Expected number of inversions after a sequence of random adjacent transpositions -an exact expression[END_REF].

Observe that x 00 = 1 -O(1/m 3 ). The lower bound on I m,n already shows that if n ≫ m 3 , then x 00 n = o(1) and I m,n ∼ m(m+1)

4

. Also, if n ∼ κm 3 for some κ > 0, then x 00 n ∼ α for some α ∈ (0, 1): then I m,n is still quadratic in m, but the upper bound shows that the ratio I m,n /m 2 will be less than 1/4. The other regimes correspond to cases where n ≪ m 3 . This section is devoted to the linear (and sub-linear) regime. We first state our results, and then comment on their meaning.

Proposition 11. Assume n ≡ n m = o(m). Then I m,n n = 1 + O(n/m). Assume n ≡ n m = Θ(m). That is, κ 1 m ≤ n m ≤ κ 2 m for two positive constants κ 1 and κ 2 . Then I m,n n = f (n/m) + O(1/m)
where

f (κ) = 1 2πκ ∞ 0 1 -exp(-8κt 2 /(1 + t 2 )) t 2 (1 + t 2 ) dt = j≥0 (-1) j (2j)! j!(j + 1)! 2 (2κ) j .
Comment. When n m is sub-linear, the inversion number equals its largest possible value, n, with high probability. For instance, I 1000,10 ≃ 9.9. When n m grows linearly with m, the expected inversion number is still linear in n, but with a ratio f (κ), where κ = n/m. This ratio decreases from f (0) = 1 to f (∞) = 0 as κ increases from 0 to ∞. The fact that f (0) = 1 is consistent with the sub-linear result.

Note that for a related continuous time chain, with inversion number D t m , it has been proved that when t = κm, the random variable D t m /t converges in probability to a function of t described in probabilistic terms [START_REF] Berestycki | Limiting behavior of the distance of a random walk[END_REF].

Proof. The starting point of both results is the following expression of I m,n , which corresponds to (5):

I m,n = 1 8(m + 1) 2 m k,j=0 c j + c k (1 -c j )(1 -c k ) 1 -1 - 4 m (1 -c j c k ) n . (20)
• Assume n m = o(m). Then 1 - 4 m (1 -c j c k ) n = 1 - 4n m (1 -c j c k ) + (1 -c j c k ) 2 O n 2 m 2 ,
uniformly in j and k. Thus

I m,n = n 2m(m + 1) 2 m k,j=0 (c j + c k )(1 -c j c k ) (1 -c j )(1 -c k ) + 1 4(m + 1) 2 m k,j=0 (1 -c j c k ) 2 (1 -c j )(1 -c k ) O n 2 m 2
(we have bounded |c j + c k | by 2) and the result follows using m k,j=0

(c j + c k )(1 -c j c k ) (1 -c j )(1 -c k ) = 2m(m + 1) 2 and m k,j=0 (1 -c j c k ) 2 (1 -c j )(1 -c k ) = (2m + 1)(m + 1) 2 . (21)
To prove these two identities, one may start from the following "basic" identity:

m j=0 1 1 -c j = (m + 1) 2 , (22)
which follows for instance from [START_REF] Eriksen | Estimating the expected reversal distance after a fixed number of reversals[END_REF].

• Assume now n m = Θ(m) and denote κ = n/m. Then

1 - 4 m (1 -c j c k ) n = exp (-4κ(1 -c j c k )) 1 + (1 -c j c k ) 2 O(n/m 2 ) .
Thus

I m,n = 1 8(m + 1) 2 m k,j=0 c j + c k (1 -c j )(1 -c k ) (1 -exp (-4κ(1 -c j c k ))) + 1 8(m + 1) 2 m k,j=0 (c j + c k )(1 -c j c k ) 2 (1 -c j )(1 -c k ) exp (-4κ(1 -c j c k )) O(n/m 2 ).
By [START_REF] Ph | The evolution of two stacks in bounded space and random walks in a triangle[END_REF], the absolute value of the second term above is bounded by

1 4(m + 1) 2 m k,j=0 (1 -c j c k ) 2 (1 -c j )(1 -c k ) O(n/m 2 ) = O(n/m) = O(1).
Recall that c j = cos (2j+1)π 2m+2 . Hence the first term in the expression of I m,n looks very much like a (double) Riemann sum, but one must be careful, as the integral

π 0 π 0 cos x + cos y (1 -cos x)(1 -cos y)
(1exp (-4κ(1cos x cos y))) dx dy diverges (the integrand behaves like 1/x around x = 0, and like 1/y around y = 0). Let us thus write

1 -exp (-4κ(1 -c j c k )) = [1 -exp (-4κ(1 -c j ))] + [1 -exp (-4κ(1 -c k ))] -[1 -exp (-4κ(1 -c j )) -exp (-4κ(1 -c k )) + exp (-4κ(1 -c j c k ))] ,
so that

I m,n = 2S 1 -S 2 + O(1) with S 1 = 1 8(m + 1) 2 m k,j=0 c j + c k (1 -c j )(1 -c k ) (1 -exp (-4κ(1 -c j )))
and

S 2 = 1 8(m + 1) 2 m k,j=0 c j + c k (1 -c j )(1 -c k ) × [1 -exp (-4κ(1 -c j )) -exp (-4κ(1 -c k )) + exp (-4κ(1 -c j c k ))] .
The first sum reads

S 1 = 1 8(m + 1) 2 m j=0 1 -exp (-4κ(1 -c j )) 1 -c j m k=0 c j + c k 1 -c k = 1 8 m j=0 1 -exp (-4κ(1 -c j )) 1 -c j (1 + c j ) - 1 8(m + 1) m j=0 1 -exp (-4κ(1 -c j )) 1 -c j , as m k=0 1 1 -c k = (m + 1) 2 and m k=0 c k 1 -c k = m(m + 1)
(see [START_REF] Flatto | Random shuffles and group representations[END_REF]). Both terms in S 1 are now bona fide Riemann sums. More precisely,

π m + 1 m j=0 1 -exp (-4κ(1 -c j )) 1 -c j (1 + c j ) = π 0 1 -exp(-4κ(1 -cos x)) 1 -cos x (1 + cos x)dx + O(1/m) uniformly in κ ∈ [κ 1 , κ 2 ]
, and similarly

1 m + 1 m j=0 1 -exp (-4κ(1 -c j )) 1 -c j = O(1).
Thus

S 1 = m 8π π 0 1 -exp(-4κ(1 -cos x)) 1 -cos x (1 + cos x)dx + O(1).
Similarly, S 2 is a (double) Riemann sum associated with a converging integral, and is O(1). Thus

I m,n = 2S 1 + O(1) = n 4πκ π 0 1 -exp(-4κ(1 -cos x)) 1 -cos x (1 + cos x) dx + O(1).
The integral can be rewritten by setting t = tan(x/2), and this gives the second result of the proposition, with the integral expression of f (κ). The expansion of f (κ) in κ is then routine. Assume n ≡ n m = Θ(m 3 ). That is, κ 1 m 3 ≤ n m ≤ κ 2 m 3 for two positive constants κ 1 and κ 2 . Then I m,n m 2 ∼ g(n/m 3 ) where

g(κ) = 1 4 - 16 π 4   j≥0 e -κπ 2 (2j+1) 2 /2 (2j + 1) 2   2 .
Comment. When n ≫ m 3 , the inversion number equals, at least at first order, its "limit" value m(m+1) 4

, which is the average number of inversions in a permutation of S m+1 taken uniformly at random. When n = Θ(m 3 ), the inversion number is still quadratic in m, but with a ratio g(κ), where κ = n/m 3 . This ratio increases from 0 to 1/4 as κ goes from 0 to ∞, which makes this result consistent with the super-cubic regime. Note that for a related continuous time chain, with inversion number D t m , it has been proved that when t = κm 3 , the random variable D t m /m 2 converges in probability to a function of t described in probabilistic terms [START_REF] Berestycki | Limiting behavior of the distance of a random walk[END_REF].

As said above, I m,n ∼ m 2 /4 as soon as n ≫ m 3 . Whether the next term in the expansion of I m,n is exact, that is, equal to m/4, depends on how n compares to m 3 log m, as shown by the following proposition. In particular, it clarifies when I m,n = m(m+1) Thus if c < 1/π 2 , there exists γ > 0 such that

I m,n ≤ m(m + 1) 4 -Θ(m 1+γ ),
while if c > 1/π 2 , there exists γ > 0 such that

I m,n = m(m + 1) 4 -O(m 1-γ ).
For the critical value c = 1/π 2 , the following refined estimate holds: if n ≡ n m ∼ 1/π 2 m 3 log m+ αm 3 + o(m 3 ), then

I m,n = m(m + 1) 4 - 16m π 4 e -απ 2 (1 + o(1)).
Proof of Proposition 12. The first result is a direct consequence of Lemma 10, given that x 00 = 1 -O(1/m 3 ). Assume now that n m = Θ(m 3 ). We start from the following expression of I m,n , which corresponds to (4):

I m,n = m(m + 1) 4 - 1 8(m + 1) 2 m k,j=0 c j + c k (1 -c j )(1 -c k ) 1 - 4 m (1 -c j c k ) n . (23)
Let M ≡ M m be an integer sequence that tends to infinity in such a way M m = o( √ m). We split the sum over j and k into two parts: j ≤ M , k ≤ M in the first part, j > M or k > M in the second part. Let us prove that the second part can be neglected. We have:

|SUM 2 | := j or k>M c j + c k (1 -c j )(1 -c k ) x jk n ≤ j or k>M 2 (1 -c j )(1 -c k ) ≤ 4 m k=0 1 1 -c k m j=M+1 1 1 -c j .
By [START_REF] Flatto | Random shuffles and group representations[END_REF], the sum over k is O(m 2 ). Moreover, the sum over j is a Riemann sum, and the function x → 1/(1cos x) is decreasing between 0 and π, so that

m j=M+1 1 1 -c j ≤ m + 1 π π (2M +1)π 2m+2 dx 1 -cos x = m + 1 π tan (2M+1)π 4m+4 = O m 2 M .
This gives

SUM 2 = O m 4 M = o (m 4 ).
Let us now focus on small values of j and k. The following estimates hold uniformly in j and k, when 0 ≤ j, k ≤ M :

c j + c k (1 -c j )(1 -c k ) = 128(m + 1) 4 (2j + 1) 2 (2k + 1) 2 π 4 1 + O M 2 m 2 , (24)
= 128m 4 (2j + 1) 2 (2k + 1) 2 π 4 (1 + o(1)) , 1 - 4 m (1 -c j c k ) n = exp - nπ 2 ((2j + 1) 2 + (2k + 1) 2 ) 2m(m + 1) 2 1 + O n M 4 m 5 (25)
= exp -nπ 2 ((2j + 1) 2 + (2k + 1) 2 ) 2m 3

(1 + o(1)) .

We have used the fact that n = O(m 3 ) and M = o( √ m). Hence,

SUM 1 := M k,j=0 c j + c k (1 -c j )(1 -c k ) x jk n = 128m 4 π 4   M j=0 e -nπ 2 (2j+1) 2 2m 3 (2j + 1) 2   2 (1 + o(1)) .
Finally, observe that

M j=0 e -nπ 2 (2j+1) 2 2m 3 (2j + 1) 2 = h n m 3 + O(e -n M 2 m 3 ) = h n m 3 (1 + o(1))
where

h(κ) = M j=0 e -κπ 2 (2j+1) 2 /2 (2j + 1) 2 .
Writing

I m,n = m(m + 1) 4 - 1 8(m + 1) 2 (SUM 1 + SUM 2 )
then yields the result.

Proof of Proposition 13. The first result follows from the second identity of Lemma 10, using

n = cm 3 log m(1 + o(1)) and log x 00 = - π 2 m 3 (1 + O(1/m)).
If c < 1/π 2 , we choose ǫ such that γ := 1-cπ 2 -ǫ > 0, and the second result follows. If c > 1/π 2 , we choose ǫ such that γ := cπ 2 -1 -ǫ > 0, and the third result follows.

Assume now n = 1/π 2 m 3 log m + αm 3 + o(m 3 ). We will show that the sum over j and k occurring in [START_REF] Janse Van Rensburg | Partially directed paths in a wedge[END_REF] is dominated by the contribution obtained for j = k = 0. This contribution reads 2c 0

(1 -c 0 ) 2 x 00 n = 128m 3 π 4 e -απ 2 (1 + o(1)).
The contribution obtained for j = k = m is

- 2c 0 (1 + c 0 ) 2 x 00 n = O(1/m).
Finally, for all other values of j and k, there holds 0 < x jk ≤ x 01 . Hence 0≤k,j≤m (k,j) =(0,0),(m,m)

c j + c k (1 -c j )(1 -c k ) x jk n ≤ 2 x 01 n 0≤k,j≤m 1 (1 -c j )(1 -c k ) = 2(m+1) 4 x 01 n = O(1/m).
We have first used the identity ( 22 Comment. By Propositions 11 and 12, one has

I m,n √ mn → √ κf (κ) if n m ∼ κm, while I m,n √ mn → g(κ) √ κ if n m ∼ κm 3 .
It can be proved that

lim κ→∞ √ κf (κ) = 2 π = lim κ→0 g(κ) √ κ ,
which makes all three regimes consistent. Note that for a related continuous time chain, with inversion number D t m , it has been proved that when m ≪ t ≪ m 3 , the random variable D t m / √ mt converges in probability to 2/π [START_REF] Berestycki | Limiting behavior of the distance of a random walk[END_REF].

Proof. The proof mixes arguments that we have already used for small and large times. As in the small time case, we start from [START_REF] Eriksson | Expected number of inversions after a sequence of random adjacent transpositions[END_REF]. As in the large time case, we split the sum over j and k into two parts: j ≤ M , k ≤ M in the first part, j > M or k > M in the second part. Here, M = M m is an integer sequence that satisfies

m 3 n ≪ M ≪ min m 3 n , m m n 1/4
. Such a sequence exists under the assumptions we have made on n.

We begin by bounding the second part. As in the large time case, one easily obtains

|SUM 2 | := j or k>M c j + c k (1 -c j )(1 -c k ) (1 -x jk n ) ≤ j or k>M 2 (1 -c j )(1 -c k ) = O m 4 M . (26)
We now focus on the case j, k ≤ M . The estimates ( 24) and ( 25) still hold uniformly in j and k. This gives this may be rewritten as

SUM 1 := M k,j=0 c j + c k (1 -c j )(1 -c k ) (1 -x jk n ) = 128m 4 π 4 M j,k=0 1 -exp -nπ 2 ((2j+1) 2 +(2k+1) 2 ) 2m 3 (2j + 1) 2 (2k + 1) 2 (1 + o(1)) .
SUM 1 = 128m 4 π 4 2S M M k=0 1 (2k + 1) 2 -S M 2 (1 + o(1)) (27)
with

S M = M j=0 1 -exp -nπ 2 (2j+1

Now recall that

I m,n = 1 8(m + 1) 2 (SUM 1 + SUM 2 ).

Combining the above estimate of SUM 1 and the bound (26) on SUM 2 gives the announced estimate of I m,n .

Perspectives

Many interesting Markov chains on groups have been studied, and it is natural to ask to which similar problems the approach used in this note could be adapted. To make this question more precise, let us underline that such problems may involve changing the dynamics of the chain, changing the statistics under consideration, or changing the underlying group.

Changing the chain. Already in the symmetric group S m+1 , a number of different dynamics have been considered. One can for instance multiply at each step by any transposition [START_REF] Berestycki | A phase transition in the random transposition random walk[END_REF][START_REF] Diaconis | Generating a random permutation with random transpositions[END_REF][START_REF] Schramm | Compositions of random transpositions[END_REF], or consider cyclic permutations by allowing multiplications by the transpositions (0, 1), (1, 2), . . . , (m -1, m) and (m, 0), or only allow transpositions with the largest element [START_REF] Flatto | Random shuffles and group representations[END_REF][START_REF] Diaconis | Comparison theorems for reversible Markov chains[END_REF], or perform block transpositions [START_REF] Bóna | The average number of block interchanges needed to sort a permutation and a recent result of Stanley[END_REF], or shuffle the elements in various ways [START_REF] Bayer | Trailing the dovetail shuffle to its lair[END_REF], etc. See for instance [START_REF] Diaconis | The cutoff phenomenon in finite Markov chains[END_REF].

Changing the statistics. We have focused in this paper on the inversion number. It is an interesting parameter from a probabilistic point of view, as it gives an indication of whether the chain may be mixed at time n or not. However, in certain biological contexts, it may be more sensible to estimate instead the natural "distance" between π (0) and π (n) , defined as the minimum number of chain steps that lead from one to the other [START_REF] Berestycki | A phase transition in the random transposition random walk[END_REF][START_REF] Berestycki | Limiting behavior of the distance of a random walk[END_REF]. For the chain studied here, this coincides with the inversion number, but for other dynamics it will be a different parameter. For instance, if we multiply by any transposition, one step suffices to go from 012 to 210, whereas the inversion number of 210 is 3. Clearly, the approach we have used here relies heavily on the possibility of describing in simple terms the evolution of the parameter under consideration (as we did by the combination of ( 6) and Lemma 3). Note that an explicit expression for the expected distance after n (non-necessarily adjacent) transpositions has been given in [START_REF] Eriksen | Estimating the expected reversal distance after a fixed number of reversals[END_REF].

Changing the group. Although many different chains in many different groups have been considered, we are primarily thinking of classical families of (finite or affine) Coxeter groups, because the inversion number admits in these groups a natural generalization (the length) which has usually a simple description [START_REF] Björner | Combinatorics of Coxeter groups[END_REF]Chap. 8]. According to [START_REF] Eriksen | Expected number of inversions after a sequence of random adjacent transpositions -an exact expression[END_REF], an explicit formula, due to Troili, is already known for the average length of the product of n generators in the affine group Ãm . Another solved case is the hypercube Z m 2 [START_REF] Diaconis | Asymptotic analysis of a random walk on a hypercube with many dimensions[END_REF].

(- 1 )

 1 k (m -2ℓ) 2⌈s/2⌉ -1 ⌈s/2⌉ + ℓ + k(m + 1)andh s,m = j∈Z (-1) j 2⌊s/2⌋ ⌊s/2⌋ + j(m + 1)

Lemma 10 .

 10 For m ≥ 3 and n ≥ 0, there holds m(m + 1) 4 (1 -x 00 n ) ≤ I m,n ≤ m

4 .

 4 Large times: cubic and beyond Proposition 12. If n ≡ n m ≫ m 3 , then I

4 +

 4 o(m): this happens beyond n ∼ 1/π 2 m 3 log m. Proposition 13. Let c > 0, and n ≡ n m ∼ cm 3 log m. For all ǫ > 0, m(m + 1) 4 -O(m 2+ǫ-cπ 2 ) = I m,n ≤ m(m + 1) 4 -Θ(m 2-ǫ-cπ 2 ).

5 .

 5 ), then n = 1/π 2 m 3 log m + αm 3 + o(m 3 ), and log x 01 = -5π 2 /m 3 (1 + o(1)). The estimate of I m,n follows. The intermediate regime Proposition 14. Assume m ≪ n m ≪ m 3 . Then, denoting n ≡ n m ,

Using 1 -

 1 e a+b = (1 -e a ) + (1 -e b ) -(1 -e a )(1 -e b )

  )

2

 2

2m 3 (2j + 1) 2 . 1 -

 321 The sum S M is close to a Riemann sum. More precisely, e -x 2x 2 dx (1 + o(1)) = √ π + o(1).In particular, S M = o(1). Returning to[START_REF] Saloff-Coste | Random walks on finite groups[END_REF], let us finally note that

Acknowledgements. The author thanks an anonymous referee, Niklas Eriksen and David Wilson for their suggestions, comments, and help with the bibliography.


Muriel Mescam

Marek Kretowski

Johanne Bezy-Wendling
email: <johanne.bezy-wendling@univ-rennes1.fr

Ltsi

Multiscale model of liver DCE-MRI towards a better understanding of tumor complexity

Keywords: Computational modeling, Liver, Blood, MRI simulation, Magnetic resonance imaging, Mathematical model, Portals, Tumors, complex system, image analysis, liver tumors

Multiscale model of liver DCE-MRI towards a better

understanding of tumor complexity.

Introduction

Combining image analysis and computational modeling is one of the possible means to understand complex living systems . The [1 ] liver, and its pathological modifications, such as tumors, is a system with very complex dynamics, especially related to its original double vascular supply. Understanding the anatomical (structural) and functional processes underlying the observations, whatever the modality, is a necessary step to a better characterization of lesions and consequently to propose the best therapy to the patient. The definition of a model is always done taking into account the way it will be used. In this study, we propose a model-based image analysis approach whose objective is to give a physiological sense to parameters extracted from Dynamic Contrast Enhanced Magnetic Resonance Imaging (DCE-MRI) of the liver. These features result from many interacting components associated to organ characteristics (tissue, vascularization) and to image formation parameters (contrast agent, sequence). We focused on the liver and its most widespread malignant tumor, which is the HepatoCellular Carcinoma (HCC). Worldwilde, it is one of the most common malignant primitive tumors; its incidence has been regularly increasing over the past years without slowing down , and the prognosis is generally bad. Several kinds

[2 ] [3 ]
of therapy, including resection, radiofrequency ablation, chemotherapy, and antiangiogenic molecules, can be proposed to patients.

Regardless of the therapeutic solution, the treatment will be better if the tumor detection and characterization are realized during the first stages of its development, when the tumor size is less than 2 or 3 cm, and when the hepatic function is still correct. MRI is a non invasive and very safe observation method, particularly well adapted to evaluate hepatic lesions because it provides a good lesion to liver contrast [START_REF] Sahani | Imaging the liver[END_REF] . The use of contrast agents, and especially those based on gadolinium, provides an estimation of the tumor hypervascularization, and ] perfusion modifications, which are the main markers of the tumoral growth and malignancy. Rapid acquisition MRI sequences allow us to follow the enhancement just after injection and then to show the differences of biodistribution of the contrast medium between normal and tumoral tissues. However, the correct interpretation of MR images of HCC, and consequently the resulting diagnosis, requires a good comprehension of the physiological phenomena involved in the contrast agent distribution, and of the physical processes related to the image formation. DCE-MRI, based on sequences of MR images acquired after contrast agent injection, is used in clinical practice to assess liver tumor malignancy. Classical pharmacokinetic models are usually fitted to perfusion curves extracted from DCE-MRI, in order to determine physiological markers such as and , two constants that are used to classify tumors. The contribution of such analysis to K ep K trans the diagnosis and classification of liver lesions has been proven and is not discussed here. However, as far as the comprehension of [START_REF] Caldeira | Liver tumor assessment with DCE-MRI[END_REF] physiological phenomena is concerned, DCE-MRI is somewhat limited regarding the amount and variety of the available data. Since the technique is based on real data, it is impossible to control every step that leads to the observation (image contrast) from the modifications occuring in the organ, to the image formation. Physiological modeling represents a non invasive way to control parameters that would be difficult or even impossible to access with real experiments. In this aim, we propose the framework presented in , where the Fig. 1 coupling between a model of the tumoral liver and a model of dynamic MRI is presented. This multiscale model of liver DCE-MRI provides data we already have a precise knowledge on, in order to relate physiological phenomena to simulated observations. The general aim here is not to replace DCE-MRI but to provide understanding of the mechanisms of production of the DCE-MRI signal, based on a detailed model of contrast agent transport in the liver. Two main steps can be distinguished in this model-based approach: the i) macroscopic model of the liver, including tissue and vessels, where the structure and the geometry of the vascular network are simulated in normal and tumoral cases (see section A), and the dynamic image formation model including transport of the contrast agent and ii) microscopic characteristics related to this transport at the capillary level (see section B). Only a brief presentation of the macroscopic model is given in the next section, compared to the more detailed description of the dynamic image formation model. Indeed, the most important and original part of the work presented here concerns the model of parenchymal distribution of MRI contrast agents. Detailed explanations on the macroscopic model and some of its applications can be found in our previous studies (

). Results are presented [START_REF] B Zy-Wendling É | A 3-D dynamic model of vascular trees[END_REF] - [START_REF] Kretowski | Simulation of biphasic CT findings in hepatic cellular carcinoma by a two-level physiological model[END_REF] in the last section: a sensitivity analysis of the model for two main physiological parameters (blood flow and permeability), and an MRI simulation of HCC in different pathophysiological conditions describing the main steps of carcinogenesis from normal liver to poorly differenciated HCC (HCCp) ( ). Fig. 2 Model based approach

Model of the hepatic tumor growth

In their very recent paper , Lloyd present a generic framework for modeling tumor growth and emphasize that two main [START_REF] Lloyd | A computational framework for modelling solid tumor growth[END_REF] et al.

components, which are tissue and vascularization, have to be integrated in the model. We already proposed a very comparable framework in 1999 . However, when considering an organ such as the liver, this kind of generic approach must be specified because of the very [START_REF] B Zy-Wendling É | A 3-D dynamic model of vascular trees[END_REF] particular irrigation of this organ and its tumors (double blood supply by the hepatic artery and the portal vein) and especially if this model is used to simulate dynamic images. This section deals with the macroscopic part of the hepatic computational model. The modeled liver consists of two main components: the tissue (which corresponds to the parenchyma) and the vascular network. The tissue modeling is first explained, taking into account the liver growth, and the possible differences in tissular characteristics. Then, in the vascular model, the structure and geometry of the three hepatic vascular trees (hepatic artery, portal vein and hepatic vein) are precisely considered, as they are essential elements that appear clearly in images and also constitute the pipes that transport blood and contrast agent molecules toward capillaries (place of exchanges between blood and tissue).

Tissue modeling

The tissue is represented by a set of Macroscopic Functional Units (MFU) evenly (but randomly) distributed inside the tridimensional organ shape. Each MFU is a small, fixed size part of the tissue. It is characterized by its spatial, relative position inside the shape, and by its class, which determines most of functional/structural properties (rhythm of mitosis/necrosis), and physiological features (e.g. blood flow rate). Several classes of MFUs can be defined to differentiate functional (or pathological) regions of tissue (normal/tumoral).

Furthermore, the MFU class can be modified over time, simulating the possible evolution of a disease (e.g. tumor). Certain parameters (such as blood flow rate) of a class are described by their distribution, and are randomly chosen for each new MFU; (e.g. This is the case for the terminal blood flow whose natural variability is modeled in such a way).

Vascular network modeling

A literature review on vascular models is given in and in . Some of them are generic, and others have been specifically

[6 ] [START_REF] Kretowski | Leondes C . Modeling for medical image analysis: framework and applications[END_REF] developed for a particular organ (such as the heart in , or the kidney in ). The tissue is supplied by a vascular network,

[12 ] [13 ] [7 ]
composed of three trees, with blood going from the arterial and portal trees to the venous one, through the MFUs. A vessel segment (part of vessel between two consecutive bifurcations) is represented by an ideal, rigid tube with fixed radius, wall thickness and position. In the model, all vessels except the capillaries are distinguishable. On the contrary, the geometry of capillaries is not considered, as these smallest vessels are hidden in the MFU (see .B.1.). Based on morphometrical investigation dealing with bigger vessels (e.g. conducted by Zamir " " §

and Chee ), it is assumed that a single vascular structure has a form of a binary tree. Blood is considered as a Newtonian fluid, with [START_REF] Zamir | Branching characteristics of human coronary arteries[END_REF] constant viscosity, and its flow is governed by Poiseuille s law (which relates the pressure difference between two extremities of a vessel

'
with the blood flow rate, length and radius of that vessel). At each bifurcation the law of matter preservation is in effect (the quantity of blood, which enters a bifurcation by a vessel, leaves it through its two descendant vessels). Another constraint deals with the decreasing vessel radii in the vascular tree, creating a relationship between the radius of a vessel and radii of its two descendants. Assuming that positions of all vessels are fixed, the vascular trees consistency is assured. A vascular tree is consistent if: it has the same blood pressure

' i)
and fixed blood flow in all terminal vessels attached to MFUs, the Poiseuille s law in each vessel, matter preservation and bifurcation ii) ' laws in each bifurcation are fulfilled. A computationally effective method for the consistency assurance in the vascular tree is described in . Newly appeared MFUs are initially not perfused by the existing vascular network. The closest vessels (in each tree) sprout towards it. [START_REF] Kretowski | Fast 3D modeling of vascular trees[END_REF] A new bifurcation is created, and its geometry is controlled by local minimization of additional blood volume needed for the MFU perfusion. In order to find out the optimal configuration, each candidate vessel temporarily creates a bifurcation perfusing the MFU and the configuration volume is calculated. Additionally, the problem of avoiding possible collisions between perfusing vessels is taken into account. It concerns both intersections of vessels coming from the same tree or from two trees. Finally, the configuration of new bifurcations from all the trees with the lowest sum of volumes is chosen to permanently perfuse the MFU. The next step is devoted to a recalculation of vessels characteristics (i.e. blood flow, pressures and radii) in the whole tree to fulfill all constraints (i.e. physical and ' physiological laws). The fast-updating method used to accelerate this time consuming operation is presented in . The development of a [START_REF] Kretowski | Fast 3D modeling of vascular trees[END_REF] pathological process (e.g. tumor) is obtained by means of the conversion mechanism: a conversion represents the period when the current class of MFU inside a chosen region (e.g. sphere) can be changed (with the given probability) to another class. The sequence enables to model various stages of pathology development, when macro-cells characteristics evolve gradually.

Some illustration of simulated vascular models can be seen on , where the vascular network of a normal liver is represented.

Fig. 3

Model of dynamic MRI

Parenchymal biodistribution of MRI contrast agents

The simulation of the propagation of an MRI contrast agent in the liver parenchyma involves a precise mathematical description of its distribution in every sub-unit of the tissue. As far as gadolinium-based extracellular contrast agents are concerned, the main components to deal with are the plasma region and its surrounding interstitial space between capillaries and hepatic cells , called the space of Disse. In

- -
this aim, a first compartmental Physiologically Based PharmacoKinetic (PBPK) model has been developed and is described in . This [START_REF] Mescam | A physiologically-based pharmacokinetic model of vascular-extravascular exchanges during liver carcinogenesis: application to MRI contrast agents[END_REF] five-compartments model considers the terminal branches of the vascular trees (hepatic arteriole, portal venule and hepatic venule) as full compartments, which communicate with the space of Disse. However, since changes in the permeability of these normally

" "
non-exchanging micro-vessels, are not well defined today, especially in the context of hepatocarcinogenesis, these micro-vessels are assumed impermeable in the present study. As a matter of fact, they only constitute the double arterio-portal input and the hepatic output of the microscopic model. The capillary network (sinusoids) is consequently the only media where exchanges occur, through its fenestrated endothelium ( ). Fig. 4 These transcapillary exchanges, occuring through the pores and fenestrations of the endothelial cells, are described by a pore model. modified by Bassingthwaighte . These functions depend on the molecule/pore radii ratio ( / ). The solute partition coefficient [START_REF] Bassingthwaighte | A practical extension of hydrodynamic theory of porous transport for hydrophilic solutes[END_REF] α α k = r r k is defined by (1 ) . Where the P clet number, , is a dimensionless number measuring the ratio of transport by convection to transport by diffusion. The é Pe total fluid and solute fluxes, respectively and through the sinusoidal membrane, are thus given by the sum of fluxes through each kind J v J s of pathway : k

One important improvement of the present model, compared to the compartmental model presented in , relies on the axial [START_REF] Mescam | A physiologically-based pharmacokinetic model of vascular-extravascular exchanges during liver carcinogenesis: application to MRI contrast agents[END_REF] distribution of its plasma and interstitial (space of Disse) regions. Such additional feature allows for concentration gradients along the capillary length and, in this sense, gives a better representation of the real transport processes, than compartmental models where well-stirred regions are filled instantaneously. The new model, presented in , is then comparable to BTEX units (Blood Tissue Fig. 4 EXchange units), defined by Bassingthwaighte et al in , where the convection-permeation-dispersion mechanisms are described in the [START_REF] Bassingthwaighte | Blood-tissue exchange via transport and transformation by capillary endothelial cells[END_REF] following system of partial differential equations (PDEs):

Where is the capillary length and is the concentration at the input of the plasma compartment, resulting from the dual L C in arterio-portal supply ( ). 9

In addition, the model allows for variations in regions volume, reflected in pressure variations (( ) and ( )). Such phenomena has Where is the lymphatic conductance. K L Mass conservation, in this system, is consequently ensured by the following system of ordinary differential equations:

Coupling the two levels, in order to obtain a multiscale model of contrast-agent propagation, necessitates the calculation of the global concentration in each functional unit of the macroscopic model (MFU). This is done by ( ), from which an averaged concentration is 15 computed.

MRI acquisition simulation

A 3D MRI simulator, named SIMRI and recently developed by Benoit-Cattin , was used to generate simulated images of the et al. [START_REF] Benoit-Cattin | The SIMRI project: a versatile and interactive MRI simulator[END_REF] liver. The simulation system is based on the Bloch equation resolution. A 3D virtual object constitutes the input of the simulator. Each Where relaxivities, , are constants that are contrast agent dependent. r 1,2

Simulation results and discussion

Simulation procedure

All the results presented in this paper were obtained following the same simulation procedure. After coupling both macro-and microscopic models, a virtual tumor is generated, based on a simulated adult liver. Propagation of a Gd-DOTA contrast agent (Dotarem , ® Guerbet, France), currently used in clinical routine, is simulated, using, as entry functions, real concentration-vs-time curves obtained from a DCE-MRI study. The arterial input function (AIF) was measured from a ROI selected in the aorta (arteries are thin and hardly visible on images), on a sequence of MR images acquired after an intravenous bolus injection of Gd-DOTA, with a time step of approximately 10 seconds. The time needed by the mixing blood-contrast agent to reach the main branch of the hepatic artery is neglected since the aorta-hepatic artery pathway is straightforward compared to the aorta-portal vein pathway. In order to account for the portal delay as accurately as possible, a second ROI was selected in a portal vein (larger than the hepatic artery). Once the agent concentration is computed in every voxel as a function of time, post-injection T , T and PD maps are created at typical arterial and portal phases (i.e. brightening tissues with lower T , thus explaining the observed contrast ( ). Later, at portal phase, Gd-DOTA concentration has 1 Fig. 5 decreased and the contrast appears hypointense compared with the normal surrounding tissue ( bottom). Fig. 5 Sensitivity analysis A large number of parameters controls each level of the multiscale model, thus enhancing the difficulty of its identification. However, it is possible and relevant to evaluate its sensitivity to some particular parameters, such as those related to the development of the HCC.

For instance, inflows and capillary permeability play an important role in the evolution of the pathology from a benign regenerative nodule to a poorly-differentiated HCC (

). First, the influence of each parameter is separately evaluated through a sensitivity analysis. The Fig. 2 output of the model, i.e. simulated MR images of the liver, is obtained for increasing values of arterial inflow, , and then for increasing F ha values of the permeability-surface area product, which is referred to by the notation in the following study. In each case, enhancement PS concentration curves are simulated at the microscopic scale and corresponding simulated bi-phasic (arterial and portal) images are depicted. A tumoral region is considered, and the ratio of the mean signal in the tumor to the mean signal in the surrounding normal

" " R S " "
tissue is compared to the ratio of Gd-DOTA concentration in both tissues. R C Concentration-vs-time curves obtained for varying inflows , are presented on the graph of . These concentrations are F ha Fig. 6 simulated from the distributed model of contrast agent propagation, separately uncoupled from the vascular model, for the classical gadolinium chelate (Gd-DOTA) ( ). In order to observe the influence of the arterial inflow only, the portal supply is inhibited, and the 15 parameters controlling the permeability are chosen so that the latter is high and the model is therefore not diffusion-limited. Concentration profiles obtained for increasing arterial flow are well aligned with the established arterialization phenomena that occurs in liver tumors, showing an enhancement peak at the so-called arterial phase. In addition, a non-linear behavior of the concentration, as a function of , F ha can be observed at the portal phase, where the concentration in the liver grows as the inflow grows, to reach a maximum value above which it decreases again. Values of the concentration in the liver, at two particular times after injection, corresponding to arterial (24 seconds) and portal (50 seconds) phases, are extracted from these curves, and the ratio / is depicted as a function of , on C tumoral C normal F ha the right graphs of , respectively. Such representation not only confirms the previous observations regarding the linearity of Fig. 7.a and b ( ) at the arterial phase, and its non-linearity during the portal phase, it also allows expectation of some particular contrast dynamics C F ha on simulated MR images. Indeed, a linear relationship exists between the contrast agent concentration and relaxation parameters ( ). 16

However, an exact quantification of the contrast enhancement is not straightforward, since the relationship between the tissue relaxations and the effective NMR signal intensity also depends on the applied sequence and on the acquisition parameters. Nevertheless, constrast dynamics can at least be explained by the concentration dynamics, as shown on , where the behavior of with matches Fig. 7. a and b R S F ha that of . These observations are illustrated by the corresponding set of simulated MR images, on which the contrast tumoral/normal R C tissue confirms the quantitative measurements. Furthermore, the nonlinearity of the contrast evolution with the inflow at the portal phase, whereas it is linear at the arterial phase, suggests an easier interpretation of real images at the arterial phase, concerning the relationship between the observation on the image and the inflow it reflects. In addition, images were simulated for three different slice thicknesses (1 mm, 4 mm and 8 mm). From the measurements presented on , it seems evident that the larger the slice, the better the contrast. For Fig. 7 very thin slices (1 mm), it can even become relatively difficult to give a correct interpretation of the relation between the contrast and the inflow. Since the intensity of the signal in the tumor region is higher than in the normal tissue, the SNR increases faster in this region, with the slice thickness, which may explain the higher contrast in this cases. In addition, in thicker slices, tumor regions are more vascularized whether normal tissue is more homogeneous whatever the slice thickness. However, the important issue here is not slice thickness per se but SNR, which can be improved by a number of means other than changing slice thickness, including changing the in-plane resolution, changing the receive coil, changing field strength, changing the number of averages.

In the following simulations, the relation between the contrast enhancement and the capillary permeability is studied. This time, results are presented for only one slice thickness. Slices of 4 mm are a good compromise between the quality of the contrast and the computation time. In addition, this is commonly used in clinical routine, as shown in qualitative observations presented in III.B. Among the numerous § physiological parameters that control the permeability, variations in the quantity of pores and fenestrations on the membrane of sinusoids are directly related to the capillarization phenomena (progressive disappearance of fenestrations), which is one of the main mechanisms involved in hepatocarcinogenesis. However, it is more common among the modeling community to deal with the permeability-surface area product, , an upper scale parameter that generally describes solute transport through a membrane barrier. The pore model described in PS the first section does not allow a direct calculation of from the pores density. As a matter of fact, an estimated value of is extracted PS ' PS after fitting concentration curves, simulated for varying values of the densities of small and large pores and fenestrations, to a simple BTEX model in which the diffusion term is replaced by ( ) in and . Simulation results showed a different S J

p s PS C p -C d equations (7) (8)
concentration-time course when setting a low input arterial flow, compared to a high one ( ). Effectively, for a low inflow, Fig. 8. a and b the concentration in the liver parenchyma is clearly not dependent on the capillary permeability at both arterial and portal phases and is thus flow-limited. Consequently, almost no variation can be noticed in the signal intensity with this parameter (

). In addition, the Fig. 9 contrast agent concentration is very close in both classes of tissue, which explains the poor contrast on corresponding simulated images at the two phases. On the contrary, at high inflow, a linear relationship exists between the concentration and the permeability, and then indirectly between the signal intensity and the permeability. This relationship is no longer linear at portal phase. Indeed, when capillaries are weakly permeable, the contrast agent is mainly confined to the vascular space. Consequently, the solute is primarily driven from the arterio-portal entrance to the hepatic venous extremity of the network, by convection and dispersion. The washout is therefore completed much faster. When the permeability is no longer negligible, the diffusion process tends to hinder this washout. However, when it becomes significant, i.e. above some threshold, the steady-state is rapidly reached and the model is no longer diffusion-limited.

Although the previous results show the evidence of the existence of relationships between phenomena observed on images and some physiological parameters that are markers of the tumoral development, the kinetics of contrast enhancement remain difficult to explain. Furthermore, in this sensitivity analysis, only one parameter was varied at a time. However, the real processes involved in the development of a pathology and its evolution are obviously much more complex and imply the variation of a set of parameters instead of a single one.

Application to the hepatocarcinogenesis

The multi-step hepatocarcinogenesis describes the different stages of the development of a liver hypervascularized tumor, from benign regenerative nodules (RN) to poorly-differentiated hepatocellular carcinoma (HCCp), through respectively low-and high-grade dysplastic nodule (DNl,h), and medium-differentiated HCC (HCCm) (

). Parallel modifications of several physiological characteristics affecting Fig. 2 both the tissue and the vascularization are involved in the evolution of the lesion from one step to the other. For instance, we have seen previously that arterial/portal inflows and capillary permeability may affect the observed contrast differently. Gadolinium-enhanced T 1 -weighted GRE MRI is specifically used to detect HCC because of the highly arterialized neovascularization of this lesion, compared to DNs and RNs. However, the same sequence (same acquisition parameters as used in previous simulations) has been used here to simulate MR images of all different grades, so that only physiological parameters are responsible for differences observed between the simulated images.

In this section, simulation results are still given at two scales: concentration-time curves at the microscopic scale ( ), and Fig. 10 the study they are extracted from. However, a general behavior of portal and arterial supplies during hepatocarcinogenesis has been established and enables us to define reasonable values . In the present study, arterial and portal flows were thus chosen taking into [START_REF] Matsui | Detection and characterization of small hepatocellular carcinoma[END_REF] account normal inflows and relative variations described on the diagram of Figure 8 in . Similarly, the number of fenestrations is [START_REF] Matsui | Detection and characterization of small hepatocellular carcinoma[END_REF] known to decrease progressively in capillary wall so values have been chosen relatively to known normal values .

[26 ]

In term of contrast tumor/normal tissue, simulated images are coherent with common observations brought on typical post-contrast T 1 -weighted (T1W) imaging. For benign lesions such as regenerative nodules (NR) and eventually low grade dysplastic nodules (DNl), arterial and portal inflows are close to normal ones, and sinusoids are still fully permeable, which explains the difficulty, not only in their characterization but also in their detection since there is no arterial enhancement in Gd-DOTA-enhanced T1W images. At early stages of the malignancy, the arterial inflow slightly diminishes, resulting in a hypointense lesion on T1W GRE MR images at arterial phase (Fig.

.DNh). On the contrary, when the nodule evolves to a carcinoma, the decrease in the portal supply is progressively compensated by an 10.b increase in the arterial one, in a so-called neoangiogenesis process. Consequently, the tumor appears hyperintense at the arterial phase and hypointense at the portal one ( ), when the concentration is higher in the normal tissue ( ). This enhancement is Fig. 10.b Fig. 10.b particularly strong for highly graded carcinoma (HCCp), for which the portal supply is completely absent in the tissue.

Although the capillarization phenomena tends to reduce the contrast enhancement at the arterial phase, as shown in the last paragraph, this enhancement is higher for the HCCp, for which, capillaries are almost impermeable. Indeed, on the graph of , we have seen Fig. 8.b that for a high arterial inflow, the concentration is high, whatever the permeability. This example illustrates rather well the difficulty of an interpretation from the single images, and the utility of a multiscale analysis.

Conclusion

The multiscale model of liver DCE-MRI presented in this paper has already shown good qualitative results regarding the coherence of contrast enhancements in pathological situations, especially in the case of a HCC. Furthermore, its ability to provide a multiscale analysis by coupling dynamic concentration curves to simulated MR images contributes to a better understanding of the physiological phenomena involved in the hepatocarcinogenesis. In the attempt to obtain more information on tumors, it is crucial for the model to describe the physiology and the dynamics of tissues that affect the contrast, as precisely and as realistically as possible. Here, the axially-distributed model generates more realistic time distribution of the contrast agent in the liver parenchyma, than commonly used compartmental models, in which each compartment is a well-stirred tank. In addition, parameters used in the current model have a direct physiological, physical or geometrical meaning. However, an essential identification step is still to be achieved, before being able to provide a useful tool for the characterization and, eventually, for the diagnosis of the tumor. This problem might be partially overcome by experimental studies on animal models of HCC. In addition, a precise study of the evolution of other parameters related to the images, such as texture parameters, which has already shown an interest in tissue characterization, is foreseen. The determination of relationships between these features and the physiological parameters related to the model would thus help in solving the inverse problem, by representing bridges between observations on real images and the underlying physiological phenomena they reflect. In this aim, an initialization of the vascular model by vessel branches segmented from acquired images would allow the creation of a patient-specific model, more efficient in the characterization of the tumor.

  From physico-chemical parameters related to the targeted molecule (radius and molecular diffusion ), and physical parameters related r D to the capillary membrane (pore radius and density ), this model enables the computation of transport parameters, such as hydraulic r k d k conductivity ( ), permeability ( ), and reflection coefficient ( ). Three types of ideally cylindrical pores, , are considered: L k 1 P k 2 σ k 3 k small pores , large pores and fenestrations . sp lp fWhere defines the capillary surface area, whereas the surface area of exchange, , is obtained from the size and density of the pore S pA k and from . is the perfusate viscosity. and are two decreasing hydrodynamic functions defined by Curry and recently k S p η F G[START_REF] Curry | Mechanics and thermodynamics of transcapillary exchange . Handbook of Physiology, The Cardivascular System[END_REF]

  exchanges are described by Kedem-Katchalsky equations , , which traduce fluid ( ), and solute ( ) [19 ] [20 ] exchanges for each type of pathway . Both mechanisms (solute and solvant) can be coupled, especially in the case of transport through k pores and fenestrations. Fluid movement, based on filtration through the vessel wall, is governed by Starling s hypothesis, and depends on ' hydrostatic ( ) and osmotic ( ) pressure gradients through the membrane ( ). p Π 4Where indexes and stand for the plasma and the space of Disse, respectively. Transport of molecules through the pores of the p d capillary wall, is described by ( ), derived from Patlak s equation and which accounts for concentration gradient through the wall. 5'[START_REF] Patlak | The flow of solute and sovent across a two-membrane system[END_REF]

  the past by Goresky, with the help of the Multiple Indicator Dilution (MID) technique .

  Wheredefines the surface area of the volume base, and enables the conversion from the fluid volume height to a pressure. The Base ρ choice for the base surface area controls changes in pressure: a small base induces high pressure variations whereas a large base reduces the pressure influence on the flux. Furthermore, high interstitial pressures , combined to low lymphatic pressures , may give rise to an evacuation of the interstial p d p l fluid (confined to the space of Disse), the lymphatic network, thus resulting in a lymphatic outflow ( ): via F l 12

  voxel is defined by three values describing the proton density PD, the spin-lattice relaxation time T and the spin-spin relaxation time T , 1 2 associated to the corresponding tissue. After application of a MRI sequence on the input, the magnetization kernel computes a set of RF signals, stored in the k-space. Addition of noise to the k-space, associated to its filtering before the reconstruction of the MR image using fast Fourier transform, provides relatively realistic images. Since each voxel of the liver model may contain both types of tissue (blood and parenchyma), the corresponding MRI parameters are computed knowing the proportion of each tissue in the voxel. The generation of new post-injection MRI maps (PD, T and T ) is then a necessary step to simulate DCE-MRI. gives the relationship between the of a voxel, and the associated contrast agent concentration. R 1,2 = T 1,2

  of contrast agent is first higher in the hepatic artery and then later in the portal vein). Finally, these 3D maps are used to simulate MR images via the simulator SIMRI. A typical T -weighted 3D Gradient Echo sequence is simulated at 3 Tesla 1 (repetition time ( ), 2.75ms, echo time ( ), 1.36ms, flip angle, 10 ). T R T E °Qualitative analysis A HCCp-like virtual tumor is first generated, based on a simulated adult liver. Results presented in correspond to images Fig. 5.a simulated after applying a T -weighted 3D Gradient Echo sequence on a 4mm-slice. Corresponding real images (obtained from the same 1 acquisition parameters) are presented in . At arterial phase, the amount of contrast agent is higher in the tumor than in normal Fig. 5.b tissue. Indeed, the arterial neovascularization, combined with an inhibition of the portal vascularization, produces an arterialization " "phenomena, known as a typical symptom in hepatic tumoral development. In consequence, the global concentration of Gd-DOTA in the tissue presents a peak that is close to the arterial one. The contrast between tumoral and normal tissues is then higher at this time, and T 1 and T values, driven by ( ), are decreased in the tumor. The T -weighted gradient echo sequence enhances the T contrast in images,

  .a associated MR images at the macroscopic one ( ). This time, though, a set of parameters (arterial, , and portal, , inflows, Fig. 10.b F ha F pv and permeability, reflected in the density of fenestrations ) is applied in correspondance to the five different stages of d f hepatocarcinogenesis. The values are given in . Exact values for each grade of the lesion remain undefined and vary according to table I

Fig. 1

 1 Fig. 1 Coupling of a physiological model of the liver and a model of dynamic MRI acquisition.

Fig. 2

 2 Fig. 2Main steps of the hepatocarcinogenesis. The best characterization of the tumor (benign nodule, malignant carcinoma) has to be reached in order to propose the most efficient therapy.

Fig. 3

 3 Fig. 3The hepatic vein is represented until the venule level (smallest vessels have a diameter around 30 50 micrometers). At the end of the growth simulation, the tree counts almost 50000 endings. Right: The three hepatic vascular trees (hepatic artery, portal vein, hepatic vein) are connected. For visualization purposes, the smallest vessels are not represented.

Fig. 4

 4 Fig. 4Axially-distributed BTEX model, describing the hepatic microvascular transport of molecules.

Fig. 5 Fig. 6 Fig. 7 Fig. 8

 5678 Fig. 5 Comparison between simulated a) and real b) Dynamic Contrast-Enhanced MR images of the liver. Acquisition conditions: 3D Gradient Echo sequence 3T ( 1.36 , 2.75 , 10 , 4 ); Contrast agent: Gd-DOTA; arterial phase. T E = ms T R = ms Flip Angle = °Slice Thickness = mm Up: portal phase. Note: arrows indicate the tumor. Bottom:

Fig. 9

 9 Fig. 9 Evaluation of the influence of the permeability, via the estimated permeability-surface product . Its effect on the model is observed at a) PS arterial and b) portal phases, at two different scales: at the macroscopic scale (left-hand graph), the mean signal is measured in two ROIs chosen in the tumoral and in the normal regions, after simulating T -weighted GRE MR images (slice thickness 4 mm); at the microscopic 1

Fig. 10

 10 Fig. 10 Evolution of a) Gd-DOTA concentration profile in the hepatic parenchyma and of b) T -weighted GRE MR images, with the grade of the 1

TABLE I

 I Some of the physiological parameters related to each stage of the hepatocarcinogenesis.

	lesion grade			NR	DNl	DNh	HCCm	HCCp
	( F ha ml.g 1 -min 1 . -	)	0.27	0.2	0.1	0.3	0.4
	( F pv ml.g 1 -min 1 . -	)	0.7	0.5	0.35	0.2	0
	( d f nb f en m / μ 2	)		11	11	5	0	0

[START_REF] Kretowski | Fast 3D modeling of vascular trees[END_REF]

[START_REF] Kretowski | Fast 3D modeling of vascular trees[END_REF]

[START_REF] Kretowski | Fast 3D modeling of vascular trees[END_REF]

[START_REF] Kretowski | Fast 3D modeling of vascular trees[END_REF]

[START_REF] Kretowski | Fast 3D modeling of vascular trees[END_REF]

[START_REF] Kretowski | Fast 3D modeling of vascular trees[END_REF]

[START_REF] Kretowski | Fast 3D modeling of vascular trees[END_REF]

[START_REF] Kretowski | Fast 3D modeling of vascular trees[END_REF]

[START_REF] Kretowski | Fast 3D modeling of vascular trees[END_REF]

[START_REF] Kretowski | Fast 3D modeling of vascular trees[END_REF]

[START_REF] Kretowski | Fast 3D modeling of vascular trees[END_REF]

[START_REF] Kretowski | Fast 3D modeling of vascular trees[END_REF]

[START_REF] Kretowski | Fast 3D modeling of vascular trees[END_REF]

Ackowledgements:

The authors are particularly grateful to Professor James B. Bassingthwaighte for his precious advices regarding the modeling of transport phenomena in biological systems, and to Professor Jean-Jacques Bellanger for his interest in numerical aspects of this work.