Anne-Laure Basdevant

Nathanaël Enriquez

Lucas Gerin

Distances in the highly supercritical percolation cluster

Keywords: first-passage percolation, supercritical percolation, TASEP

On the supercritical percolation cluster with parameter p, the distances between two distant points of the axis are asymptotically increased by a factor 1 + 1-p 2 + o(1 -p) with respect to the usual distance. The proof is based on an apparently new connection with the TASEP (totally asymmetric simple exclusion process).

Introduction

First passage percolation is a model introduced in the 60's by Hammersley and Welsh [START_REF] Hammersley | First-passage percolation, subadditive processes, stochastic networks, and generalized renewal theory[END_REF] which asks the question of the minimal distance D(x) between the origin 0 and a distant point x of Z 2 , when edges have i.i.d. positive finite lengths. One can prove by subadditivity arguments that in every direction such distances grow linearly: for each y ∈ Z 2 \{0}, D(ny)/n converges almost surely to a constant µ(y). The particular value of µ((1, 0)) =: µ is called the time constant. It is unknown except in the trivial case of deterministic edge lengths. We refer the reader to [START_REF] Kesten | Aspects of first passage percolation[END_REF] for an introduction on first passage percolation.

In this article we study the extreme case where edges have lengths 1 with probability p ∈ (0, 1), +∞ with probability 1 -p. Then, the distance D(x) coincides with the distance between the origin and x in the graph induced by bond percolation on Z 2 : each edge of Z 2 is open with probability p and closed with probability 1 -p. When p > 1/2, it is known (see [START_REF] Grimmett | Percolation[END_REF] Chap.1 for an introduction to bond percolation) that there exists almost surely a unique infinite connected component of open edges. We write x ↔ y if x, y belong to the same connected component, and x ↔ ∞ if x is in the infinite cluster.

Gärtner and Molchanov ([4] Lemma 2.8) were the first to rigorously prove that if 0 and x belong to the infinite component, D(x) is of order x. Garet and Marchand ( [START_REF] Garet | Asymptotic shape for the chemical distance and first passage percolation on the infinite Bernoulli cluster[END_REF], Th.3.2) improved this result and showed that, even if the subadditivity argument fails, the limit still holds: for each y ∈ Z 2 \{0}, there exists a constant µ(y) such that, on the event {0 ↔ ∞}, we have a.s. The aim of the present paper is to give an asymptotics of the time constant when p is close to one. One clearly has D(n, 0) ≥ n. On the other side, among the n edges of the segment joining 0 to (n, 0), about n(1 -p) of them are closed but with high probability the two extremities of each such edge can be joined by a path of length three. This naive approach suggests an upper bound of 1 + 2(1 -p) for the time constant. To our knowledge, the best known upper bound comes from Corollary 6.4 of [START_REF] Garet | Asymptotic shape for the chemical distance and first passage percolation on the infinite Bernoulli cluster[END_REF] and is equal to 1 + (1 -p). We obtain in the present paper a sharp asymptotics of the time constant when p goes to one. Theorem 1. On the event {0 ↔ ∞}, we have a.s.

1 µ p def = lim n→∞ 0↔(n,0) D(n, 0) n = 1 + 1 -p 2 + o(1 -p).
Our result says that the graph distance and the L 1 distance asymptotically differ by a factor 1 + (1 -p)/2. Note that Garet and Marchand ([3] Cor.6.4) observed that these two distances coincide in all the directions inside a cone containing the axis {y = x}. The angle of this cone is characterized by the asymptotic speed of oriented percolation of parameter p, studied by Durrett [START_REF] Durrett | Oriented percolation in two dimensions[END_REF].

The key ingredient of the proof relies on a correspondance between the synchronous totally asymmetric simple exclusion process (TASEP) on an interval and the graph distance on the percolation cluster inside an infinite strip.

First bounds on µ p

We denote by P p the product measure on the set of edges of Z 2 of length 1 under which each edge is open independently with probability p. Since p > 1/2, we have P p {0 ↔ ∞} > 0, so we can also define Pp the probability P p conditioned on the event {0 ↔ ∞}.

Pp {A} = P p {A ∩ {0 ↔ ∞}} P p {0 ↔ ∞} .
When no confusion is possible, we will omit the subscript p.

The origin is in the infinite cluster, unless there is a path of closed edges in the dual lattice surrounding 0. In the whole paper, we take p close enough to one so that this occurs with high probability (to fix ideas, with probability greater than 1 -2(1 -p) 4 ).

Because of the conditioning, it is not possible to apply directly subadditive arguments to the sequence D(n, 0). To overcome this difficulty, we adapt the ideas of [START_REF] Garet | Asymptotic shape for the chemical distance and first passage percolation on the infinite Bernoulli cluster[END_REF] and consider the sequence of points of the axis which lie in the infinite cluster. This enables us to derive bounds on µ p .

For n ≥ 1, let (T n (k), 0) be the k-th intersection of the infinite cluster and the set {(in, 0), i ∈ N}.

Proposition 1. We have for all n ≥ 1,

lim k→∞ E D(T 1 (k), 0) k 1l {0↔∞} = µ p ≤ E D(T n (1), 0) n 1l {0↔∞} .
1 We use notations introduced in [START_REF] Garet | Asymptotic shape for the chemical distance and first passage percolation on the infinite Bernoulli cluster[END_REF]: the subscript 0 ↔ (n, 0) means that we only take n's for which (n, 0) is in the infinite component. More precisely, if (T 1 (n), 0) stands for the n-th point of the half-line N × {0} belonging to the infinite component, then lim n→∞ 0↔(n,0)

D(n,0) n := limn→∞ D(T 1 (n),0) T 1 (n) .
Proof. This proposition is mainly a direct consequence of Lemma 3.1 of [START_REF] Garet | Asymptotic shape for the chemical distance and first passage percolation on the infinite Bernoulli cluster[END_REF]. Indeed, this lemma states that, for all n ≥ 1, there exists a constant f such that lim k→∞ D(T n (k), 0) nk = f P a.s. and in L 1 ( P).

Moreover, by subadditivity, we have nf ≤ E (D(T n (1), 0)| 0 ↔ ∞) and besides

µ p = P{0 ↔ ∞}f.
Combining these two facts, we get the upper bound. For the left equality, we now use the L 1 convergence in [START_REF] Durrett | Oriented percolation in two dimensions[END_REF] with n = 1. This gives

lim k→∞ E D(T 1 (k), 0) k 0 ↔ ∞ = f = µ p P{0 ↔ ∞} .
For p close to 1, the upper bound can be simplified using the following proposition.

Proposition 2. For all δ > 1, there exists c δ > 0 such that for all n ≥ 10 and p ∈ (5/6, 1), we have

E(D(T n (1), 0)1l {0↔∞} ) ≤ E(D(n, 0)1l {0↔(n,0)↔∞} ) + c δ (1 -p) 2 n δ .
To prove this proposition, we first show two lemmas.

Lemma 1. For all δ > 1, there exists a constant C δ > 0 such that, for n ≥ 1 and for p ∈ (5/6, 1), we have

E(D 2 (n, 0)1l {0↔(n,0)↔∞} ) ≤ (C δ n δ ) 2 .
Proof. We have

E(D 2 (n, 0)1l {0↔(n,0)↔∞} ) ≤ n 2δ + ∞ i=n 2δ P{D(n, 0) ≥ √ i, 0 ↔ (n, 0) ↔ ∞}.
Fix some q ∈ (0, 1 2δ ) such that q + 1 2δ ≤ 1 2 and for i ≥ n 2δ , let Γ i be the box [-i

1 2δ , n + i 1 2δ ] × [-i q 6 , i q 6 ]. A self avoiding path in Γ i has less than |Γ i | = i q 3 (n + 2i 1 2δ ) ≤ √ i steps. Thus P{D(n, 0) ≥ √ i, 0 ↔ (n, 0) ↔ ∞} ≤ P{0 (n, 0) in Γ i , 0 ↔ (n, 0) ↔ ∞}.
The event {0

(n, 0) in Γ i , 0 ↔ (n, 0) ↔ ∞} implies the existence in the dual of a path of closed edges starting from the border of Γ i , and having at least i q 3 steps since it must disconnect 0 and (n, 0). Using that

|∂Γ i | ≤ 6i 1 2δ , we get P{D(n, 0) ≥ √ i, 0 ↔ (n, 0) ↔ ∞} ≤ 6i 1 2δ (3(1 -p)) i q 3 .
This yields, for p ∈ (5/6, 1),

E(D 2 (n, 0)1l {0↔∞} ) ≤ n 2δ + 6 ∞ i=n 2δ i 1 2δ 1 2 i q 3 ≤ (C δ n δ ) 2 .
Lemma 2. Recall that T n (1) denotes the first point among {n, 2n, 3n, . . . } which is in the infinite component. There exists C > 0 such that for any p ∈ (5/6, 1), any n ≥ 10 and j ≥ 2

P {T n (1) ≥ jn} ≤ C(1 -p) 4+ √ j-2 . ( 2
)
Proof. For j = 2, the left-hand side in ( 2) is simply

P{(n, 0) is disconnected } ≤ C(1 -p) 4 . For 3 ≤ j ≤ 10, the event {T n (1) ≥ jn} is included in {(n, 0) ∞ and (2n, 0)
∞}. These two points are either disconnected by two different paths, or by the same path (which then has length ≥ 2n). The latter case has a much smaller probability, and then P{T n (1) ≥ jn} ≤ C ′ (1 -p) 8 , and (2) holds.

We now do the case j ≥ 10. For each integer i, let A i be the event that (ni, 0) is disconnected, and let

A i,r = {(ni, 0) is disconnected by a path included in [ni -r; ni + r] × [-r; r]} .
One plainly has that each A i,r ⊂ A i and that A i,r , A j,r are independent as soon as the corresponding boxes are disjoint, namely if 2r ≤ n|j -i|. For j ≥ 2, set J = ⌊ √ j⌋, we write

P{T n (1) ≥ jn} = P{A 1 , A 2 , . . . , A j-1 } ≤ P{A J , A 2J , A (J-1)×J } ≤ J-1 ℓ=1 P{A ℓJ,nJ/2 } + P{∃ℓ ≤ J -1; A J \ A ℓJ,nJ/2 } ≤ (1 -p) 4(J-1) + J i≥nJ i(3(1 -p)) i ≤ (1 -p) 4(J-1) + 4nJ 2 (3(1 -p)) nJ ≤ C(1 -p) 4+ √ j-2 ,
since a path disconnecting (nℓJ, 0) which is not included in the box [nℓJ -nJ 2 ; nℓJ + nJ 2 ] × [-nJ 2 ; nJ 2 ] has at least nJ edges. We are now able to prove Proposition 2.

Proof of Proposition 2. We write, using Lemmas 1 and 2:

E(D(T n (1), 0)1l {0↔∞} ) = ∞ j=1 E(D(jn, 0)1l {Tn(1)=jn} 1l {0↔∞} ) ≤ E(D(n, 0)1l {0↔(n,0)↔∞} ) + ∞ j=2 E(D(jn, 0) 2 1l {0↔(jn,0)↔∞} ) 1/2 P{T n (1) = jn} 1/2 ≤ E(D(n, 0)1l {0↔(n,0)↔∞} ) + ∞ j=2 C δ (jn) δ P{T n (1) ≥ jn} 1/2 ≤ E(D(n, 0)1l {0↔(n,0)↔∞} ) + C ′ δ n δ (1 -p) 2 ∞ j=2 j δ 1 6 (j-2) 1/4 ≤ E(D(n, 0)1l {0↔(n,0)↔∞} ) + c δ (1 -p) 2 n δ .

Percolation on a strip and TASEP

As a first step towards our main result, we shall study distances in percolation on an infinite strip. We will reduce this problem to the analysis of a finite particle system, this allows explicit computations from which will result the bounds in Theorem 1.

Here is the context we will deal with in the whole section. Fix an integer K and ε ∈ (0, 1). Let Z K be the infinite strip Z × [[-K, K]], with three kinds of edges :

• Vertical edges {(i, j) → (i, j + 1), i ∈ Z, j ∈ [[-K, K -1]]}; • Horizontal edges {(i, j) → (i + 1, j), i ∈ Z, j ∈ [[-K, K]]}; • Diagonal edges {(i, j) → (i + 1, j + 1) and (i, j) → (i + 1, j -1)}.
We now consider a random subgraph of Z K equipped with distances: Cross Model . (i) Vertical and horizontal edges have length 1, whereas diagonal edges have length 2.

(ii) Diagonal and vertical edges are open.

(iii) Each horizontal edge is open (resp. closed) independently with probability 1 -ε (resp. ε).

For i ≥ 0 and j ∈ [[-K, K]], let D K,d (i, j) be the distance between (0, 0) and (i, j) inside Z K in the Cross Model (see an example in Fig. 1), the 'd' stands for the addition of diagonal edges. Sometimes we also need to consider the distance between two vertices x, y of Z K , which will be denoted by D K,d (x → y).

Since vertical and diagonal edges are open, every point in Z K is connected in the Cross Model to 0, hence D K,d (i, j) is finite for every i, j. We also set The main goal of this section is to estimate D K,d (n, 0), when K, n are large. To do so, we introduce a particle system associated to the process (D K,d i ) i≥0 . Let us consider the state space {•, •} 2K (identified to {1, 0} 2K ), and denote its elements in the form (y -K+1 , y -K+2 , . . . , y 0 , y 1 , . . . , y K ).

D K,d i = D K,d (i, j), j ∈ [[-K, K]] . By construction, we have along each vertical edge |D K,d (i, j) -D K,d (i, j + 1)| = 1.
Let (Y i ) i≥0 be the process with values in {•, •} 2K defined as follows :

∀j ∈ [[-K + 1, K]], Y j i = • = 1 if D K,d (i, j) = D K,d (i, j -1) -1. • = 0 if D K,d (i, j) = D K,d (i, j -1) + 1.
Let say that the site j is occupied by a particle at time i if Y j i = • and empty otherwise. • at time t + 1, for each j, a particle at position j = -K + 1, . . . , K -1 moves one step forward if the site j + 1 is empty at time t, with probability α and independently from the other particles.

• at time t+1, a particle enters the system at position -K +1 if site -K +1 is empty at time t, with probability β.

• at time t + 1, if there were a particle at position K at time t, it exits the system with probability γ. Proof of Proposition 3. Let us note that since all the vertical edges are open, the optimal path from 0 to (i, j) ∈ Z K , i ≥ 0 never does a step from right to left. Moreover, the vector D K,d i+1 depends only on D K,d i and on the edges

{(i, j) → (i + 1, j), j ∈ [[-K, K]]}, hence it is Markov.
Let us now prove that the displacement of particles follows the rules of TASEP. We detail the case in which there is a particle on the edge (i, j -1) → (i, j) but no particle on the edge (i, j) → (i, j + 1), i.e. at time i there is a particle in position j and no particle in position j + 1.

This means that if

D K,d (i, j) = ℓ then D K,d (i, j -1) = D K,d (i, j + 1) = ℓ + 1.
whether the horizontal edges (i, j -1) → (i + 1, j -1) and (i, j + 1) Now, D K,d (i + 1, j) depends only on the edge (i, j) → (i + 1, j): it is equal to ℓ + 1 if this edge is open and the particle lying at j stays put. If the edge (i, j) → (i + 1, j) is closed, D K,d (i + 1, j) = ℓ + 3, which corresponds for the particle lying at j to a move to j + 1.

→ (i + 1, j + 1) are open or closed, we have D K,d (i + 1, j -1) = D K,d (i + 1, j + 1) = ℓ + 2 (this
We leave the cases in which a particle is followed by an other particle and in which an empty edge is followed by an empty edge, which are similar, to the reader.

We do the bottom boundary case (when a particle may enter the system, see the figure below). Assume that there is no particle at time i in position -K + 1. This means that if we set ℓ = D(i, -K) then D(i,

-K + 1) = ℓ + 1. If the edge (i, -K) → (i + 1, -K) is open then D(i + 1, -K) = ℓ + 1
and there is no particle at time i + 1 in position -K + 1. Otherwise D(i + 1, -K) = ℓ + 3 and a particle appears at time i + 1 in -K + 1. The right-boundary case (when a particle exits) is similar.

We have thus seen in the proof that the knowledge of (Y i ) i≥0 fully determines the metric D K,d using the following recursive identity:

D K,d (i + 1, j) = D K,d (i, j) + 1 + 2.1l {Y j i =•,Y j+1 i =•,Y j i+1 =•,Y j+1 i+1 =•} . Let ν K,ε be the stationary measure of the synchronous TASEP (Y i ) i≥0 . Let ν K,ε (•, •) def = ν K,ε (y 0 = •, y 1 = •).
Proposition 4. We have the following asymptotics for the distances on Z K in the Cross Model:

lim n→∞ 1 n E(D K,d (n, 0)) = 1 + 2εν K,ε (•, •).
Proof of Proposition 4. We have seen in the proof of Proposition 3 that

D K,d (i + 1, 0) = D K,d (i, 0) + 1,
unless a particle has moved at time i from position 0 to 1, in which case D K,d (i+ 1, 0) = D K,d (i, 0) + 3. This shows that

E(D K,d (i+1, 0)) = E(D K,d (i, 0))+1+2P Y 0 i = •, Y 1 i = •, {(i, 0) → (i + 1, 0) is closed} , thus E(D K,d (n, 0)) = n + 2ε n-1 j=0 P Y 0 j = •, Y 1 j = • ,
which gives, by the Markov chain ergodic theorem, the proof of the proposition.

We thus need an estimate of ν K,ε (•, •). It turns out that the stationnary measure of the synchronous TASEP was studied in detail by Evans, Rajewsky and Speer [START_REF] Evans | Exact solution of a cellular automaton for traffic[END_REF] using a matrix ansatz. Proposition 5. The following identities relative to ν K,ε (•, •) hold:

(i) For all K, ε, ν K,ε (•, •) = A ε (K) εA ε (K) + A ε (K + 1)
,

with A ε (K) = 1 K K k=1 K k K k+1 (1 -ε) k . (ii) lim K→∞ ν K,ε (•, •) = 1 - √ 1 -ε 2ε . (iii) For all α > 0, lim ε→0 ν ε -α ,ε (•, •) = 1 4 .
Proof. The first two assertions are consequences of (4.24), (8.21) and (10.13) of [START_REF] Evans | Exact solution of a cellular automaton for traffic[END_REF]. For (iii), we are led to prove that lim ε→0

Aε(ε -α +1)

Aε(ǫ -α ) = 4. Denote a(K, k) def = K k K k + 1 (1 -ε) k .
The ratio

a(K, k + 1) a(K, k) = (1 -ε) (K -k)(K -k -1) (k + 1)(k + 2) (3) is asymptotically equal to (1 -ε)( K-k k ) 2
for large values of k and K -k. Therefore the sequence (a(K, k)) 1≤k≤K increases from 1 to some k max (K) and decreases after. Moreover, k max (K) ∼ K( 12 -ε 8 ). For K > 0 and β ∈ (0, 1), let us decompose the sum KA ε (K) into

K( 1 2 -ε β )-1 k=1 a(K, k) + K 2 k=K( 1 2 -ε β ) a(K, k) + K( 1 2 +ε β ) k= K 2 +1 a(K, k) + K k=K( 1 2 +ε β )+1 a(K, k) and denote by A - β (K), B - β (K), B + β (K), A + β (K)
, the four successive sums. Note first, from the analysis of the ratio (3), that A - β (K) and A + β (K) are both subgeometric sums with rate 1 -Cε β . This implies,

A + β (K) ≤ Cε -β a(K, K( 1 2 + ε β )) A - β (K) ≤ Cε -β a(K, K( 1 2 -ε β )).
In the same time, from the variations of the sequence (a(K, k)) 1≤k≤K , we deduce

B + β (K) ≥ Kε β a(K, K( 1 2 + ε β )) B - β (K) ≥ Kε β a(K, K( 1 2 -ε β )).
Therefore, as soon as the K ≫ ε -2β , the ratio

KAε(K) B - β (K)+B + β (K) goes to 1. Now, for all k, a(K + 1, k) a(K, k) = (K + 1) 2 (K -k)(K -k + 1)
.

This ratio is between 4(1 -cε β ) and 4(1

+ Cε β ) when k is between K( 1 2 -ε β ) and K( 1 2 + ε β
). This implies that, for all β > 0, the ratio

B - β (K+1)+B + β (K+1) B - β (K)+B + β (K)
converges to 4 as ε goes to 0 and K goes to ∞. Consequently, the same occurs also for (K+1)Aε(K+1)

KAε(K)
, as soon as K ≫ ε -2β . We choose β = min 1, α 4 to prove the result.

1. i is maximal, 2. among them, j is minimal.

In this way, there is no vertical edges neither to the right of E nor below it. There are three cases according to whether the other edge joining (i, j -1) in P goes to (i -1, j -1), to (i -1, j) or to (i -1, j -2). In each of these three cases we can do a substitution that removes the vertical edge or moves it to the left:

(i,j) (i,

j) (i,j) (i,j) (i,j) (i,j-1) (i,j) (i,j-1) (i,j-1) (i,j-1) (i,j-1) (i,j-1)

This substitution does not change the length of the path, it is thus still optimal. After iterating the process, there may remain vertical edges only on the first column.

Proposition 7. Consider a standard percolation on Z K where each horizontal and vertical edge is open with probability

1 -ε. Denote D K the associated distance. Let A = A(K, n, ε) be the event "in each square of area 1 of [[0; n]] × [[-K; K]] i.e.
of vertices {(i, j), (i, j + 1), (i + 1, j + 1), (i, j + 1)}, at most one edge is closed". Then, for each n, K,

E(D K (n, 0) | A) ≤ E(D K,d (n, 0)) + 3K.
Proof. On the event A, adding diagonal edges with length 2 does not decrease the length of optimal paths since either the path (i, j) → (i, j +1) → (i+1, j +1) or (i, j) → (i + 1, j) → (i + 1, j + 1) is open and moreover, there is at most a distance 3 between the two extremities of any edge of [[0

; n]] × [[-K; K]].
Thanks to Lemma 3, there is always for the Cross Model a path of minimal length using only horizontal and diagonal edges except maybe on the at most K vertical steps along the first column {0} × [[-K, K]]. Thus, on the event A, we get

D K (n, 0) ≤ D K,d (n, 0) + 3K.
This yields

E(D K (n, 0) | A) ≤ E(D K,d (n, 0) | A) + 3K.
Let us note now that A is an increasing event and -D K,d (n, 0) is an increasing random variable. Thus, the FKG inequality (see [START_REF] Grimmett | Percolation[END_REF]) yields

E(-D K,d (n, 0)1 A ) ≥ E(-D K,d (n, 0))P{A}
which can be rewritten

E(D K,d (n, 0) | A) ≤ E(D K,d (n, 0)).

The lower bound

We now return to the original model and consider percolation on Z 2 where each edge is closed independently with probability ε = 1 -p. Recall that D(i, j) denotes the distance between the origin and (i, j) ∈ Z 2 . Using Proposition 1, in order to prove the lower bound in Theorem 1, we just need to show that

lim k→∞ 1 k E(D(T 1 (k), 0)1l {0↔∞} ) ≥ 1 + ε 2 + o(ε).
Let D d (x) be the distance between the origin and x in Z 2 , with diagonal and vertical edges all open. Adding edges decreases the distances: D d (x) is always smaller than D(x). Thus, we have,

E(D(T 1 (k), 0)1l {0↔∞} ) ≥ E(D d (T 1 (k), 0)1l {0↔∞} ) ≥ E(D d (k, 0)1l {0↔∞} ),
since T 1 (k) ≥ k and the sequence (D d (n, 0)) is increasing (thanks to vertical edges). Note that every minimal path from the origin to (k, 0) (with diagonal and vertical edges open) has a length less than 2k and lies in the strip Z k . Therefore

E(D d (k, 0)1l {0↔∞} ) = E(D k,d (k, 0)1l {0↔∞} ). Using Proposition 6, we have E(D k,d (k, 0)) ≥ k(1 + 2εν k,ε (•, •)). Besides, E(D k,d (k, 0)(1 -1l {0↔∞} )) ≤ 2kP{0 ∞} ≤ 4kε 4 .
We get

E(D(T 1 (k), 0)1l {0↔∞} ) ≥ k(1 + 2εν k,ε (•, •) -4ε 4 ).
Letting k tends to infinity, we obtain with Proposition 5 the desired lower bound.

The upper bound: a short path

As in the previous section, we consider standard percolation on Z 2 where each edge is closed with probability ε = 1 -p. Recall that Proposition 2 states that for any δ > 1 there exists a c δ > 0 such that, for all ε small enough and n ≥ 10,

µ p ≤ E(D(n, 0)1l {0↔(n,0)↔∞} ) n + c δ ε 2 n δ-1 . (4)
So to prove Theorem 1, it is sufficient to find n = n(ε) ≥ 10 such that εn δ-1 = o(1) and

E(D(n, 0)1l {0↔(n,0)↔∞} ) ≤ n(1 + ε 2 + o(ε)).
Fix ℓ > 1 and K, n > 2ℓ and define the boxes (C i ) i≥0 by

C i = [ℓ, n -ℓ] × [(2i -1)K + 1, (2i + 1)K].
Let E i denote the event "in each square of area 1 of C i i.e. of vertices {(i, j), (i, j+ 1), (i + 1, j + 1), (i, j + 1)}, at most one edge is closed". Let us note that the events (E i ) i≥0 are independent and have the same probability. Moreover, we have

P {E c i } ≤ |C i |P{more than 2 edges are closed in a given square} ≤ 2Kn(6ε 2 + 4ε 3 + ε 4 ) ≤ 22Knε 2 . Hence, if I = inf{i ≥ 0, E i }, I + 1 is geometrically distributed with parameter P{E 0 }. C 0 C i K 0 ε i ε ′ i n c i c ′ i ℓ Figure 3:
The sketch of an almost optimal path, on the event B i ∩ {I = i}.

Let c i , c ′ i denote the vertical boundaries of C i i.e. c i def = {ℓ} × [(2i -1)K + 1, (2i + 1)K] and c ′ i def = {n -ℓ} × [(2i -1)K + 1, (2i + 1)K] and ε i , ε ′ i the boxes defined by ε i = [-ℓ, ℓ] × [-K, (2i + 1)K], ε ′ i = [n -ℓ, n + ℓ] × [-K, (2i + 1)K].
Let B i be the event

B i = {0 ↔ c i in ε i } ∩ {(n, 0) ↔ c ′ i in ε ′ i }.
We have

E(D(n, 0)1l {0↔(n,0)↔∞} ) ≤ ∞ i=0 E(1l {I=i}∩Bi D(n, 0))+ ∞ i=0
E(D(n, 0)1l {I=i}∩B c i ∩{0↔(n,0)↔∞} ).

(5) Let χ 1 (i) (resp. χ ′ 1 (i)) denote the distance between 0 and the segment c i inside the box ε i (resp. (n, 0) and the segment c ′ i inside the box ε ′ i ). Denote also χ 2 (i) the distance between (ℓ, 2iK) and (n -ℓ, 2iK) inside the box C i . We have D(n, 0) ≤ χ 1 (i)+χ ′ 1 (i)+χ 2 (i)+max Besides, we have

P{B c i ∩ {I = i} ∩ {0 ↔ (n, 0) ↔ ∞}} ≤ 2P{{0 c i in ε i } ∩ {0 ↔ ∂ε i } ∩ {I = i}} = 2P{{0 c i in ε i } ∩ {0 ↔ ∂ε i }}P{I = i},
since the two first events only depend on the edges inside the box ε i and the event {I = i} only depends on the edges inside the box C i . If 0 ↔ ∂ε i but is not connected to c i in the box ε i , then there exists in the dual graph a path of closed edges with at least 2ℓ edges. Hence P{{0 c i in ε i } ∩ {0 ↔ ∂ε i }} ≤ 2K(i + 1)2ℓ(3ε) 2ℓ .

Figure 1 :

 1 Figure 1: A configuration of percolation in Z K for K = 2 and the associated distances D K,d . Note the importance of diagonal edges: D K,d (3, 0) = 5 instead of 7 if there was none.

Figure 2 :Definition 1 .

 21 Figure 2: The same configuration with particles.

Proposition 3 (

 3 The particles follow a TASEP). The processes (D K,d i ) i≥0 and (Y i ) i≥0 are Markov chains. Moreover, (Y i ) i≥0 has the law of discrete time synchronous TASEP on [[-K + 1, K]] with jump rate ε and exit and entry rate ε.

  is because the two diagonal edges starting from (i, j) are open), see the left part of the figure below:

x∈ciD

  (x → (ℓ, 2iK))+max y∈c ′ i D(y → (n-ℓ, 2iK)).

( 6 )

 6 Moreover, on the event {I = i} ∩ B i , the r.h.s. of this inequality is finite and we haveχ 1 (i) ≤ |ε i | = 4ℓ(i + 1)K χ ′ 1 (i) ≤ |ε ′ i | = 4ℓ(i + 1)K (7) max x∈ci D(x → (ℓ, 2iK)) ≤ 3K max y∈c ′ i D(y → (n -ℓ, 2iK)) ≤ 3K. (8)To bound χ 2 (i), let us note that χ 2 (i) only depends on the edges inside the box C i . Thus, we haveE(χ 2 (i)1l {I=i}∩Bi ) ≤ E(χ 2 (i)1l {I=i} ) = E(χ 2 (i) | I = i)P{I = i} = E(χ 2 (0) | I = 0)P{I = i}.Moreover, the event {I = 0} coincides with the event A(K, n -2ℓ, ε) defined in Proposition 7. This yields, using Proposition 6 and Proposition 7,E(χ 2 (0) | I = 0) ≤ E(D K,d (n -2ℓ, 0)) + 3K ≤ n(1 + 2εν K,ε (•, •)) + 5K. (9)Combining (6),(7),(8),(9), we get∞ i=0 E(1l {I=i}∩Bi D(n, 0)) ≤ n(1 + 2εν K,ε (•, •)) + 5K + K(8ℓE(I + 1) + 6). (10)For the second term in the right hand side of (5), using Hölder inequality, we getE(D(n, 0)1l {I=i}∩B c i ∩{0↔(n,0)↔∞} ) ≤ P{B c i ∩ {I = i} ∩ {0 ↔ (n, 0) ↔ ∞}} 1/2 E(D 2 (n, 0)1l {0↔(n,0)↔∞} ) 1/2 . (11)Lemma 1 yields E(D 2 (n, 0)1l {0↔(n,0)↔∞} ) 1/2 ≤ C δ n δ .

We will need further the following bounds on E(D K,d (n, 0)). Proposition 6. For n ≥ 0, we have

Proof. By subadditivity of (E(D K,d (n, 0))) n≥0 , the sequence E(D K,d (n, 0))/n is decreasing. Together with Proposition 4, this proves the left inequality.

For the right inequality, the idea is to start the Markov chain (Y i ) i≥0 from its stationary distribution. Let Ỹ0 = ( Ỹ j 0 , j ∈ [[-K + 1, K]]) with law ν K,ε . Take DK,d (0, 0) = 0 and define inductively DK,d (0, j) for each j such that DK,d (0, j) = DK,d (0, j -1) -

) is a realization of the chain starting from DK,d 0 then, by Proposition 3, we have

By construction of the chain, there is a (random) J such that

where D K,d is, as before, the true distance after percolation in Z K . Using triangular inequality, we get

Taking expectation gives the expected bound.

If we gather the results obtained so far in this section, we roughly obtain that for the Cross Model,

In order to apply this result to usual percolation on strips (with horizontal and vertical edges open with probability 1 -ε), we have to prove that distances in both models differ very little. Fortunately, it happens that distances in both models are quite similar as soon as there are no contiguous closed edges, which is the case with high probability in any fixed rectangle, when ε is small. We first choose in the Cross Model a particular path among all the minimal paths: Lemma 3. For the Cross Model, for each K, n, there exists a path from (0, 0) to (n, 0) of minimal length that only goes through horizontal and diagonal edges, except possibly through the vertical edges of the first column

Proof. Consider one of the minimal paths from (0, 0) to (n, 0), we will modify this path P to get one which satisfies the desired property. Among all the vertical edges of P, let E = (i, j -1) → (i, j) be the one for which Plugging this into (11) gives

(12) Recall now that I + 1 is a geometric random variable with parameter P{E 0 } ≥ 1 -22Knε 2 . Thus, there exists C < ∞ such that, for any K, n, ε are such that 22Knε 2 ≤ 1 2 , we have

Combining (4), ( 5), ( 10) and (12), we get, for ε ∈ (0, 1 6 ), K, n ≥ 2ℓ and such that 22Knε 2 ≤ 1 2

where C δ,ℓ is a constant depending only on ℓ and δ > 1. Take now ℓ = 3, n = ε -3/2 , K = ε -1/4 and δ = 3/2. We get

We conclude by using Proposition 5.


Msc G Guérin

PhD D Ambard

Professor Pascal Swider
email: swider@cict.fr

Cells, growth factors and bioactive surface properties in a mechanobiological model of implant healing

Keywords: implant fixation, mechanobiology, growth factor, osteoblast, transport in porous media, biomaterial implant fixation, mechanobiology, growth factor, osteoblast, transport in porous media

Interface conditions are of prime importance for implant fixation in the early post-operative period and modelling of specific biochemical interactions at implant surface is still missing. We hypothesized that updating osteoblast adhesion properties and growth factor source in an active zone located at the implant surface was relevant to model biochemical interactions of implant with its environment. We proposed an innovative set of diffusive-convective-reactive equations which relevant parameters were the cell decay factor, the cell motility, and the growth factor balance.

Initial comparison with histomorphometic results from a stable PMMA canine implant model provided an encouraging base to implement a numerical sensitivity analysis to evaluate the role of three types of bioactive surfaces: acid-etched titanium, coarse grit blasted acidetched titanium and coarse grit blasted acid-etched titanium with RGDS peptide. We found that cell diffusion decrease (acid-etched + RGDS peptide vs. PMMA), and increase of local growth factor fraction (PMMA vs. acid-etched + RGDS peptide), significantly improved the amount of mineralized tissue on the implant surface. When the variation of structural fraction to cell motility and growth factor synthesis was investigated, an envelope pattern with an optimum was obtained but this could be exceeded for strong surface modifications and/or for high growth factor concentrations. The model also confirmed that implant bioactive properties should play a limited role to reduce heterogeneity of new-formed tissue. In conclusion, we suggested that our innovative theoretical approach was relevant to investigate implant fixation and could potentially help in reduction of implant revision.

Introduction

Conditions favoring periprosthetic tissue healing in the early post-operative period is an open clinical problem. Low performance of implant fixation is generally associated with a low mineralization or a heterogeneous ossification of the new-formed tissue. Three parameters have a significant influence on this process: the surgical technique (Hahn et al. 1988;Morshed et al. 2007), the post-operative mechanical stimuli and the biochemical factors amongst which physico-chemical properties of implant surface are essential [START_REF] Anderson | Biological responses to materials[END_REF]Colnot et al.2007;Belmont et al. 2008;Kromer et al. 2008). Clinical studies and animal studies investigated the role of implant surface properties in osteointegration [START_REF] Overgaard | Calcium phosphate coatings for fixation of bone implants[END_REF][START_REF] Svehla | Morphometric and mechanical evaluation of titanium implant integration: Comparison of five surface structures[END_REF]Buser et al., 2007;[START_REF] Duyck | Effect of intermittent loading and surface roughness on peri-implant bone formation in a bone chamber model[END_REF][START_REF] Qu | The initial attachment and subsequent behavior regulation of osteoblasts by dental implant surface modification[END_REF]Schwartz et al., 2007).

Cell adhesion involves short-term events such as physicochemical linkages between cells and substrates and long-term events leading to the regulation of cell expression involving extracellular proteins, cell membrane proteins and cytoskeleton proteins [START_REF] Anselme | Osteoblast adhesion on biomaterials[END_REF]. Cytoskeleton proteins mediate interaction between cell membrane proteins and nucleus via actin filaments, and regulate gene expression [START_REF] Sinha | Regulation of human osteoblast integrin expression by orthopaedic implant materials[END_REF]. Extracellular matrix proteins have adhesion properties due to RGD sequence (arginine-glycine-aspartate) specific to fixation of cell membrane receptors [START_REF] Grzesik | Bone matrix RGD glycoproteins: immunolocalization and interaction with human primary osteoblastic bone cells in vitro[END_REF]. Integrin type receptors of cell membrane proteins induce cell-substrate adhesion establishing focal contacts. Cell migration on substrate can be explained by adhesive interaction involving integrins that locally generate cell traction (Wu et al., 1995). In case of implants, material properties and surface modification clearly affect cell adhesion and proliferation, and growth factors production [START_REF] Kieswetter | Surface roughness modulates the local production of growth factors and cytokines by osteoblast-like MG-63 cells[END_REF][START_REF] Kilpadi | Hydroxylapatite binds more serum proteins, purified integrins, and osteoblast precursor cells than titanium or steel[END_REF][START_REF] Colnot | Molecular analysis of healing at a bone-implant interface[END_REF][START_REF] Popat | Influence of engineered titania nanotubular surfaces on bone cells[END_REF][START_REF] Anderson | Biological responses to materials[END_REF]Raush-fan et al., 2008). It is observed that sensitivity to surface topology is significant [START_REF] Pessková | The influence of implant surface properties on cell adhesion and proliferation[END_REF][START_REF] Neidlinger-Wilke | Cyclic stretching of human osteoblasts affects proliferation and metabolism: a new experimental method and its application[END_REF][START_REF] Fermor | Primary human osteoblast proliferation and prostaglandin E 2 release in response to mechanical strain in vitro[END_REF] and generally, osteoblastic proliferation decreases when surface roughness increases [START_REF] Martin | Effect of titanium surface roughness on proliferation, differentiation, and protein synthesis of human osteoblast-like cells (MG63)[END_REF][START_REF] Rosa | Rat bone marrow cell response to titanium and titanium alloy with different surface roughness[END_REF][START_REF] Mustafa | Determining optimal surface roughness of TiO(2) blasted titanium implant material for attachment, proliferation and differentiation of cells derived from human mandibular alveolar bone[END_REF].

Cells produce a significant range of growth factors amongst with TGF-β1 is the major regulator of bone formation. It involves unique mitogenic action on osteoblasts and it shows ability to stimulate matrix formation via autocrine, paracrine and endocrine modes [START_REF] Roberts | Transforming growth factor-β[END_REF]. Here also, substrate properties modify the amount of TGF-β1 produced by osteoblasts. It is shown that the production rate is higher for rough surface [START_REF] Schwartz | Underlying mechanisms at the bone-surface interface during regeneration[END_REF] or in the presence of oxide (Albrektsson et al., 1983), while the presence of metal particles is affecting TGF-β1 synthesis (Wang et al., 1997).

Relevant biomathematical models of migration and differentiation have been proposed in skin and fracture healing, angiogenesis or bioreactors (Tracqui, 1999;Tranqui et al., 2000;Bailon-Plaza and Van Der Meulen, 2001;Sengers et al., 2007). Most computational models of implant fixation were initially based on pure mechanical approaches [START_REF] Ramamurti | Factors influencing stability at the interface between a porous surface and cancellous bone: a finite element analysis of a canine in vivo micromotion experiment[END_REF][START_REF] Prendergast | Finite element models in tissue mechanics and orthopaedic implant design[END_REF], Viceconti et al., 2000;[START_REF] Mesfar | Modeling of biomedical interfaces with nonlinear friction properties[END_REF][START_REF] Simon | Influence of the stiffness of bone defect implants on the mechanical conditions at the interface; a finite element analysis with contact[END_REF]2003) and more recently, biochemical aspects have been added [START_REF] Fernandes | A contact model with ingrowth control for bone remodelling around cementless stems[END_REF][START_REF] Puthumanapully | Do multi-layer beads on porous coated implants influence bone ingrowth? A finite element study[END_REF], Geris et al. 2008a & b;[START_REF] Checa | A mechanobiological model for tissue differentiation that includes angiogenesis: a lattice-based modeling approach[END_REF]. Reliability of numerical predictions has rarely been correlated to in-vivo or ex-vivo data. We proposed an original framework to model implant healing unifying mechanical and biomathematical approaches while considering the biological tissue as a multiphasic convective-diffusive-reactive medium (Ambard et al., 2004a&b, 2005a&b, 2006(Ambard et al., 2004a&b, 2005a&b, , 2007(Ambard et al., 2004a&b, 2005a&b, , 2009)). The reliability of our theoretical approach was supported by ex-vivo data from canine implant models [START_REF] Søballe | Tissue ingrowth into titanium and hydroxyapatite-coated implants during stable and unstable mechanical conditions[END_REF][START_REF] Vestermark | Mechanical interface conditions affect morphology and cellular activity of sclerotic bone rims forming around[END_REF]).

It appears that interface conditions are of prime importance for bone apposition and implant anchorage. However, modelling of specific biochemical interactions at implant surface is still missing and in our knowledge, no predictive model explicitly involved the biochemical interactions between the implant surface properties and the surrounding biological tissue. We hypothesized that local modification of cell adhesion properties and growth factor source was a way to model interactions between the implant and its environment. To proceed, we firstly proposed a set of new mechanobiological governing equations. Second, we evaluated the model relevance against histologic data from a PMMA canine implant model [START_REF] Vestermark | Mechanical interface conditions affect morphology and cellular activity of sclerotic bone rims forming around[END_REF]. Around this reference solution, we thirdly implemented a numerical sensitivity analysis to quantified the potential role of three types of modified surfaces: acid-etched titanium, coarse grit blasted acid-etched titanium and coarse grit blasted acid-etched titanium with RGDS peptide.

Material & Methods

Governing equations:

The association of porous media mechanics to biomathematics allowed the diffusiveconvective-reactive governing equations (1) to be derived. Non-linear coefficients L, C, D, Ω described the local variation, the convection, the diffusion, and the source terms. Output measure u was the structural or mineralized fraction φ s , the fluid fraction φ f , the growth factor concentration C g (TGF-β1) and the osteoblast concentration C o . At the scale of tissue volume unit, volumes of cell phase and growth factor phase were neglected compare to the volume of structural and fluid phases. Fluid fraction φ f also described the porosity of the newly formed tissue since φ s + φ f = 1. Phases were incompressible with no substrate strain in quasi-static isothermal evolving. Structural porosity, fluid flow and growth factors conditioned the cellular behavior: proliferation, chemotactic and haptotactic migrations of cells, and mineral fraction apposition.

u L D u C u t Ω ∂ = Δ -∇ + ∂ (1)
The evolving of osteoblast concentration was described by cell conservative equation (2) where diffusion coefficient D o and convection coefficients χ o and h o , respectively allowed modelling random cell migration [START_REF] Lauffenburger | Measurement of leukocyte and chemotaxis parameters with a linear under-agarose migration assay[END_REF][START_REF] Dee | Osteoblast population migration characteristics on substrates modified with immobilized adhesive peptides[END_REF], chemotactic active migration and haptotactic active migration. Chemotactic flux of cells was dependant upon the growth factor gradient g C ∇ whereas haptotactic flux was conditioned by the gradient of structural matrix s φ ∇ . The cell source Ω o represented the proliferation process with the logistic law (3) where α o and C om respectively were the rate of cell divisions and the proliferation threshold. The role of available growth factors in cell proliferation was specifically taken into account using algebraic difference between the current concentration C g and the physiological threshold C gp [START_REF] Centrella | Transforming growth factor-β gene family members and bone[END_REF]Alliston and Derynck, 2000).

( ) ( ) ( ) ο φ φ χ φ φ φ Ω ∂ = ∇ -∇ - ∇ + ∂ f o o f o o f g o f o s C D C C h C t (2) 3 ( ) ( Ω α φ = - × o o f o om o g gp C C C C C ) - (3)
The growth factor concentration per tissue volume unit was expressed by conservative equation (4). Replacing the initial diffusion coefficient D g by the effective diffusion coefficient [START_REF] Maxwell | A treatise on electricity and magnetism[END_REF][START_REF] Neale | Prediction of transport processes within porous media: Diffusive flow processes within an homogeneous swarm of spherical particles[END_REF] allowed modelling local retention of growth factors in the primary gap due to surgery. At Darcy scale, fluid flux relative to structural phase was relied to the pore pressure p and permeability κ s dependant from φ f [START_REF] Cowin | Bone mechanics handbook[END_REF]. ( ) ( )

φ φ φ κ Ω ∂ = ∇ - ∇ + ∂ f g ge f g f g s C D C C p t g with 1 2
(3 )

φ φ - = - ge f g f D D (4)
Equation ( 5) governed the growth factor source Ω g . The endocrine mode was governed by the growth factor concentration imposed as initial condition in the model. To account modes of production, equation ( 5) involved fluid fraction (or tissue porosity) φ f , osteoblast concentration C o , and growth factor concentrations (C g , C gp ). Coefficient α g was the coefficient of growth factor synthesis and d g was the decay factor taking into account the consumption of growth factor due to osteoblast activity [START_REF] Stein | Mechanisms of action of skeletal Growth factors in osteoblasts[END_REF].

2

(

Ω α φ φ = - × - ) g g f o g f o g gp C d C C C (5)
The synthesis of mineral matrix was modelled in equation ( 6) using the source Ω s involving the synergetic action of osteoblast C o and growth factors C g [START_REF] Meinel | Localized insulin-like growth factors I delivery to enhance new bone formation[END_REF]. It was promoted by the fluid fraction φ f and the coefficient of matrix synthesis α s .

2

(

Ω α φ = × - ) s s f o g gp C C C (6)

Active zone of implant

The periprosthetic zone of interest was divided into three sub-domains as shown in Figure 1a. The peripheral domain was the host trabecular bone; the boundaries of intermediate domain were the implant surface and the drill hole. This corresponded to the immediate postoperative gap. We considered a third domain corresponding to the zone of influence of the implant or the active zone. This was a layer of a few microns thickness located at the implant surface, where implant bioactive properties modified the growth factor source and the cell flux. Cell adhesion was surface-dependant thanks to diffusion coefficient D o in equation ( 2). Low value of D o meant significant adhesive properties of the surface.

The local influence on growth factor phase was managed adapting conservative equations ( 5) and ( 6). First, the decay factor d g , expressed by equation 7a, was obtained assuming a negligible production mode (α g =0), a negligible pressure gradient ( 0) ∇ = p , and nil growth factor gradients ( . Second, the synthesis factor α 0) ∇ = g C g was calculated using equation (7b) with d g previously determined and according to the same assumptions in gradients.

( ) [ ( ) φ - ∂ = - × - ∂ g g f o g C d C C t ] g p C (a) 1 [ ( α φ ) ] ∂ = + × - ∂ g g g fo g o C d C C C C t g p (b) (7)
Finally, the active zone was characterized by the cell diffusion coefficient (or motility) D o , the synthesis coefficient α g and the cellular decay factor d g .

Application to the experimental canine model and sensitivity analysis

The reference implant was a PMMA stable canine implant previously evaluated in-vivo [START_REF] Vestermark | Mechanical interface conditions affect morphology and cellular activity of sclerotic bone rims forming around[END_REF]) (Fig. 1b) . Boundary conditions and peri-implant tissue formation showed a polar symmetry with a variable level of mineralization in the radial direction r. The set of governing equations previously presented was translated into a custom made axisymmetric finite element formulation. Nine nodes quadratic elements were used with structural fraction φ s , growth factor concentration C g , osteoblast concentration C o and fluid pressure p as nodal variables. The associated computational model was implemented with the Python programming language (Python Software Foundation®) using Pysparse module to solve the spatio-temporal finite element matrix system.

The model is described in Figure 1c. The domain was radially defined by the implant radius (r i = 3.25 mm), the drill-hole (r d = 4.1 mm) and the surrounding trabecular bone (r s = 7 mm ≈ 2 × r i ). We assumed a 15 microns thick layer for the active zone considering implant surface roughness comprises between 0 and 5 microns [START_REF] Kim | Microrough titanium surface affects biologic response in MG63 osteoblast-like cells[END_REF][START_REF] Marinucci | Effect of titanium surface roughness on human osteoblast proliferation and gene expression in vitro[END_REF][START_REF] Meirelles | Increased bone formation to unstable nano rough titanium implants[END_REF] and mean osteoblast size of 10 microns [START_REF] Buser | Enhanced Bone Apposition to a Chemically Modified SLA Titanium Surface[END_REF]. The meshing was made of 201 nodes and 1206 DOF with steps evolving from 15 microns for the active zone, 60 microns for the intermediary zone and varying from 60 microns to 315 microns up to the surrounding trabecular bone.

As boundary conditions, fluid flux, cell flux and growth factor flux were nil at boundaries of the investigated domain shown in Figure 1c; the porosity φ f of the host trabecular bone was 50 %. As temporal initial conditions, osteoblast concentration of 1000 cell/mm 3 was homogenously distributed into the surrounding bone and post-operative gap, and growth factor concentration of 0.15 ng/mm 3 was homogenous into the gap. According to in-vivo experiment, the healing process was evaluated up to 8 weeks post-operatively. Model data are summarized in Table 1.

As shown previously, the implant active zone was characterized by the triplet (D o , d g , α g ).

This triplet in its wholeness was unavailable in the literature, especially including our specific constraints: osteoblastic lineage, TGF-β 1 and orthopedic implant in-vivo. Only data coming from cell assays could be used. First, decay factor d g was obtained from equation (7a) using a finite difference scheme in time imposing a decay of 5% each 100 minutes (Wakefield et al. 1990). This parameter was assumed intrinsic to the cell phase, and it was kept independent from the substrate because of a lack of available experimental data in the literature. Second, cell diffusion coefficient D o was extracted from (Dee et al.1999). Third, the synthesis coefficient α g was obtained from equations (7a) and (7b) using a finite difference resolution in time and TGF-β1 concentrations measured in presence of several types of substrates. PMMA, acid-etched titanium, coarse grit blasted acid-etched titanium, and coarse grit blasted acidetched titanium with RGDS peptide (Raush-fan et al., 2008). According to in-vitro experiments, the structural fraction was nil (φ s = 0) inducing a porosity φ f fixed to unity.

Since accessible by histomorphometry in case of PMMA implant, the distribution pattern of structural (or mineralized) fraction φ s was targeted to evaluate the reliability of the computational model. After this step, a numerical sensitivity analysis was conducted involving the surface-dependant doublet (D o , α g ) in case of PMMA and surface-modified titanium substrates.

Results

Computations were established using a decay factor d g of 8.5×10 -11 mm 3 /cell.s, obtained from the resolution of equation (7a). Because of lack of available data in literature, we assumed that d g was identical for the investigated bioactive surfaces, since intrinsic to the osteoblastic phase and established with the same growth factor concentration. However, parameters (D o , α g ) obtained from equation 7b were surface-dependant and summarized in Table 2.

The ex-vivo and computed radial distribution pattern of structural fraction φ s with PMMA implant were plotted in Figure 2. In digitized histological slices, the summation in concentric zones of structural fraction pixels was divided by the pixel overall summation of the zone of interest to derive the radial distribution of structural fraction. Mean value was established on a base of four slices. It can be observed in Figure 2 that magnitude and radial pattern were well predicted by the computational model.

Distribution pattern for the other three types of bioactive surfaces were also plotted in Figure 2. We observed significant differences since φ s at the implant surface, increased from 59 % for PMMA to 87% for the coarse grit blasted acid-etched with RGDS peptide. Differences were attenuated close from the drill hole and non-significant in the host bone whatever the implant constitutive material. We varied the active zone thickness (15 microns ± 20 %) and it confirmed the local effect of this parameter.

The numerical sensitivity analysis plotted in Figure 3 first confirmed that the decrease of diffusion coefficient D o globally increased the structural fraction φ s into the active zone. To go further, it appeared that the variation pattern of φ s to doublet (D o , α g ) could be described by a generic monotonic function increasing up to an optimum and then decreasing. This was clear for the acid-etched grit-blasted titanium with the optimum (point A) obtained for α g ≈ 16×10 - 11 ng/cell.s, and for the acid-etched grit-blasted titanium + RGDS which optimum was obtained for α g ≈ 38×10 -11 ng/cell.s (point B). For the acid-etched titanium, the optimum was not reached for the range of α g , and or the PMMA, the structural fraction was not influenced by α g . We evaluated a theoretical active surface with a particularly low diffusion coefficient D o of 1.25 1.25×10 -7 mm 2 /s, and we obtained a constant decrease of structural fraction; the optimum could potentially be reached for lower values of α g .

Discussion

Periprosthetic healing is an endochondral process involving complex biological events, and how biological medium reacts with the implant was underlying our approach. The mechanical aspect belonged to the domain of reactive transport in porous medium. We proposed an original set of diffusive-convective-reactive governing equations to investigate in time and space, the role of surface bioactive properties on cell flux, growth factor balance and finally structural (or mineral) fraction formation.

We proposed to model the bioactive properties of the implant using the triplet (D o , d g , α g ) locally updated in an active zone of few microns. Osteoblastic flux was governed by the diffusion coefficient D o that conditioned the random motility of cells due to adhesion sites and porosity of the surface. The haptotactic flux was related to local gradient of structural fraction φ s . The capability of implant surface to modify the balance of TGF-β1 was taken into account the cell decay factor d g and the synthesis coefficient α g . Although conservative equations ( 3) and ( 4) allowed it, effects on chemotaxis and cell proliferation were not directly updated by χ o and α o respectively, but the chemotactic flux was driven by the local growth factor gradient relied to α g , d g and local porosity φ f . This choice was due to a lack of available data in relevant literature. For the same reason, we were unable to update the local effect of bioactive surface on cell proliferation (α o in equation ( 3)), and structural (or mineral) matrix synthesis (α s in equation ( 6)).

We considered the mechanobiological interactions between osteoblasts and the implant by using the scale of cell population and growth factor phase. A multi-scale approach could potentially be relevant. This zone could be modeled at the cell scale with addition of local physico-chemical parameters. However, this approach should involve significant knowledge of surface properties and local cell-substrate interactions. This could be a perspective knowing that implementation of bio-mathematical governing equations at this scale should be challenging.

The relevance of our governing equations was initially evaluated in an in-vivo application, which required a custom-made finite element resolution. This application was a stable PMMA canine implant [START_REF] Vestermark | Mechanical interface conditions affect morphology and cellular activity of sclerotic bone rims forming around[END_REF]) whose distribution patterns of periprosthetic tissue were available by histomorphometry. About structural fraction φ s , the computational model showed good results in term of magnitudes and spatial distribution. This allowed establishing the basis of a numerical sensitivity analysis where the role of three type of bioactive surface was evaluated: acid-etched titanium, coarse grit blasted acid-etched titanium, and coarse grit blasted acid-etched titanium with RGDS peptide.

The numerical model could have shown numerical instabilities because of strong non-linear coupling. Previously, this risk has been evaluated by implementing a parametric sensitivity analysis (Ambard et al. 2004). Despite large variations of input parameters, up to ± 75 %, no instability was detected. The modeling of the active zone involving the cell motility, the coefficient of growth factor synthesis and the cellular decay, modified the diffusive terms and the reactive terms of governing equations. This could potentially provoke instabilities. After multiple calculations, we checked that no instability of the numerical scheme was induced by these significant perturbations.

The decrease of cell diffusion, acid-etched + RGDS peptide vs PMMA, and the increase of growth factor synthesis, PMMA vs acid-etched + RGDS peptide, significantly improved the amount of structural fraction φ s (or mineralized tissue) on the implant surface. Lower D o induced a longer stay of osteoblasts on the implant surface and higher α g increased concentration of TGF-β1. In consequence, osteoblast attraction by chemotaxis was enhanced and mineralized tissue formation was increased.

When the variation pattern of φ s to doublet (D o , α g ) was investigated, an envelope pattern involving an optimum seemed to be obtained (point A and B in Figure 3). However, when this optimum was exceeded especially for increased value of growth factor synthesis, we obtained a decrease of mineral fraction formation in the active zone. This could be explained by an untimely amount of cells and an early haptotactic migration towards the surrounding bone where the porosity gradient remained high. This result was corroborated by the literature regarding osseoinduction of surface modified substrate, and suggesting that osseointegration could not be optimal for strong surface modification and/or high dose of growth factors [START_REF] Svehla | Morphometric and mechanical evaluation of titanium implant integration: Comparison of five surface structures[END_REF][START_REF] Broderick | Depressed bone mineralization following high dose TGF-β1 application in an orthopedic implant model[END_REF]Buser et al., 2007;[START_REF] Colnot | Molecular analysis of healing at a bone-implant interface[END_REF][START_REF] Qu | The initial attachment and subsequent behavior regulation of osteoblasts by dental implant surface modification[END_REF].

We found that the distribution pattern of new-formed mineralized tissue was not significantly modified in the zones of post-operative gap and host trabecular bone. It confirmed that implant surface properties should play a limited role to reduce heterogeneity of new-formed tissue, which is unfavourable to the long-term survival of the arthroplasty. Our theoretical model involved a significant interest to help in the understanding of complex events in implant fixation. It involved a certain potential to enhance mixed technical-clinical strategy to reduce revision, avoiding empirical strategy. However, the robustness of a predictive model is strongly dependant upon the reliability of input data and because our model relied on interdisciplinary fields, its range of uncertainties was enlarged. Before targeting transfers in clinical setting, we aim at implementing a systems biology approach by which biological questions will be addressed through integrating in-vitro experiments with functionalized bioactive substrates in iterative cycles with computational modelling and theory.

Acknoledgement: J.E. Bechtold PhD and Prof. K. Søballe MD PhD.

Albrektsson, T., Brånemark, P.I., Hansson, H.A., Kasemo, B., Larsson, K., Lundström, I., McQueen, D.H., Skalak, R., 1983. The interface zone of inorganic implants in vivo: Titanium implants in bone. Annals of Biomedical Engineering 11, 1-27. Alliston, T.N., Derynck, R., 2000. Transforming growth factor-β in skeletal development and maintenance. In: Canalis, E. (Ed.), Skeletal Growth Factors. TPA, Philadelphia, 233-249. Ambard, D., Pedrono, A, Swider, P., 2004. A predictive numerical model of the periprosthetic tissue formation surrounding a stable implant, 50 th Annual Meeting of Orthopaedic Research Society, USA. Ambard, D., Bechtold J.E., [START_REF] Vestermark | Mechanical interface conditions affect morphology and cellular activity of sclerotic bone rims forming around[END_REF]. Prediction of the periprosthetic tissue mineralization in a stable implant coupling cell diffusion in a porous media, 6th Intl Symposium on Computer Methods in Biomechanics and Biomedical Engineering, Spain. Ambard, D., Bechtold J.E., [START_REF] Vestermark | Mechanical interface conditions affect morphology and cellular activity of sclerotic bone rims forming around[END_REF]

Captions for Tables and Illustrations

Table 1 -Parameters of the reference computational model (PMMA implant). Data were from relevant litterature: Wakefield et al.1990;[START_REF] Maheshwari | Deconstructing (and reconstructing) cell migration[END_REF][START_REF] Dee | Osteoblast population migration characteristics on substrates modified with immobilized adhesive peptides[END_REF]Tranqui et al. 2000;[START_REF] Vestermark | Mechanical interface conditions affect morphology and cellular activity of sclerotic bone rims forming around[END_REF]Bailon-Plaza et al 2001;Raush-fan et al. 2008. The experimental pattern concerned the PMMA implant and it was obtained from histomorphometry. The active zone of 15 microns thick was located at the implant surface and sensitivity analysis involved titanium substrate with three types of bioactive surface treatment: acid-etched, coarse-grit-blasted acid-etched, coarse-grit-blasted acid-etched with RGDS peptide.

Mechanobiology and implant surface properties

Guérin G., Ambard D., Swider P.

Parameters

  1

Figure 1 -

 1 Figure 1 -(a) The periprosthetic zone of interest comprises three sub-domains: the host trabecular bone, the post-operative gap and the active zone of few microns thick where implant bioactive properties modified growth factor sources and cell flux. (b) Canine experimental device: stable PMMA implant. (c) Axisymetric model of the canine implant device showing the radial location of zones of interest: host trabecular bone, post-operative gap, implant active zone.

Figure 2 -

 2 Figure 2 -Radial distribution patterns of structural fraction φ s for different types of implants.

  Geometry r i = 3.25 mm r s = 7 mm r d = 4.1 mm active zone : 15 μm Osteoblast C o = 1000 cell/mm 3 C om = 10000 cell/mm 3 α o = 1.9 ×10 -10 mm 6 /cell.ng.s α s = 6 × 10 -9 mm 6 /cell ng s h o = 1.8 × 10 -6 mm 5 /s.kg χ ο = 3 × 10 -5 mm 5 /s.ng D o = 2.5 × 10 -7 mm 2 /s Growth factor C g = 0.15 ng/mm 3 C gp =0.05 ng/mm 3 D g = 4.8 × 10 -6 mm 2 /s α g = 5 × 10 -12 ng/cell s d g = 8.5 × 10 -11 mm 3 /cell s

Figure 1

 1 Figure 1 -(a) The periprosthetic zone of interest comprises three sub-domains: the host trabecular bone, the post-operative gap and the active zone of few microns thick where implant bioactive properties modified growth factor sources and cell flux. (b) Canine experimental device: stable PMMA implant. (c) Axisymetric model of the canine implant device showing the radial location of zones of interest: host trabecular bone, post-operative gap, implant active zone.

Figure 3 -

 3 Figure3-Influence of bioactive surface treatment on the structural fraction φ s in the active zone located at the implant surface. Surface dependant parameters were the diffusion coefficient D o (10 -7 mm 2 /s) and the coefficient of growth factor production in autocrine mode α g (10 -12 ng/cell s). The cell decay factor was constant (d g = 8.5 ×10 -11 mm 3 /cell s). Materials were PMMA and titanium substrate with three types of bioactive surface treatment: acidetched, coarse-grit-blasted acid-etched, coarse-grit-blasted acid-etched with RGDS peptide.

  . Prediction of the periprosthetic tissue mineralization in a stable implant coupling cell diffusion in a porous media. 6 th Intl Symposium on Computer Methods in Biomechanics and Biomedical Engineering. experimental loaded implants. J Orthop Res. 22(3): 647-52. Viceconti, M., Muccini, R., Bernakiewicz, M., Baleani, M., Cristofolini, L.2000. Largesliding contact elements accurately predict levels of bone-implant micromotion relevant to osseointegration. J Biomech. 33(12):1611-8. Wakefield, L.M., Winokur, T.S., Hollands, R.S., Christopherson, K., Levinson, A.D., Sporn, M.B., 1990. Recombinant latent transforming growth factor β 1 has a longer plasma halflife in rats than active transforming growth factor β 1, and a different tissue distribution. Journal of Clinical Investigation 86, 1976-1984. Wang, J.Y., Wicklund, B.H., Gustilo, R.B., Tsukayama, D.T., 1997. Prosthetic metals interfere with the functions of human osteoblast cells in vitro. Clinical Orthopaedics and Related Research 339, 216-226. Wu, C.Y., Fields, A.J., Kapteijn, B.A.E., McDonald, J.A., 1995. The role of a4b1 integrin in cell motility and fibronectin matrix assembly. Journal of Cell Science 108, 821-829.

Table 2 -

 2 Bioactive properties of active zone according to four type of substrate. Values were computed from equation (7a) and (7b) using data from relevant literature:Wakefield et al.1990; Raush-fan et al. 2008. Coefficients for other tissue were α g = 5×10 -12 ng/cell s and D o = 2.5×10 -7 mm 2 /s.

Table 1 -

 1 Parameters of the reference computational model (PMMA implant). Data were from relevant litterature:Wakefield et al.1990;[START_REF] Maheshwari | Deconstructing (and reconstructing) cell migration[END_REF][START_REF] Dee | Osteoblast population migration characteristics on substrates modified with immobilized adhesive peptides[END_REF] Tranqui et al. 2000;[START_REF] Vestermark | Mechanical interface conditions affect morphology and cellular activity of sclerotic bone rims forming around[END_REF] Bailon-Plaza et al 2001; Raush-fan et al. 2008.

			active zone of implant
				coarse grit
	PMMA	acid-etched titanium	coarse grit blasted acid-etched titanium	blasted acid-etched titanium
				+ RGDS peptide
	d g (10 -11 mm 3 /cell s) 8.5 5 α g (10 -12 ng/cell s) D o (10 -7 mm 2 /s) 2	8.5 8.5 1.6	8.5 21 1.6	8.5 21 1.45

Table 2 -

 2 Bioactive properties of active zone according to four type of substrate. Values were computed from equation (7a) and (7b) using data from relevant literature:Wakefield et al.1990; Raush-fan et al. 2008. Coefficients for other tissue were α g = 5×10 -12 ng/cell s and D o = 2.5×10 -7 mm 2 /s.Mechanobiology and implant surface propertiesGuérin G.,Ambard D., Swider P.

Mechanobiology and implant surface properties Guérin G., Ambard D., Swider P.