Emmanuel Boissard

Serge Cohen

Thibault Espinasse

James Norris
email: j.r.norris@statslab.cam.ac.uk

Diffusivity of a random walk on random walks

Keywords: Random walk, Graph, Central limit theorem AMS classification (2000): 05C81, 60F05

We consider a random walk

with the constraint that each coordinate of the walk is at distance one from the following one. In this paper, we show that this random walk is slowed down by a variance factor σ 2 K = 2 K+2 with respect to the case of the classical simple random walk without constraint.

1 Presentation of the random walk

Let Z (1) n , • • • , Z (K+1) n
∈ Z K+1 denote the heights of K + 1 simple random walks on Z, conditioned on satisfying

∀n ∈ N, ∀i ∈ [1, K] , Z (i+1) n -Z (i) n = 1. (1)
More precisely, the random walk is a Markov chain on the state space of K-step walks in Z S K = {(z (1) , . . . , z

(K+1) ) ∈ Z K+1 , ∀i ∈ [1, K] , |z (i+1) -z (i) | = 1}
where the next step from z ∈ S K is selected uniformly among the neighbours of z in the usual lattice Z K+1 that belong to S K . In other words, we consider K +1 simple random walks on the lattice Z coupled under a shape condition.

As in the case of a simple random walk, the rescaled trajectory of a walker, say Z (1) , will converge in law to a Brownian motion. However, it is interesting to note that the constraint between each coordinate only slow down the walk by decreasing its variance.

Since it is classical to illustrate for our students the simple random walk as the motion of a drunk man, we can illustrate the previous mathematical fact by considering the random walk as the motion of a chain of prisoners. It should convince even non mathematicians that the motion of the walk is slowed by the constraint. However it seems very hard to guess the variance from this comparison ! More precisely, denote ∀t ≥ 0, ∀n ∈ N, ξ

(n) t = Z (1) ⌊nt⌋ √ n ,
where ⌊x⌋ is the integer part of the real number x.

THEOREM 1. The rescaled random walk ξ

(n) t t≥0
converges in law, as n goes to infinity, to a Brownian motion with variance

σ 2 K = 2 K + 2 .
Convergence to the Brownian motion is the usual invariance principle : the noteworthy statement here is that it is possible to give an explicit expression for the limit diffusivity of the process, and that its expression is particularly simple.

Our object of interest, the motion of Z (1) , is a non-Markovian process that falls into the class of random walks with internal structure. Related questions of limit diffusivity for random walks conditioned to respect some geometric shape have been studied in the literature, under the name of "spider random walks" or "molecular spiders", see [START_REF] Antal | Molecular Spiders in One Dimension[END_REF]. The computation of the limit diffusion coefficient is also a central aim there, although the model and methods are different.

Our initial motivation was however more remote. Actually we first addressed this question starting from combinatorial problems related to 6-Vertex model in relation with the Razumov Stroganov conjecture (See [START_REF] Cantini | Proof of the Razumov-Stroganov conjecture[END_REF] for instance). The problem can also be related to random graph-homomorphisms (See [START_REF] Benjamini | Random graphhomomorphisms and logarithmic degree[END_REF]) or the Square Ice Model (See [START_REF] Lieb | Residual entropy of square ice[END_REF] where Z (1) evolves on a torus.)

Roughly speaking we can say that, in the literature we read, the authors consider questions related to the uniform distribution on a sequence of finite graphs G K and wonder about various asymptotics when K → +∞. Later in the article the evolution of Z n will be described as the simple random walk on a graph G K . Hence on the one hand our problem is a very simplified version of the problems stated above, On the other hand we were surprised to have such a simple formula for σ 2 K = 2 K+2 , which is true for all K and not only for K → +∞. We thought in the beginning that the proof of this fact should be simple but it turns out that, although elementary, the tools used to obtain the result are more sophisticated than expected. It is the aim of this note to show these tools.

To prove the theorem, we will look for a decomposition

Z (1) n = M n + f (Z n-1 , Z n ),
where (M n ) n∈N is a martingale, and where f : Z K+1 × Z K+1 → R is a bounded function. We will then show that the following limit exists :

σ 2 K = lim n→∞ 1 n E M 2 n .
and that it is indeed the desired diffusivity. The path to this conclusion is akin to classical results for Central Limit Theorems for Markov chains (E.g. [START_REF] Meyn | Markov chains and stochastic stability[END_REF]).

We will use another equivalent, albeit more geometric, point of view on this decomposition. We split the chain in two parts : on the one hand, the motion of one of the walkers, and on the other hand, the relative positions of the walkers (which we call the "shape" of the chain at a given time). The latter part is a Markov chain over the state space {-1, 1} K and our quantity of interest is (almost) an observable of this chain. Computing the martingale decomposition that we wish for amounts to decomposing a discrete vector field over this new state space into a divergence-free part (corresponding to the martingale part) and a gradient part (corresponding to the function f ). Owing to a particular geometric property of this vector field, for which we coin the term "stationarity", it is indeed possible to perform this calculation explicitly.

The random walk with constraint

Let us denote

∀n ∈ N, ∀i ∈ [1, K] , Y (i) n = Z (i+1) n -Z (i) n ,
and

Y n = Y (1) n , • • • , Y (K) n ∈ {-1, 1} K .
Here Y n describes the shape of Z

(1)

n , • • • , Z (K+1) n ∈ Z K+1 , i.e. the position of each Z (i)
n relatively to the previous one, and belongs to V K = {-1, 1} K , whereas Z

(1) n can be seen as the height of the first walker. Obviously, the evolution of the chain of walkers may be described by the variables

Z (1) n , Y n n∈N .
For a convenient analysis, we will represent our process as the simple random walk on a (multi)-graph G K , which we define below.

Set

V K = {-1, 1} K : the multi-graph G K is given as a triplet (V K , E + K , E - K ) where E + K , E - K ⊂ V K × V K
are two edge sets, called respectively the set of "positive" and "negative" edges. A couple (a, b) ∈ V K × V K , a = b, belongs to E + K if the vector ba ∈ {-2, 0, 2} K has nonzero entries of alternating signs, with the first one negative. Moreover, E + K also contains a loop from each a ∈ V K to itself, noted (a, a) + .

Likewise, E - K contains those couples (a, b) ∈ V K × V K , a = b, such that ba ∈ {-2, 0, 2} K has nonzero entries of alternating sign, with the first one positive, and self-loops noted (a, a) -for each a ∈ V K .

Set

E K = E + K ∪ E - K .
Finally, we consider the following function on E K .

DEFINITION 1. Let A : E K → {1, -1} be the function that takes the value 1 on E + K and -1 on E - K . Note that we have, for any a ∈ V K , A((a, a) ± ) = ±1.

PROPOSITION 1. Let (W K n ) n≥0 be the simple random walk on G K . The processes Y n , Z (1) n n∈N and W K n , n j=1 A(W K j-1 , W K j ) n∈N have the same distribution.
Proof. It is sufficient to prove that the two Markov chains (taking values in

{-1, 1} K × Z) Y n , Z (1)
n n∈N and W K n , n j=1 A(W K j-1 , W K j )
n∈N have the same transition matrix.

We claim that

L Y n+1 -Y n , Z (1)
n+1 -Z (1) n |Y n , Z (1)
n = L   W n+1 -W n , A(W K n , W K n+1 )|W n , n j=1 A(W K j-1 , W K j )   . Recall that S K = z ∈ Z K+1 , ∀i ∈ [1, K], z (i+1) -z (i) = 1 . Denote A = (z 1 , z 2 ) ∈ S 2 K , ∀i ∈ [1, K + 1], z (i) 1 -z (i) 2 = 1 .
Since (Z n ) n∈N and (W n ) n∈N are both simple random walks, the corresponding transition matrices are given respectively by

P Z z 1 ,z 2 = 1 Card {z ∈ S, (z 1 , z) ∈ A} 1 1 (z 1 ,z 2 )∈A P W w 1 ,w 2 = 1 Card w ∈ {-1, 1} K , (w 1 , w) ∈ E K 1 1 (w 1 ,w 2 )∈E K .
Denote now

δ : S → {-1, 1} K (z (1) , . . . , z (K+1) ) → (z (2) -z (1) , . . . , z (K+1) -z (K) ), so that Y n = δ(Z n ), n ∈ N.
To prove the proposition, it is sufficient to prove that

(z 1 , z 2 ) ∈ A ⇔ z (1) 2 -z (1) 1 = 1 and δ(z 2 ) -δ(z 1 ) ∈ E ǫ K ,
where ǫ = -sgn(z

1 ) (and sgn is the sign function).

First note that, if ǫ = ±1, and

if ∀i ∈ [1, K + 1], z (i) 1 = z (i) 2 -ǫ, then (z 1 , z 2 ) ∈ A. Moreover, we have in this case δ(z 1 ) = δ(z 2 ), and (δ(z 1 ), δ(z 2 )) ∈ E ǫ K as claimed. On the other hand, if δ(z 1 ) = δ(z 2 ) and z (1) 2 -z (1) 1 = ǫ ∈ {-1, 1}, then ∀i ∈ [1, K + 1], z (i) 1 = z (i) 2 + ǫ. For the chain, it means that, if Y n+1 = Y n then Z (1) n+1 -Z (1) n is +1, -1 with equal probability 1 2 independently of Y n , Y n+1 , Z (1)
n . Assume now that z 1 = z 2 . We will prove that (δ(z 1 ), δ(z 1 ))

∈ E ǫ K , where ǫ = z (1) 1 -z (1)
2 , i.e. δ(z 2 )δ(z 1 ) ∈ {-2, 0, 2} K , with nonzero entries of alternating signs, and a first one of the sign of ǫ.

Indeed, there is a first index

τ 1 = inf{i ∈ [1, K] , such that δ(z 1 ) (i) = δ(z 2 ) (i) }.
There are two possible cases

δ(z 1 ) (τ 1 ) = +1, δ(z 2 ) (τ 1 ) = -1, z (τ 1 ) 2 = z (τ 1 ) 1 + 1, z (τ 1 +1) 2 = z (τ 1 +1) 1 -1, δ(z 1 ) (τ 1 ) = -1, δ(z 2 ) (τ 1 ) = +1, z (τ 1 ) 2 = z (τ 1 ) 1 -1, z (τ 1 +1) 2 = z (τ 1 +1) 1 + 1.
In the two cases we have

A(δ(z 1 ), δ(z 2 )) = z (τ 1 ) 2 -z (τ 1 ) 1 = z (1) 2 -z (1)
1 . Furthermore, we have

δ(z 2 ) (τ 1 ) -δ(z 1 ) (τ 1 ) = -2A(δ(z 1 ), δ(z 2 )). Then if τ 1 < K let us define τ 2 = inf{i ∈ [τ 1 + 1, K] , such that δ(z 2 ) (i) = δ(z 1 ) (i) },
where we set by convention τ 2 = K +1, if the condition defining the infimum is never satisfied.

Using the same arguments, we get

δ(z 2 ) (τ 2 ) -δ(z 1 ) (τ 2 ) = 2A(δ(z 1 ), δ(z 2 )).
By induction one can define

τ j = inf{i ∈ [τ j-1 + 1, K] , such that δ(z 1 ) (i) = δ(z 2 ) (i) }, until τ j-1 ≥ K. We have δ(z 2 ) (τ j ) -δ(z 1 ) (τ j ) = 2(-1) j A(δ(z 1 ), δ(z 2 )).
By definition, we get (δ

(z 1 ), δ(z 2 )) ∈ E ǫ K , where ǫ = z (1) 1 -z (1) 2 . On the other hand, if (δ(z 1 ), δ(z 2 )) ∈ E ǫ K , where ǫ = z (1) 1 -z (1) 2 ∈ {-1, 1}, then one can recover explicitly (z 1 , z 2 ) from the definition of δ. Moreover, the condition ∀i ∈ [1, K + 1], z (i) 2 -z (i) 1
= 1 is implied by the previous arguments (following the definitions of the τ j ).

We may denote, Y n e → Y n+1 when e ∈ E + (Z

(1) n+1 = Z (1) n + 1) and Y n e ← Y n+1 when e ∈ E -(Z (1) n+1 = Z (1)
n -1). For a general a ∈ V K the previous enumeration of its neighbors is surprisingly complicated but we can provide some simple examples.

For instance if a = (1, . . . , 1) ∈ V K , has only K + 2 neighbors:

a ← a a → a ∀i ∈ [1, K], a → (1, • • • , -1, • • • , 1), ( 2
)
where the -1 is in the i-th position.

Note now that the graph G K can also be described inductively : there are only six following possibilities for (Z

(1) n , Z (2) n , Z (1) n+1 , Z (2)
n+1 ), described below: In the figure, we have used the concatenation notation : given a string a, the string 1a, resp. (-1)a, is obtained by adding a 1, resp. -1 in front of a. Looking only at the 3 cases such that Z

(1)

n+1 -Z (1)
n = 1, we can deduce the construction of G + K+1 from G + K :

V K+1 = {-1, 1} K+1 E + K+1 = (1a, 1b), (a, b) ∈ E + K ∪ ((-1)a, (-1)b), (a, b) ∈ E + K ∪ (1a, (-1)b), (a, b) ∈ E - K Figure 1
shows the first two graphs G + 1 and G + 2 . Note that each edge of G 1 gives 3 edges for G 2 , one on each facet {1a, a ∈ V K }, and {(-1)a, a ∈ V K } and one crossing from the facet {1a, a ∈ V K } to the facet {(-1)a, a ∈ V K }.

We obtain the cardinality D K of E K (as a multigraph) by induction: We will also make use of the number δ K of edges of the form (1a, (-1)b) ∈ E + K which can be computed by induction:

D K = 2.3 K .
δ K = 3 K-1 .
Let us now describe vector fields on this graph.

Vector fields on graphs

In the previous section a function A has been defined on edges of G K . We will consider here A as a vector field on G K .

Definitions

Let us first recall some classical definitions.

DEFINITION 2 (Vector fields). A vector field on

G K = (V K , E K ) is a func- tion S : E K → R such that ∀(a, b) ∈ E K , a = b, S(a, b) = -S(b, a). (3)
and such that, for any a ∈ V K , S((a, a) + ) = -S((a, a) -).

DEFINITION 3 (Gradient vector fields). We say that the vector field S on G K is a gradient vector field if there exists a function f on the vertices of

G K such that for each edge (a, b) ∈ E K , S(a, b) = f (b) -f (a).
The gradient vector field associated with f is denoted by ∇f .

DEFINITION 4 (Divergence and divergence-free vector fields). The divergence of a vector field S at point x ∈ G K is defined by

(∇ • S)(x) = (x,y)∈E K S(x, y).
We say that a vector field S is divergence-free if its divergence vanishes at all points.

We can endow the set of vector fields with a scalar product

S, S ′ E K = 1 D K e∈E K S(e)S ′ (e).
Please note that the sum runs over all edges, including loops (a, a) ± . Denote also, for any subsets φ, ψ of V K ,

J φ,ψ (S) = a∈φ,b∈ψ,(a,b)∈E K S(a, b),
the flux of S going from φ to ψ. Note that a divergence free field B verifies ∀φ ⊂ V K , J φ,V K \φ (B) = 0.

Hodge decomposition of vector fields on graphs

In analogy with the case of vector fields in Euclidean spaces, we can decompose any vector field on G K into the sum of a gradient vector field and a divergence-free field. The following proposition is well-known.

PROPOSITION 2. Let S be a vector field on G K . There exist a unique gradient vector field ∇f and a unique divergence-free field B such that

S = ∇f + B. Moreover ∇f, B E K = 0
The last identity simply means that gradient fields and divergence-free fields are orthogonal complements of each other in the vector space of vector fields over G K .

In our case we are interested in stationary vector fields.

DEFINITION 5 (Stationary vector field). A subgraph G of the complete graph on {-1, 1} K is stationary, if the following holds. For u, u

′ ∈ {-1, 1} K and v ∈ R k such that u + v, u ′ + v ∈ {-1, 1} K , if (u, u + v) is an edge of G, then (u ′ , u ′ + v) is an edge of G.
A vector field S defined on a subgraph of the complete graph on {-1, 1} K with a stationary domain is stationary if for all (u, v) edges of G, S(u, v) only depends on uv (where {-1, 1} K is embedded in R K in an obvious way).

REMARKS 1. Note that, thanks to the construction of

E + K , it is stationary. Moreover if (u, u + v) ∈ E + K and u ′ , u ′ + v ∈ V K , then (u ′ , u ′ + v) ∈ E + K
and thus the vector field A taking values 1, resp. -1, on E + K , resp. E - K , is stationary. So we may expect the gradient vector field ∇f in the Hodge decomposition of A to be stationary. Unfortunately if S is a stationary vector field on G K and its decomposition is

S = ∇f + B
as per Proposition 2, then ∇f is not always stationary.

Nevertheless it turns out that the gradient vector field ∇f in the Hodge decomposition of A is actually stationary as it will be shown in the next section.

Hodge decomposition of A

Let us recall the Definition 1 the vector field A on G K is such that

∀e ∈ E + K , A(e) = 1 ∀e ∈ E - K , A(e) = -1.
In this section our aim is to compute a function f such that

(∇ • A)(a) = (∇ • (∇f ))(a) ∀a ∈ V K . ( 4
)
One can first remark that ∀a ∈ V K

(∇ • A)(a) = Card{b ∈ V K , (a, b) ∈ E + K } -Card{b ∈ V K , (a, b) ∈ E - K }.
(5) In the previous equation we used the notation Card for cardinality of sets. We will introduce various notations related to other cardinalities

α(a) = Card{b ∈ V K , (a, b) ∈ E + K } (6) ᾱ(a) = Card{b ∈ V K , (a, b) ∈ E - K }. (7)
Then we need also to define for 0 Let us consider the function

≤ k ≤ K, α k (a)
f 1 on V K such that f 1 (a) = α 1 (a) -ᾱ1 (a) = K i=1 a i . ( 8
)
The last equation is trivial since α 1 (a) is the number of the a i s equal to 1 and ᾱ1 (a) the number of a i s equal to -1. Obviously we also get

α 1 (a) + ᾱ1 (a) = K (9)
for any vertex a in V K . Let us remark that for any function

f on V K (∇ • (∇f ))(a) = (a,b)∈E K f (b) -f (a).
Since f 1 yields the sum of the digits of any vertex, we first observe that if (a, b) ∈ E K and the number of digits i ∈ {1, . . . , K} such that a i = b i is even then f 1 (b)f 1 (a) = 0. If the number of digits i ∈ {1, . . . , K} such that

a i = b i is odd and (a, b) ∈ E + K then f 1 (b) -f 1 (a) = -2. One can then deduce that (∇ • (∇f 1 ))(a) = -2(α od (a) -ᾱod (a)). ( 10
)
It turns out that if we consider the function f 2 on V K such that

f 2 (a) = α 2 (a) -ᾱ2 (a), (11) then
(∇ • (∇f 2 ))(a) = -(K + 2)(α ev (a) -ᾱev (a)). ( 12
)
We will prove (12) by induction on K. To do that we split (∇ • (∇f 2 )) into the sum of two functions

φ(a) = (a,b)∈E + K f 2 (b) -f 2 (a) (13) φ(a) = (a,b)∈E - K f 2 (b) -f 2 (a). ( 14
)
To proceed the induction argument we remark that for any vertex a in

V K f 1 (1a) = f 1 (a) + 1 (15) f 1 (-1a) = f 1 (a) -1 (16) f 2 (1a) = f 2 (a) + ᾱ1 (a) (17) f 2 (-1a) = f 2 (a) -α 1 (a). (18)
We will then compute φ(1a), φ(-1a), φ(1a), φ(-1a). The easiest computation is

φ(1a) = (1a,b)∈E - K+1 f 2 (b) -f 2 (1a) = (a,c)∈E - K f 2 (1c) -f 2 (1a), since (1a, -1c) ∈ E + K+1 . Then φ(1a) = (a,c)∈E - K f 2 (c) -f 2 (a) + (a,c)∈E - K ᾱ1 (c) -ᾱ1 (a).
To evaluate

One also get in the same way

φ(-1a) = φ(a) -α od (a). ( 20
)
The induction is a bit more involved for

φ(1a) = (1a,b)∈E + K+1 f 2 (b) -f 2 (1a) = (1a,1c)∈E + K+1 f 2 (1c) -f 2 (1a) + (1a,-1c)∈E + K+1 f 2 (-1c) -f 2 (1a).
Because of (17)

(1a,1c)∈E + K+1 f 2 (1c) -f 2 (1a) = (a,c)∈E + K f 2 (c) -f 2 (a) + (a,c)∈E + K ᾱ1 (c) -ᾱ1 (a)
= φ(a) + α od (a).

Then

(1a,-1c)∈E + K+1 f 2 (-1c) -f 2 (1a) = (a,c)∈E - K f 2 (c) -f 2 (a) - (a,c)∈E - K α 1 (c) + ᾱ1 (a) = φ(a) - (a,c)∈E - K α 1 (c) + ᾱ1 (c) + (a,c)∈E + K ᾱ1 (c) -ᾱ1 (a) = φ(a) -K ᾱ(a) -ᾱod (a). Hence φ(1a) = φ(a) + φ(a) + α od (a) -ᾱod (a) -K ᾱ(a). (21)
Similarly we get

φ(-1a) = φ(a) + φ(a) -α od (a) + ᾱod (a) + K ᾱ(a). (22)
We can now evaluate the functions φ, φ. Proof. We will only sketch the proof performed by an easy induction on K for φ, computations are similar for φ. Let us assume that (23), (24) hold for K, we have to compute φ(1a) and φ(-1a) and check that they fulfill (23), (24) for K + 1. Because of ( 21)

LEMMA 1. ∀a ∈ V K , φ (
φ(1a) = φ(a) + φ(a) + α od (a) -ᾱod (a) -K ᾱev (a) -K ᾱod (a).
One can check that ᾱev (1a) = ᾱev (a), α ev (1a) = α ev (a) + ᾱod (a), α od (1a) = α od (a) + ᾱev (a) + 1, ᾱod (1a) = ᾱod (a). Using (23) for K, we get (23) for K + 1 and φ(1a). The computations for φ(-1a) are left to the reader.

By summing ( 23) and (24) we get (12), and we deduce that if we take

f = - 1 2 f 1 + 1 K + 2 f 2 + K 2 , (25)
A -∇f is divergence free. Please note that obviously the additive constant in ( 25) is arbitrary but it yields the following convenient expression of f (a) in terms of the digits of a ∈ V K .

LEMMA 2.

f (a) = i, a i =-1 F K i , (26)
where

F K i = 3 + 2(K -i) K + 2 . ( 27
)
Moreover ∇f is a stationary gradient vector field.

Proof. Since we already know that f 1 (a) = K i=1 a i , it is enough to show by induction on K that

f 2 (a) = 1 2 K i=1 a i (1 + K -2i). ( 28
)
Obviously the formula is true for K = 1. Because of ( 17) and (18),

f 2 (1a) + f 2 (-1a) = 2f 2 (a) -f 1 (a) f 2 (1a) -f 2 (-1a) = K.
Hence (28) is proved for K + 1. Owing to f 1 (a) = K i=1 a i , and (28), it is obvious that ∇f 1 , ∇f are stationary and consequently ∇f is a stationary gradient vector field.

Proof of theorem 1

We denote by (∇f, B) the decomposition of A as per Proposition 2.

Back to the original problem, we recall that ∀n ≥ 0, Z

n+1 -Z (1) n = B(Y n , Y n+1 ) + ∇f (Y n , Y n+1 ). (1)
Let us denote

M n = Z (1) n -f (Y n ). Let F n = σ((Y k , Z (1)
k ), k ≤ n), we have

E[M n+1 -M n |F n ] = (∇ • B)(Y n ) = 0,
and (M n ) n∈N is a martingale.

We now sketch out how to apply the Central Limit Theorem for Markov chains. Let e n = (Y n , Y n+1 ), (e n ), n ≥ 0 is a Markov chain on E K . Then our quantity of interest Z [START_REF] Benjamini | Random graphhomomorphisms and logarithmic degree[END_REF] n is an additive observable of the process (e n ), n ≥ 0, as Z (1) n -Z

(1)

0 = n-1 k=0 Z (1) k+1 -Z (1)
k = n-1 k=0
A(e k ).

The Central Limit Theorem for Markov chains (see e.g. [START_REF] Meyn | Markov chains and stochastic stability[END_REF]) shows that ( 1

√ n Z (1)
⌊nt⌋ ) t≥0 converges as n → +∞ to a Brownian motion, with variance given by

σ 2 K = lim n→+∞ 1 n E Z (1) n 2 = lim n→+∞ 1 n E M 2 n + lim n→+∞ 1 n E f (Y n ) 2 + lim n→+∞ 2 n E [(f (Y n )M n )] .
Let us first compute lim n→+∞

1 n E M 2 n . We can remark that M 2 n - n-1 i=0 B(Y i , Y i+1 ) 2 is a F n martingale, hence lim n→+∞ 1 n E M 2 n = lim n→+∞ 1 n E n-1 i=0 B(Y i , Y i+1 )
If µ denotes the invariant measure for the random walk (Y n ) n≥0 , by ergodicity, we get

lim n→+∞ 1 n E n-1 i=0 B(Y i , Y i+1 ) 2 = E µ (B(Y 0 , Y 1 )) 2 .
Under µ the distribution of (Y 0 , Y 1 ) is uniform on E K because Z 1 is uniformly chosen among all neighbours of Z 0 , hence

lim n→+∞ 1 n E(M 2 n ) = B 2 , (29)
Using that f does not depend on n, we obtain

lim n→+∞ 1 n E f (Y n ) 2 = 0,
and by Cauchy Schwarz inequality and (29)

lim n→+∞ 2 n E [(f (Y n )M n )] = 0.
So we get

σ 2 K = lim n→+∞ 1 n E n-1 i=0 B(Y i , Y i+1 ) 2 ,
and by ergodicity,

σ 2 K = E µ (B(Y 0 , Y 1 )) 2 .
Since, under µ, the distribution of (Y 0 , Y 1 ) is uniform on

E K σ 2 K = B 2 .
By orthogonality of B and ∇f , σ 2 K = A 2 -A, ∇f . Thus, it remains to compute A, ∇f (since A = 1, by definition). At this point we use the fact that ∇f is a stationary field, in the sense of Definition 5. Then if we denote by a = (1, . . . , 1), z i = (1, . . . , -1, 1, . . . , 1) (where -1 is in i-th position) and, for 1 ≤ i ≤ K, we have by ( 26)

F i = ∇f (a, z i ), Now we compute A, ∇f as a function of F 1 , the value ∇f on the edge (1, • • • , 1) → (-1, 1, • • • , 1). By (3) A, ∇f = 2 D K (a,b)∈E + K (∇f )(a, b). If φ = {1a ′ , a ′ ∈ V K-1 } and ψ = {(-1)b ′ , b ′ ∈ V K-1 }, then using the relation between E K-1 and E K , A, ∇f = 2 D K (a,b)∈E + K ∩φ (∇f )(a, b)+ 2 D K (a,b)∈E + K ∩ψ (∇f )(a, b)+ 2 D K J φ,ψ (∇f ).
By stationarity of ∇f,

A, ∇f = 4 D K (a,b)∈E + K ∩φ (∇f )(a, b) + 2 D K J φ,ψ (∇f ) = 4 D K (a,b)∈E + K ∩φ (∇f )(a, b) + 4 D K J φ,ψ (∇f ) - 2 D K J φ,ψ (∇f ).
Then, because of the definition of φ and ψ,

A, ∇f = 4 D K (a ′ ,b ′ )∈E + K-1 (∇f )(1a ′ , 1b ′ ) + (∇f )(1b ′ , (-1)a ′ ) - 2 D K J φ,V K \φ (∇f ), = 4 D K (a ′ , b ′ )∈E + K-1 f (1b ′ ) -f (1a ′ ) + f ((-1)a ′ ) -f (1b ′ ) - 2 D K J φ,V K \φ (A).
Then by stationarity of ∇f

A, ∇f = 4 CardE + K-1 D K F 1 -2 δ K D K = 2 K + 2 .
REMARKS 2. In the proof, the way we guessed ( 27) is a bit mysterious.

Assuming that ∇f is a stationary gradient vector field, the family (F i ) i∈ [1,K] can be computed considering the system of equations given by ∀i ∈ [1, K], J M i ,N i = 0, where M i = (a j ) j∈[1,K] ∈ V K , a i = 1 and N i = (a j ) j∈[1,K] ∈ V K , a i = -1 . This leads to the following system:

          3 K-1 -3 K-2 . . . -3 -1 -3 K-2
. . . 3 K-1 -3 K-2 -1 -3 . . . -3 K-2 3 K-1           F =        3 K-1 3 K-2 . . . 3 1       
.

The unique solution is given by : ∀j ≤ K, F j = 3 + 2(Kj) K + 2 .

Even if the guess is correct, we did not find another way as the techniques used in the Section 3.3 to show that ∇f is a stationary gradient vector field.

Some further questions

In this final section we briefly outline some related problems.

• It is possible to make sense of the process when K is infinite. Several questions arise : what happens to the process of one marked walker ? Is there a scaling limit under equilibrium for the "shape"process (Z (k) -Z (1) ) 1≤k ?

Another natural step would be to let K grow with N in a suitable way, so as to get a scaling limit for the two-parameter process (Z

(k) n - Z (1)
n ) 1≤k≤K+1,n∈Z .

• One may also ask about different quantities, such as the diameter of the set of walkers under the invariant measure for the entire walk.

• One may also consider random walkers conditioned on satisfying different shape constraints, and on graphs more general than Z. As a starting example, what happens if we work on a torus, i.e. if we force also |Z

(K+1) n -Z (1)
n | = 1 ? The "shape" chain changes in this case and it is no longer irreducible over {-1, 1} K+1 (one may check that the number of -1 symbols is fixed, and that this enumerates the recurrence classes). It is interesting to point out that this setup is the one chosen by E. Lieb for the computation of the "six-vertex constant" in [START_REF] Lieb | Residual entropy of square ice[END_REF].

Acknowledgments

We warmly thank Charles Bordenave for fruitful conversations and useful comments.

Figure 1 :

 1 Figure 1: Construction of G + 1 and G + 2 .

  the number of the vertices b such that (a, b) ∈ E + K with k digits i ∈ {1, . . . , K} such that a i = b i . Similarly ᾱk (a) is the number of the vertices b such that (a, b) ∈ E - K and k digits i ∈ {1, . . . , K} such that a i = b i . Then we consider α ev (a) = k even α k (a) and α od (a) = k odd α k (a). Similarly ᾱev (a) = k even ᾱk (a) and ᾱod (a) = k odd ᾱk (a).

  (a,c)∈E - K ᾱ1 (c) -ᾱ1 (a) we use that ᾱ1 (a) is the number of digits equal to -1 in a. If (a, c) ∈ E - K and if the number of a i = c i is even then ᾱ1 (c) = ᾱ1 (a). If this number is odd then ᾱ1 (c) = ᾱ1 (a) -1. Therefore (a,c)∈E - K ᾱ1 (c) -ᾱ1 (a) = -ᾱod (a). Hence φ(1a) = φ(a) -ᾱod (a).

  a) = ᾱev (a) -(K + 1)α ev (a) + α od (a) + ᾱod (a) (23) φ(a) = -α ev (a) + (K + 1) ᾱev (a) -ᾱod (a)α od (a) (24)

  3 K-1 . . .

	. . .	. . .	. . .	-9 . . .	-3 . . .
	-3	-9


Geneviève Almouzni

Aline V Probst
email: aline.probst@univ-bpclermont.fr

Heterochromatin maintenance and establishment: lessons from the mouse pericentromere

Defined as a chromatin structure that remains condensed throughout the cell cycle heterochromatin is generally transcriptionally silent and is characterized by a specific molecular signature. Constitutive heterochromatin at the pericentromere is a conserved feature throughout evolution, which impacts genome stability. Here, we will summarize recent advances in our understanding of the dynamics of mouse pericentric heterochromatin during the cell cycle and development. Comparison with heterochromatin maintenance in fission yeast will enable discussions of the common basic principles and various mechanisms exploited in the distinct organisms.

Introduction

Initially discovered based on differential staining with DNA dyes, heterochromatin is defined as chromosomal domains that, in contrast to euchromatin, remain condensed throughout the cell cycle [START_REF] Heitz | Das Heterochromatin der Moose[END_REF] . Studies in different model organisms have identified heterochromatin regions close to centromeres and telomeres, which make up constitutive heterochromatin 2 . These forms of heterochromatin are critical for genome maintenance and are considered to be vital for repression of transposable elements. Especially, heterochromatin that forms at repetitive elements at the pericentromeres has been characterized extensively and emerged as the paradigm for constitutive heterochromatin. Here, we will summarize our knowledge on pericentric heterochromatin establishment and dynamics during murine development. We will also compare the features of pericentric heterochromatin in fission yeast and mice to illustrate possible common principles of heterochromatin maintenance.

While formation of a repressive chromatin environment at pericentric repeats is evolutionarily conserved, in different organisms distinct combinations of chromatin marks characterize pericentric heterochromatin at the molecular level. These chromatin marks include DNA methylation, specific histone variants, hypoacetylation and repressive histone modifications, in particular H3K9 methylation, which constitutes a binding site for the Heterochromatin protein 1 (HP1) family. A shared element, whose exact importance remains to be elucidated, is the clustering of pericentric heterochromatin into chromocenters in interphase nuclei (Figure 1). This spatial arrangement, which can compartmentalize heterochromatin away from the rest of the genome, might serve to create a local environment where silencing factors could possibly concentrate and help maintain the heterochromatin status. In fission yeast, the three chromosomes bring their centromeres together in a single chromocenter [START_REF] Funabiki | Cell cycle-dependent specific positioning and clustering of centromeres and telomeres in fission yeast[END_REF] while in mouse as well as in certain plants such as Arabidopsis, the arrangement of the chromosomes into subgroups in interphase gives rise to several chromocenters [START_REF] Guenatri | Mouse centric and pericentric satellite repeats form distinct functional heterochromatin[END_REF][START_REF] Fransz | Interphase chromosomes in Arabidopsis are organized as well defined chromocenters from which euchromatin loops emanate[END_REF] .

Interestingly, chromocenter organization is globally maintained even through S-phase [START_REF] Quivy | A CAF-1 dependent pool of HP1 during heterochromatin duplication[END_REF] and dissociation in individual chromosomes is only observed during mitosis. However, this organization can differ in distinct cell types and developmental stages as has been documented in mouse [START_REF] Mayer | Common themes and cell type specific variations of higher order chromatin arrangements in the mouse[END_REF][START_REF] Terranova | The reorganisation of constitutive heterochromatin in differentiating muscle requires HDAC activity[END_REF][START_REF] Solovei | Nuclear architecture of rod photoreceptor cells adapts to vision in mammalian evolution[END_REF][START_REF] Haaf | Higher order nuclear structure in mammalian sperm revealed by in situ hybridization and extended chromatin fibers[END_REF][START_REF] Probst | Structural differences in centromeric heterochromatin are spatially reconciled on fertilisation in the mouse zygote[END_REF] . This demonstrates the fact that chromocenter organization is dynamic during development; it is established early during development, and depending on the cell type and the state of differentiation different numbers of chromocenters are formed.

Maintenance of heterochromatin marks throughout the cell cycle

Studies in S. pombe have been instrumental in gaining insight into the mechanisms in place in order to ensure the maintenance of heterochromatic features during the cell cycle. At fission yeast centromeres, heterochromatin forms at repetitive DNA elements, referred to as pericentric dg and dh repeats (Figure 1) and is essential for proper loading of the centromere specific H3-variant Cenp-A and kinetochore assembly [START_REF] Folco | Heterochromatin and RNAi are required to establish CENP-A chromatin at centromeres[END_REF] . Nucleosomes at the pericentric repeats are hypoacetylated and enriched in H3K9me2, a histone mark catalyzed by the histone methyltransferase Clr4 and bound by chromodomain-containing proteins Chp1 and the HP1 homologues Swi6 and Chp2, which have the capacity to multimerize. Clr4 itself can bind H3K9me2 via its proper chromodomain and methylate neighboring nucleosomes thereby creating additional binding sites for HP1 proteins. Hypoacetylation of histones is ensured by histone de-acetylases recruited by the HP1 homologs [START_REF] Yamada | The nucleation and maintenance of heterochromatin by a histone deacetylase in fission yeast[END_REF] . Together these chromatin marks limit Pol II occupancy and ensure the transcriptionally repressed state at pericentric heterochromatin. Interestingly, expression of Forward and Reverse pericentric transcripts is differentially regulated: Forward transcripts accumulate both in strains deficient in Swi6 and RNAi pathway components, while Reverse transcripts accumulate only in RNAi mutants. This suggests that expression of the Forward strand is under transcriptional and the Reverse strand under post-transcriptional control [START_REF] Volpe | Regulation of heterochromatic silencing and histone H3 lysine-9 methylation by RNAi[END_REF] . However, also in wild type cells, an exquisite equilibrium exists between transcriptional silencing and activating activity. Indeed, low-level, bidirectional transcription of the pericentric repeats is actually required to ensure the stable maintenance of molecular heterochromatin marks during the cell cycle. In wild type cells, transcription of pericentric repeats peaks in S-phase. Transcription during this specific time window of the cell cycle could possibly occur as a consequence of the phospho-methyl switch in mitosis which leads to Swi6 delocalization and the dilution of repressive histone marks at the daughter strands [START_REF] Kloc | RNA Interference Guides Histone Modification during the S Phase of Chromosomal Replication[END_REF][START_REF] Chen | Cell cycle control of centromeric repeat transcription and heterochromatin assembly[END_REF] , Figure 2A. Pericentric transcripts are thought to form double stranded RNA or hairpin structures [START_REF] Djupedal | Analysis of small RNA in fission yeast; centromeric siRNAs are potentially generated through a structured RNA[END_REF] , which are further processed by the RNAi machinery [START_REF] Volpe | Regulation of heterochromatic silencing and histone H3 lysine-9 methylation by RNAi[END_REF] .

The resulting small RNAs (siRNAs) that accumulate transiently in S-phase then target the RNA-induced transcriptional silencing (RITS) complex [START_REF] Verdel | RNAi-mediated targeting of heterochromatin by the RITS complex[END_REF] and Clr4 to re-establish heterochromatin features. The RNA-directed RNA polymerase complex (RDRC) ensures amplification of the siRNA signal. However, in the absence of the RDRC, small degradation products of pericentric transcripts can mediate H3K9 methylation [START_REF] Halic | Dicer-independent primal RNAs trigger RNAi and heterochromatin formation[END_REF] . Taken together, these data suggest that the continuous presence of heterochromatin at pericentric repeats actually requires re-establishment of heterochromatin features after each round of replication.

Collectively, work in fission yeast has revealed a critical role for pericentric transcripts and the RNAi-dependent pathway in cell cycle maintenance of heterochromatin in locations where heterochromatin was already established.

In mice, the AT-rich major satellites are arranged in arrays up to several megabases in length 20 . These highly methylated repeats that are organized into hypoacetylated chromatin show high H3K9 methylation that depends on the Suv39h1/2 methyltransferases and are bound by the three isoforms of HP1. Spreading of heterochromatin features is thought to exploit a selfsustaining loop mechanism through HP1 self-association and interaction with the H3K9HMTase, Suv39h1/2, which can add more H3K9me [START_REF] Aagaard | Functional mammalian homologues of the Drosophila PEV-modifier Su(var)3-9 encode centromereassociated proteins which complex with the heterochromatin component M31[END_REF][START_REF] Maison | HP1 and the dynamics of heterochromatin maintenance[END_REF] , and the capacity of HP1 to recruit de novo methyltransferase activity [START_REF] Lehnertz | Suv39h-mediated histone H3 lysine 9 methylation directs DNA methylation to major satellite repeats at pericentric heterochromatin[END_REF] . Observations of reduced DNA methylation at pericentric repeats observed in Suv39h mutant cells supports this model [START_REF] Lehnertz | Suv39h-mediated histone H3 lysine 9 methylation directs DNA methylation to major satellite repeats at pericentric heterochromatin[END_REF] . In both S. pombe and mouse pericentromeres, H3K9HMTase, working as the writer, together with the reader, HP1/Swi6, propagate the heterochromatin state. However, a role for RNA-based mechanisms in the maintenance of heterochromatin features has remained difficult to establish in mammals [START_REF] Kanellopoulou | Dicerdeficient mouse embryonic stem cells are defective in differentiation and centromeric silencing[END_REF][START_REF] Murchison | Characterization of Dicer-deficient murine embryonic stem cells[END_REF][START_REF] Cobb | T cell lineage choice and differentiation in the absence of the RNase III enzyme Dicer[END_REF] . In mammals, pericentric repeats are transcribed [START_REF] Lehnertz | Suv39h-mediated histone H3 lysine 9 methylation directs DNA methylation to major satellite repeats at pericentric heterochromatin[END_REF][START_REF] Rudert | Transcripts from opposite strands of gamma satellite DNA are differentially expressed during mouse development[END_REF] and regulated in a cellcycle dependent manner [START_REF] Lu | Proliferation-dependent and cell cycle regulated transcription of mouse pericentric heterochromatin[END_REF] , yet, during the cell cycle, transcripts accumulate in mitosis and in G1/S. In contrast, in fission yeast, transcription ceases in mid S-phase when pericentric heterochromatin replicates, and therefore, the analogy to fission yeast cannot be extended to include the exact timing during the cell cycle. Furthermore, if RNAi-based mechanisms play a role, as reported for the transposon-control in oocytes [START_REF] Stein | RNAi: mammalian oocytes do it without RNA-dependent RNA polymerase[END_REF] the mechanism possibly operates without amplification of the siRNA signal, since a properly defined RNA-directed RNA polymerase still remains to be characterized. It is also possible that maintenance of molecular heterochromatin marks through chromosome duplication may be ensured by the mutual reinforcement of mechanisms that maintain DNA methylation and inheritance of modified histones as discussed [START_REF] Probst | Epigenetic inheritance during the cell cycle[END_REF] . In this context, the parental histone mark would be sufficient to reinitiate the whole process of heterochromatinization by recruiting HP1 and in turn the H3K9me methyltransferase. Furthermore, the maintenance methyltransferase Dnmt1 together with the SRA-domain protein Np95 31 could also ensure methylation maintenance and recruit histone deacetylase activity necessary to remove acetylation marks on the newly incorporated histones [START_REF] Fuks | DNA methyltransferase Dnmt1 associates with histone deacetylase activity[END_REF][START_REF] Unoki | ICBP90, an E2F-1 target, recruits HDAC1 and binds to methyl-CpG through its SRA domain[END_REF] . Np95 shows affinity to hemi-methylated DNA [START_REF] Arita | Recognition of hemimethylated DNA by the SRA protein UHRF1 by a base-flipping mechanism[END_REF][START_REF] Avvakumov | Structural basis for recognition of hemi-methylated DNA by the SRA domain of human UHRF1[END_REF] and interacts with H3 tails methylated at position K9 [START_REF] Rottach | The multidomain protein Np95 connects DNA methylation and histone modification[END_REF] , thereby linking DNA and histone modifications to ensure Np95's function during replication of pericentric domains.

Notably, the chromatin assembly machinery at the replication fork contributes in a decisive manner to heterochromatin maintenance. The Chromatin Assembly Factor 1 (CAF-1) deposits new H3-H4 in a DNA-synthesis coupled manner [START_REF] Verreault | Nucleosome assembly by a complex of CAF-1 and acetylated histones H3/H4[END_REF] . New histone H3 bound to CAF-1 mono-methylated by SetDB1 38 is a favoured substrate for Suv39h, which can then establish H3K9me3 and HP1 binding [START_REF] Loyola | PTMs on H3 variants before chromatin assembly potentiate their final epigenetic state[END_REF] . How CAF-1 possibly contributes to HP1 inheritance at the replication fork by transferring parental HP1 to the replicated daughter chromatin remains to be established [START_REF] Quivy | A CAF-1 dependent pool of HP1 during heterochromatin duplication[END_REF][START_REF] Quivy | The HP1-p150/CAF-1 interaction is required for pericentric heterochromatin replication and S-phase progression in mouse cells[END_REF] , Figure 2B. The recently identified chromatin-remodelling factor, SMARCAD1, which interacts with PCNA, is likely to participate in these dynamic events connecting histone methylation and HP1 recruitment [START_REF] Rowbotham | Maintenance of Silent Chromatin through Replication Requires SWI/SNF-like Chromatin Remodeler SMARCAD1[END_REF] . In summary, in mammalian pericentric regions, a combination of factors act in concert to ensure maintenance of heterochromatin: DNA methylation maintenance and histone de-acetylation are coupled at the replication fork, as well as histone deposition with further post-translational modifications.

Thus, given the robustness of these multiple layers involved in the cell cycle maintenance of heterochromatin in mammals, an RNA-dependent mechanism may either: not be necessary at all, or may not be possible to reveal under normal conditions. Interestingly, plant heterochromatin, similar to mammalian heterochromatin, is enriched in both DNA methylation and H3K9methylation. These two chromatin marks are tightly interconnected; histone methylation depends on DNA methylation [START_REF] Soppe | DNA methylation controls histone H3 lysine 9 methylation and heterochromatin assembly in Arabidopsis[END_REF][START_REF] Tariq | Erasure of CpG methylation in Arabidopsis alters patterns of histone H3 methylation in heterochromatin[END_REF] and vice versa [START_REF] Jackson | Control of CpNpG DNA methylation by the KRYPTONITE histone H3 methyltransferase[END_REF] . DNA methylation is critical to silence transposable elements and repeated sequences [START_REF] Rigal | A "mille-feuille" of silencing: Epigenetic control of transposable elements[END_REF] and its maintenance is ensured by active mechanisms operating at the replication fork. However, plants lack HP1-like proteins at heterochromatin that could allow spreading of heterochromatin marks in a self-reinforcing loop as is proposed to occur in fission yeast and mammals. Instead, small RNA molecules processed from satellite transcripts or transposable elements do reinforce heterochromatin at pericentromeres in Arabidopsis [START_REF] May | Differential regulation of strand-specific transcripts from Arabidopsis centromeric satellite repeats[END_REF] . While small RNA molecules do not target histone methylation activity as in yeast, they direct de novo DNA methylation activity to homologous sequences through a mechanism termed RNA-directed DNA methylation (RdDM, [START_REF] Wassenegger | RNA-directed de novo methylation of genomic sequences in plants[END_REF] . In plants, it remains to be established whether transcription of pericentric repeats is particularly increased following passage of the replication fork in order to reinforce DNA methylation maintenance. These examples illustrate that while most of the heterochromatin features are shared among different species, certain heterochromatin marks are specific to a given organism, as is the way their heterochromatin maintenance is coordinated with other chromatin modifications.

Pericentric heterochromatin dynamics during reproductive development

During mammalian germ cell development and gametogenesis genome organization undergoes important changes, and this is clearly visible at pericentric heterochromatin domains. While in mature oocytes, pericentric satellite repeats retain most of the somatic heterochromatin marks with the exception of HP1α, in sperm, the majority of histones are replaced by protamines or sperm-specific histone variants [START_REF] Govin | Pericentric heterochromatin reprogramming by new histone variants during mouse spermiogenesis[END_REF] , therefore requiring de novo establishment of most heterochromatic features. Upon fertilization, protamines and the sperm-specific histone variants are removed [START_REF] Wu | Testis-specific histone variants H2AL1/2 rapidly disappear from paternal heterochromatin after fertilization[END_REF] and a nucleosomal organization is restored.

As histone assembly occurs in a DNA-synthesis independent manner, only H3.3 and not the replicative H3 variants are incorporated [START_REF] Van Der Heijden | Asymmetry in histone H3 variants and lysine methylation between paternal and maternal chromatin of the early mouse zygote[END_REF] . This accentuates the chromatin differences at pericentric satellite repeats between the two parental genomes as paternal pericentric domains also lack H3K9me3, and H4K20me3 51 (Figure 3). Instead, during the first cleavage stages, alternative repressive chromatin factors, the Polycomb group binding proteins, are enriched at paternal pericentric satellites [START_REF] Puschendorf | PRC1 and Suv39h specify parental asymmetry at constitutive heterochromatin in early mouse embryos[END_REF] . In the zygote, the paternal genome undergoes active global DNA demethylation [START_REF] Mayer | Demethylation of the zygotic paternal genome[END_REF] , however some residual DNA methylation is retained at pericentric repeats [START_REF] Santos | Dynamic reprogramming of DNA methylation in the early mouse embryo[END_REF] where it could play a role in the establishment of other heterochromatin features.

During the early cleavage stages, the spatial organization of pericentric satellite repeats is particularly dynamic; the pericentric domains reorganize from the gamete-specific organization to form rings around the nucleolar precursor bodies [START_REF] Probst | Structural differences in centromeric heterochromatin are spatially reconciled on fertilisation in the mouse zygote[END_REF] . These ring-structures progressively reorganize into chromocenters during the 2-cell stage. During this particular developmental time window, transcripts from pericentric major repeats undergo a burst in transcription reaching levels far above the amount observed in cycling embryonic fibroblasts 55 . While the repeats were known to be transcribed bi-directionally [START_REF] Rudert | Transcripts from opposite strands of gamma satellite DNA are differentially expressed during mouse development[END_REF] , unexpectedly, their transcription and spatial localization is controlled in a strand-specific manner during early embryonic development [START_REF] Probst | A strandspecific burst in transcription of pericentric satellites is required for chromocenter formation and early mouse development[END_REF] . The Forward strand is transcribed during the S-phase of the second cell cycle and peaks in late S/G2 followed by upregulation of the Reverse strand.

Transcription from both strands is then strongly downregulated at the 4-cell stage, when chromocenters have formed, which are faithfully propagated during the following cleavage stages (Figure 3). Interference with transcripts during the 2-cell stage results in embryos that fail to form chromocenters and arrest their development in the G2 phase of the second cell cycle [START_REF] Probst | A strandspecific burst in transcription of pericentric satellites is required for chromocenter formation and early mouse development[END_REF] . These experiments illustrate a critical necessity for pericentric transcription during heterochromatin establishment in early embryonic development. How transcription of the two strands is regulated in time and space and which promoters and transcription factors are necessary to drive satellite transcription specifically during this developmental transition are still unresolved questions. Differences in timing of transcription and differential regulation of subsets of repeats by distinct silencing pathways have been observed in fission yeast [START_REF] Volpe | Regulation of heterochromatic silencing and histone H3 lysine-9 methylation by RNAi[END_REF][START_REF] Kloc | RNA Interference Guides Histone Modification during the S Phase of Chromosomal Replication[END_REF] and plants [START_REF] May | Differential regulation of strand-specific transcripts from Arabidopsis centromeric satellite repeats[END_REF] . Investigating the strand-specific expression of major satellite sequences in mutant backgrounds for chromatin modifiers may help to better understand their mode of regulation.

The exact mechanism of how pericentric transcripts affect chromocenter formation and development as well as the role of the two respective transcripts remain to be dissected, however, recent data have allowed speculation about their mode of action. Fertilized eggs expressing a mutated form of H3.3 (K27R) arrest at the 2-cell stage and show HP1β mislocalization. This arrest, in addition to proper HP1β localization can be rescued by injection of double-stranded, but not single stranded major RNA, implying processing of double stranded RNA molecules and a role for small RNA based mechanisms [START_REF] Santenard | Heterochromatin formation in the mouse embryo requires critical residues of the histone variant H3.3[END_REF] . Alternatively, pericentric transcripts could function in recruiting heterochromatin components or modifiers, similar to the long non-coding transcript Air, involved in imprinting that interacts with the histone methyltransferase G9a [START_REF] Nagano | The Air noncoding RNA epigenetically silences transcription by targeting G9a to chromatin[END_REF] . An intriguing candidate to interact with pericentric transcripts is HP1, which has RNA-binding capacity [START_REF] Muchardt | Coordinated methyl and RNA binding is required for heterochromatin localization of mammalian HP1alpha[END_REF][START_REF] Maison | Higher-order structure in pericentric heterochromatin involves a distinct pattern of histone modification and an RNA component[END_REF] . HP1 indeed binds major satellite RNA, however, it does so only when sumoylated and it also shows specificity for Forward transcripts.

Interestingly, sumoylation of HP1 enhances its capacity to be recruited de novo to pericentric satellites lacking H3K9me3 and HP1 [START_REF] Maison | SUMOylation promotes de novo targeting of HP1alpha to pericentric heterochromatin[END_REF] .

In short, the specific organization of DNA in the sperm and the loss of somatic heterochromatin marks entail a challenge for the early embryo to cope with two distinct states of pericentric heterochromatin in the same nucleus. Early development therefore necessitates setting most somatic heterochromatin marks at paternal and requires consolidation of maternal heterochromatin domains; it may also be one of the rare situations in which de novo formation of constitutive heterochromatin can be studied in mammals. The demonstration of a critical role of pericentric satellite transcripts may be one of the first steps towards a better understanding of the mechanisms implicated in mammalian heterochromatin establishment. , in both mouse and Arabidopsis, heterochromatin organizes in several chromocenters [START_REF] Guenatri | Mouse centric and pericentric satellite repeats form distinct functional heterochromatin[END_REF][START_REF] Fransz | Interphase chromosomes in Arabidopsis are organized as well defined chromocenters from which euchromatin loops emanate[END_REF] .

Pericentric heterochromatin nucleates at repetitive elements. In S. pombe centromeres, these consist of the outer dg and dh repeats present on each chromosome. In Arabidopsis, heterochromatin forms on 180bp satellite repeats interspersed with Athila retrotransposons and their remnants as well as various types of transposons that flank these repeats [START_REF] May | Differential regulation of strand-specific transcripts from Arabidopsis centromeric satellite repeats[END_REF][START_REF] Kumekawa | The size and sequence organization of the centromeric region of arabidopsis thaliana chromosome 5[END_REF] . In mouse, the pericentric domains consist of large uninterrupted arrays of AT-rich 234bp-long major satellite repeats [START_REF] Vissel | Mouse major (gamma) satellite DNA is highly conserved and organized into extremely long tandem arrays: implications for recombination between nonhomologous chromosomes[END_REF] , which can extend to more than 2 Mb per chromosome.

Figure 2:

Pericentric heterochromatin maintenance during the cell cycle in S. pombe and mouse

A.

Fission yeast spends most of its lifetime in G2 phase, during which time, H3K9me2, at pericentric repeats constitutes a binding site for the chromodomain protein Swi6, the HP1 homolog. Phosphorylation of serine 10 of H3 (phospho-methyl switch) during mitosis (M) induces loss of Swi6 binding. This, together with the subsequent dilution of repressive histone modifications upon replication of pericentric repeats in early S-phase creates a small time window during which Pol II gains access and pericentric transcripts accumulate. These transcripts are processed by the RNAi machinery into small RNAs (siRNAs) that target the methyltransferase Clr4 (the homolog of mammalian SUV39H) to pericentric repeats. After deacetylation of the newly incorporated histones, Clr4 then remethylates K9, Swi6 is recruited, and transcriptional silencing is reestablished.

B.

In G2, pericentric heterochromatin is enriched in DNA methylation and H3K9me3 that is bound by HP1 in mice. While HP1 binding is affected by the phospho-methyl switch in mitosis [START_REF] Fischle | Regulation of HP1-chromatin binding by histone H3 methylation and phosphorylation[END_REF][START_REF] Hirota | Histone H3 serine 10 phosphorylation by Aurora B causes HP1 dissociation from heterochromatin[END_REF] , the extent of HP1 disruption and the mechanisms of its restoration in G1 phase are as yet unknown. Transcripts from pericentric satellites have been observed in mitosis and in G1/early S-phase, however major satellite transcription is downregulated during mid Sphase when pericentric domains replicate [START_REF] Lu | Proliferation-dependent and cell cycle regulated transcription of mouse pericentric heterochromatin[END_REF] . Coupled to DNA replication, histone de-acetylase activity is recruited by the DNA clamp PCNA and the DNA methyltransfease Dnmt1 [START_REF] Fuks | DNA methyltransferase Dnmt1 associates with histone deacetylase activity[END_REF] to deacetylate the newly incorporated histones. Dnmt1 restores DNA methylation. PCNA further recruits the Chromatin Assembly Factor 1 (CAF1). CAF1 can assemble H3 already monomethylated at K9 by SetDB1 [START_REF] Loyola | The HP1alpha-CAF1-SetDB1-containing complex provides H3K9me1 for Suv39-mediated K9me3 in pericentric heterochromatin[END_REF] into nucleosomes, but it also plays an important role in HP1 dynamics by transferring parental HP1 to daughter strands [START_REF] Quivy | The HP1-p150/CAF-1 interaction is required for pericentric heterochromatin replication and S-phase progression in mouse cells[END_REF] .

Figure 3.

Heterochromatin organization and transcription during early cleavage stages in mouse

Mouse sperm DNA tightly compacts around a core of pericentric repeats (blue) in the centre of the sperm [START_REF] Haaf | Higher order nuclear structure in mammalian sperm revealed by in situ hybridization and extended chromatin fibers[END_REF] . After fertilization, the egg completes meiosis II and extrudes one haploid genome in form of the polar body (PB). Two independent pronuclei form, in which pericentric domains reorganize into ring-structures around the nucleolus-like bodies [START_REF] Probst | Structural differences in centromeric heterochromatin are spatially reconciled on fertilisation in the mouse zygote[END_REF] . Maternal (pink)

and paternal (blue) pericentric domains are epigenetically distinct in the zygote [START_REF] Probst | Structural differences in centromeric heterochromatin are spatially reconciled on fertilisation in the mouse zygote[END_REF][START_REF] Santos | Dynamic chromatin modifications characterise the first cell cycle in mouse embryos[END_REF][START_REF] Puschendorf | PRC1 and Suv39h specify parental asymmetry at constitutive heterochromatin in early mouse embryos[END_REF] and retain asymmetry up to the 8-cell stage [START_REF] Puschendorf | PRC1 and Suv39h specify parental asymmetry at constitutive heterochromatin in early mouse embryos[END_REF][START_REF] Merico | Epigenomic differentiation in mouse preimplantation nuclei of biparental, parthenote and cloned embryos[END_REF] . Paternal pericentric repeats lack H3K9me3 and

Figure 1 :

 1 Figure 1: Pericentric heterochromatin clustering

  H4K20me3 and either do not, or only weakly bind HP1β[START_REF] Santos | Dynamic chromatin modifications characterise the first cell cycle in mouse embryos[END_REF][START_REF] Puschendorf | PRC1 and Suv39h specify parental asymmetry at constitutive heterochromatin in early mouse embryos[END_REF] . During the 2-cell stage chromocenters progressively form and are then propagated during further cleavage stages.The first transcription of the zygotic genome (minor zygotic gene activation, ZGA) takes place at the end of first cell cycle, followed by the major ZGA in the 2-cell stage. The formation of chromocenters during the 2-cell stage coincides with a transcriptional burst of pericentric satellites ( 55 see text).

Acknowledgments:

We apologize for not quoting all of our colleagues for their contributions due to space limitations. We thank C. Maison for Figure 1A and Alysia vandenBerg for critical reading of the manuscript. G.A. is supported by la Ligue Nationale contre le Cancer (Equipe labellisée Ligue 2010), PIC Programs, the European Commission Network of Excellence EpiGeneSys (HEALTH-F4-2010-257082), the European Commission ITN FP7-PEOPLE-2007-215148 "Image DDR" and FP7-PEOPLE-2008-238176 "Nucleosome 4D", ERC Advanced Grant 2009-AdG_20090506 "Eccentric", the European Commission large-scale integrating project FP7_HEALTH-2010-259743 "MODHEP", ANR "ECenS" ANR-09-BLAN-0257-01, ANR "ChromaTin" ANR-10-BLAN-1326-03 and INCa "GepiG".


Stanis Law

A Raczyński

Satoru Fukayama

Emmanuel Vincent

Melody harmonisation with interpolated probabilistic models

Keywords: Dynamic Bayesian Networks, melody harmonisation, model interpolation

Automatic melody harmonisation aims to create a matching chordal accompaniment to a given monophonic melody. Several methods have been proposed to this aim, which are generally based on musicological expertise or on unsupervised probabilistic modelling. Among the latter category of methods, most systems use the generative hidden Markov model (HMM), in which the chords are the hidden states and the melody is the observed output. Relations to other variables, such as the tonality and scale or the metric structure, are handled by training multiple HMMs or are often simply ignored. In this paper, we propose a means of combining multiple probabilistic models of various musical variables into a versatile harmonisation system by means of model interpolation. The result is a joint model belonging to the class of discriminative models, which in recent years have proven to be capable of outperforming generative models in many tasks. We first evaluate our models in terms of their normalized negative log-likelihood, or cross-entropy. We observe that log-linear interpolation offers lower cross-entropy than linear interpolation and that combining several models by means of log-linear interpolation lowers the cross-entropy compared to the best of the component models. We then perform a series of harmonisation experiments and show that the proposed log-linearly interpolated model offers higher chord root accuracy than a reference musicological rule-based harmoniser by up to 5% absolute.

Harmonisation de mélodies par interpolation de modèles probabilistes

Résumé : L'harmonisation automatique de mélodies vise à créer une suite d'accords accompagnant une mélodie donnée. Plusieurs méthodes ont été proposées dans ce but, généralement basées sur des règles musicologiques expertes ou sur une modélisation probabiliste non supervisée. Parmi cette dernière catégorie de méthodes, la plupart utilisent un modèle de Markov caché (MMC) génératif, dont les états cachés sont les accords et l'observation la mélodie. Les dépendances aux autres variables telles que la tonalité ou la structure métrique sont modélisées par des MMCs multiples ou simplement ignorées. Dans ce papier, nous proposons un moyen de combiner plusieurs modèles probabilistes de différentes variables musicales par le biais d'une interpolation de modèles. Cela aboutit à un modèle combiné appartenant à la catégorie des modèles discriminants, dont il a été démontré ces dernières années qu'ils dépassent la performance des modèles génératifs pour de nombreuses tâches.

Nous évaluons d'abord nos modèles en terme de l'opposé de leur log-vraisemblance normalisée, ou entropie croisée. Nous observons que l'interpolation log-linéaire diminue l'entropie croisée par rapport à l'interpolation linéaire et que la combinaison de plusieurs modèles par interpolation log-linéaire diminue l'entropie croisée par rapport au meilleur modèle individuel. Nous effectuons ensuite un ensemble d'expériences d'harmonisation et montrons que le modèle par interpolation log-linéaire proposé améliore la précision d'estimation de la fondamentale des accords de 5% dans l'absolu par rapport un algorithme d'harmonisation de référence basé sur des règles musicologiques expertes.

Mots-clés : Réseaux bayésiens dynamiques, harmonisation de mélodies, interpolation de modèles probabilistes

Introduction

Automatic melody harmonisation is the process of determining the most musically suitable chordal accompaniment to a given monophonic melody. It is an important part of musical composition and so it is a common exercise in all music composition classes, where it typically involves determining a four-part harmony, i.e., the movement of four voices, namely: soprano, alto, tenor and bass. The task can be either to find the harmonization for a given melody performed by the soprano voice, or to find the three other voices for a given bass line (unfigured bass), often with some chordal information given (figured bass). In this work, however, we focus on harmonisation in the more narrow senseof generating a sequence of background chords matching a given melody, which can be played on supporting instruments, e.g., on a guitar or a piano. It is simpler than four-part harmonisation because one does not need to determine the exact movements of the voices-it does not include inversions, added and removed tones, etc., but only the root pitch (C, C♯, etc.) and the chord type (major, minor, etc.). The melody together with the chord labels are typically referred to as lead sheets, an example of which is shown in Fig. 1.

Harmonisation is a necessary step in most algorithmic composition methods, which typically involve generating a melodic line first, and afterwards supplementing it with an accompaniment. For example, in their Orpheus automatic composition system, Fukayama et al. compose a melody based on constraints resulting from the tonal accent in the lyrics and from basic musicological rules, and then compose the chordal accompaniment using a collection of accompaniment patterns [START_REF] Fukayama | Automatic song composition from the lyrics exploiting prosody of the Japanese language[END_REF]. Harmonisation has also been explored as an easy way to create polyphonic ring-tones from simple melodies in the i-Ring ring-tone harmonisation system from [START_REF] Lee | i-Ring: A system for humming transcription and chord generation[END_REF]. Furthermore, automatic harmonisation has recently received significant commercial interest as an easy way for non-musicians to create well-sounding music based on simple melodies. The most well-known implementations are the MySong software developed in cooperation with Microsoft [START_REF] Simon | MySong: automatic accompaniment generation for vocal melodies[END_REF] and the commercial software package Band-in-a-Box (BIAB) [START_REF]Band-in-a-box[END_REF]. Both are designed as tools for non-professional musicians, or even non-musicians, to create songs with instrumental accompaniment by singing a melody into a microphone. All of the above methods perform a lead sheet-like harmonisation.

The history of automatic harmonisation is much older, however. Some of the earliest attempts at harmonisation were made by Steels, who proposed a heuristic search approach in [START_REF] Steels | Learning the craft of musical composition[END_REF], and by Ebcioglu, who proposed a rule-based system targeting Bach's chorales [START_REF] Ebcioglu | An expert system for chorale harmonization[END_REF][START_REF] Ebcioglu | An expert system for harmonizing four-part chorales[END_REF]. Another rule-based harmoniser based on musicological expertise-called Harmonic Analyzer -was more recently developed by Temperley and Sleator in [START_REF] Temperley | Modeling meter and harmony: A preferencerule approach[END_REF][START_REF] Temperley | Harmonic Analyzer[END_REF]. Other harmonisation methods have also been explored. Phon-Amnuaisuk and Wiggins developed a prototype system based on genetic algorithms with knowledge-rich structures in [START_REF] Phon-Amnuaisuk | The four-part harmonisation problem: a comparison between genetic algorithms and a rule-based system[END_REF]. Constraint satisfaction systems were investigated by Pachet and Roy in [START_REF] Pachet | Musical harmonization with constraints: A survey[END_REF] and sequential neural networks by Gang et al. in [START_REF] Gang | Tuning a neural network for harmonizing melodies in real-time[END_REF]. Neural networks were also used by Cunha and Ramalho in their hybrid neural-musicological real-time chord prediction system in [START_REF] Cunha | An intelligent hybrid model for chord prediction[END_REF]. Another hybrid harmonisation system that combines Hidden Markov Models (HMMs) with a set of heuristic rules for rapid training was proposed by Chuan and Chew in [START_REF] Chuan | A hybrid system for automatic generation of style-specific accompaniment[END_REF]. Among these, only [START_REF] Temperley | Harmonic Analyzer[END_REF], [START_REF] Gang | Tuning a neural network for harmonizing melodies in real-time[END_REF], [START_REF] Cunha | An intelligent hybrid model for chord prediction[END_REF] and [START_REF] Chuan | A hybrid system for automatic generation of style-specific accompaniment[END_REF] are lead sheet-like harmonisation systems.

Among the most flexible systems are those based on unsupervised probabilistic modelling, which typically utilise HMMs. Single-HMM approaches include the MySong software and an implementation based on MySong from [START_REF] Chuan | A comparison of statistical and rule-based models for stylespecific harmonization[END_REF], as well as the i-Ring ring-tone harmonisation system from [START_REF] Lee | i-Ring: A system for humming transcription and chord generation[END_REF]. A slightly more complicated four-part harmonisation approach using a dual HMM (one for harmonisation and another one for ornamentation) was proposed by Allan and Williams in [START_REF] Allan | Harmonising chorales by probabilistic inference[END_REF]. Later, in [START_REF] Paiement | Probabilistic melodic harmonization[END_REF] Paiment et al. proposed a very sophisticated, multi-level graphical model for modelling chord movement with respect to the melody, which is capable of modelling long-term relationships between chords, as opposed to HMMs that are only capable of modelling short-term dependencies. However, being a non-dynamic graphical network, their model is limited to fixed-length songs (of exactly 16 measures). By contrast with musicological rule-based methods, which require careful formulation and application of harmonisation rules for particular genres (e.g., classical or jazz), these probabilistic methods aim to automatically infer those rules from a corpus of example data and are therefore applicable to all genres, even when musicological expertise is not available.

All of the above probabilistic harmonisation systems are however limited to modelling the relation between a single hidden layer (chord sequence) and a single observed layer (melody), without any explicit mechanism to include models of other relevant musical quantities, such as the key and the current tonal centre, the rhythm and the musical accent, or the genre, period and the composer. One can use multiple HMMs corresponding to, e.g., different genres or keys as in [START_REF] Simon | MySong: automatic accompaniment generation for vocal melodies[END_REF], but with many variables this approach quickly suffers from over-fitting. In this paper we propose to build versatile chord models by interpolating between a collection of simpler sub-models using linear or log-linear interpolation.

This paper is organised as follows. Section 2 explains the proposed modelling and training approach and Section 3 gives details about the particular submodels used in our experiments. The experimental set-up and the results are described in Section 4. Finally, the conclusion is given in Section 5.

General approach

When figuring out the accompaniment, one needs to keep in mind the basic rules of tonal music: the tonality (everything in tonal music happens with respect to the tonic), the chord progressions (certain chord progressions are more natural and pleasant, e.g., progressions corresponding to the circle of fifths, progressions descending by thirds and common cadences), and of course harmonic compatibility with the melody. In the generative, HMM-based systems, the current chord is the underlying state C t , while the observation is the melody M t . Chord progression is modelled with a Markov chain P(C t |C t-1 ) and the observed melody by a multinomial distribution conditioned on the system's state P(M t |C t ) (see the top of Fig. 2).

Model structure

We propose a more flexible way of developing probabilistic harmonisation models, in which the time-varying tonality T t , as well as other musical variables can be explicitly taken into account. We propose a discriminative model, in which the chords are modelled conditionally to all other variables (see bottom of Fig. 2):

P(C t |C 1:t-1 , X 1:t ), (1)
where 1 : t denotes a range of time indices from the first to the current time frame and X is a set of other musical variables, such as the melody, the tonal centre, the metrical accent, the style or genre, etc. Discriminative models have been found to outperform generative models in many fields [START_REF] Jebara | Machine Learning: Discriminative and Generative[END_REF], such as speech recognition [START_REF] Rathinavelu | The trended HMM with discriminative training for phonetic classification[END_REF][START_REF] Woodland | Large scale discriminative training of hidden Markov models for speech recognition[END_REF], machine translation [START_REF] Och | Discriminative training and maximum entropy models for statistical machine translation[END_REF], text classification [START_REF] Rennie | Improving multiclass text classification with the support vector machine[END_REF] or genomics [START_REF] Jaakkola | A discriminative framework for detecting remote protein homologies[END_REF]. This conditional multinomial distribution has too many parameters to be used in practice, hence we approximate it by interpolating between multiple sub-models P i involving a different subset of conditioning variables A i,t ⊂ {C 1:t-1 , X 1:t }. The interpolation can be linear,

P(C t |C 1:t-1 , X 1:t ) = I i=1 a i P i (C t |A i,t ), (2)
with

I i=0 a i = 1, (3)
or log-linear,

P(C t |C 1:t-1 , X 1:t ) = Z -1 I i=1 P i (C t |A i,t ) bi , ( 4
)
where I is the number of sub-models, a i ≥ 0 and b i ≥ 0 for i = 1, . . . , I are the interpolation coefficients and

Z = Ct I i=1 P i (C t |A i,t ) bi (5)
is a normalizing factor depending on A i,t . For example, we will consider in the following:

A 1,t = {C t-1 }, A 2,t = {T t } and A 3,t = {M t }.
Linear [START_REF] Jelinek | Interpolated estimation of Markov source parameters from sparse data[END_REF] and log-linear [START_REF] Klakow | Log-linear interpolation of language models[END_REF] interpolation have been previously used in the context of natural (spoken) language modelling to combine models with different temporal spans (n-grams with different values of n). Here we have generalized this approach to interpolate between sub-models conditioned on different musical variables.

Smoothing

Although the sub-models P i now have fewer parameters, over-fitting issues may still arise due to data sparsity, so the above equations are not directly usable. In order to address these issues, each of the sub-models must be smoothed [START_REF] Zhai | A study of smoothing methods for language models applied to information retrieval[END_REF]. In this study we perform smoothing by combining each sub-model with the prior chord distribution and a uniform distribution. In the case of linear interpolation, the smoothing is applied to the interpolated model:

P(C t |C 1:t-1 , X 1:t ) = I i=1 a i P i (C t |A i,t ) + αP(C t ) + β, (6)
with

α + β + I i=1 a i = 1. (7)
In the case of the log-linear interpolation, each model is smoothed separately before combining them:

P(C t |C 1:t-1 , X 1:t ) = Z -1 I i=1 (γ i P i (C t |A i,t ) + δ i P(C t ) + ǫ i ) bi , (8)
with γ i + δ i + ǫ i = 1 (9)
for all i and

Z = Ct I i=1 (γ i P i (C t |A i,t ) + δ i P(C t ) + ǫ i ) bi . (10)

Training

The proposed models are trained on two disjoint sets of example data called training set and validation set. The sub-models P i and the prior distribution P(C t ) are first trained in the maximum likelihood (ML) sense on the training set by counting occurrences [START_REF] Zhai | A study of smoothing methods for language models applied to information retrieval[END_REF]. The interpolation coefficients a i or b i and the smoothing coefficients α and β or γ i , δ i and ǫ i are then jointly trained on the validation set according to one of two possible training objectives. Classical generative training is achieved by estimating the interpolation and smoothing coefficients in the ML sense on the validation set. Because the loglikelihood is convex [START_REF] Klakow | Log-linear interpolation of language models[END_REF], any optimization algorithm can be used. In the following, we have used a non-negatively constrained limited-memory Broyden-Fletcher-Goldfarb-Shanno (BFGS) method (a quasi-Newton optimisation), built into the GNU R environment as the optim() function [START_REF] Development | R: A Language and Environment for Statistical Computing[END_REF].

Even though the likelihood is a common evaluation measure for statistical models, higher likelihood does not always translate into better performance for the considered application. For example, it is well known in the field of automatic speech recognition that the likelihood of language models can sometimes be improved without effect on the word error rate [START_REF] Chen | Evaluation metrics for language models[END_REF]. For that reason, we alternatively propose to perform discriminative training by estimating the interpolation and smoothing coefficients so as to maximize the chord root note accuracy on the validation set, which is the main evaluation metric used in the harmonisation experiments in Subsection 4.2. Because this metric is not differentiable, gradient-based methods cannot be used so we have used the following multi-step brute-force search. First, all smoothing coefficients are fixed to α i = 0.1 and β i = 0.5 (values chosen experimentally) and the interpolation coefficients are optimised by testing all combinations of values between 0 and 1 in 0.1 steps (11 distinct values). Then, the smoothing coefficient pairs are optimised separately and sequentially in the same range. Finally, the interpolation coefficients are fine-tuned around the original optimum (±20%, 11 values) using the newly trained smoothing coefficients.

Sub-models

As a proof of concept, we have developed three sub-models that model the three most important aspects of chords: chord progressions, relation to the tonality and relation to the melody:

P 1 = P(C t |C t-1 ), (11)
P 2 = P(C t |T t ), (12)
P 3 = P(C t |M t ). (13)
To train these three sub-models, we have used a collection of around 2000 lead sheets encoded in the MusicXML format that are freely available on the Wikifonia web page [START_REF]Wikifonia Foundation[END_REF]. We converted each lead sheet into three sequences of symbols representing tonality, melody and chords. In order to do so, we first partitioned the input into regular time frames of length 1/3, 1/2, 1, 2, 4, 8 or 16 beats. For every frame, we defined the melody variable M t as the unordered list of pitch classes of the 12-tone chromatic scale appearing in the melody, which can take 2 12 = 4096 distinct values. This is similar to the melody encoding in [START_REF] Chuan | A comparison of statistical and rule-based models for stylespecific harmonization[END_REF] and [START_REF] Raphael | Functional harmonic analysis using probabilistic models[END_REF]. The tonality T t was encoded as one of 24 different key labels resulting from the combination of 12 tonics (C, C♯, D, D♯, E, F, F♯, G, G♯, A, A♯, B) and 2 modes (major or minor). The chord C t was labelled by one of 13 root pitch classes (C, C♯, D, D♯, E, F, F♯, G, G♯, A, A♯, B or "none" for non-chords) and one of 27 chord types (major, minor, dominant, diminished, half-diminished, augmented, power, suspended-second, suspended-fourth, major-sixth, minor-sixth, major-seventh, minor-seventh, dominant-seventh, diminished-seventh, augmented-seventh, major-ninth, minor-ninth, dominant-ninth, augmented-ninth, minor-eleventh, dominant-eleventh, major-minor, minor-major, major-thirteenth, dominant-thirteenth or "none" for non-chords), resulting in N = 351 distinct chord labels in total. The chord C 0 before the beginning of the song was assumed by convention to be "none". In the case of a time frame containing more than one key or chord, the longest lasting key and chord labels within that frame were selected.

The distribution of tonalities in the training set is shown in Fig. 3. The C-major key appears to be dominant, which will have an impact on the design of the models in order not to bias them toward that particular tonality, as explained in Subsection 3.2.

Chord prior

The chord prior P(C t ) is used for smoothing and as a reference model in the evaluation. The prior chord distribution trained on the training dataset is presented in Fig. 4. Note that the major, dominant, minor, minor-and majorseventh chords make up for the vast majority of the chords in the dataset. Also, due to the dominance of the C-major key in the training set, the pitch classes C, F and G have visibly higher probabilities than the other root pitch classes.

Chord bigram model

The chord progression model is built under the Markov assumption, resulting in a bigram model P 1 (C t |C t-1 ). Although longer chordal dependencies can perform better, this has already been studied by others [START_REF] Paiement | Probabilistic melodic harmonization[END_REF][START_REF] Scholz | Robust modeling of musical chord sequences using probabilistic n-grams[END_REF] and is not the focus of this paper. In order to avoid problems with data sparsity, the model was trained with state tying: probabilities of all relative chord transitions were tied together, so for example the probability of transition from C-major to G-minor (7 semitones) is identical to that of transition from G-major to D-minor (also 7 semitones). This is motivated by the observation that in tonal music songs can be freely transposed between all keys without any loss of musical correctness [START_REF] Papadopoulos | Large-scale study of chord estimation algorithms based on chroma representation and hmm[END_REF], and by the dominance of the C-major key in the training set (see Fig. 3), which would otherwise result in a biased chord distribution towards the common chords of that key.

The resulting conditional chord distribution can be observed in Fig. 5. The distribution for a 1-beat analysis frame (top of Fig. 5) is very concentrated towards the previous chord (G-major to G-major transition), because chords typically last for at least few beats. On the other hand, using a 16-beat analysis frame makes the bigram probabilities being more evenly distributed, though still dominated by the transition to the same chord.

C C# D D# E F F# G G# A A# B none 0 0.

Tonality model

In the tonality model P 2 (C t |T t ), for the same reasons as explained in Subsection 3.2, the chords corresponding to the same scale degree in different keys were tied together. In other words, observing, e.g., a dominant major chord in one key increases the probability of dominant major chords in all keys. The resulting tonality model distribution is depicted in Fig. 6. Notice that knowing the current tonality to be, in this case, T t = C-major increases the dominance of the most common degrees: the dominant (V), subdominant (IV), supertonic (ii), but mostly the tonic (I).

Melody model

Finally, for the same reasons again, state tying was used for the melody model P 3 (C t |M t ) as well. Note patterns with the same content relative to the chord root were given identical probabilities, e.g., the unordered note combination (C,G) in the chord of C-major is equally probable as the note combination (D♯,A♯) in the chord of D♯-major. The resulting melody model distribution is shown in Fig. 7. Note that having more melodic information, i.e., more notes in a frame (here C, E and G in one frame in the bottom part of the plot), makes the chord distribution significantly sparser. This means that a system with longer frames will have a more informative melody model, because a single

RR n°8110

Root frame will contain more melody notes.

Interpolation coefficients

The values of the generatively trained interpolation coefficients a i and b i are presented in Fig. 8. The coefficients for linear and log-linear interpolation follow a similar trend: the bigram model is given a progressively lower weight as the frame length increases, while the melody model behaves in an opposite manner. Indeed, for large frames the bigram model becomes less informative due to the lack of need of modelling the chord duration, while the melody model becomes more informative due to the larger melodic context.

Evaluation

Due to a large range of goals in the existing literature, we can observe a variety of ways of evaluating harmonisation algorithms: through theoretic evaluation of the modelling power by means of cross-entropies [START_REF] Allan | Harmonising chorales by probabilistic inference[END_REF][START_REF] Paiement | Probabilistic melodic harmonization[END_REF], through comparison of generated chord sequences with ground-truth chord annotations [START_REF] Chuan | A comparison of statistical and rule-based models for stylespecific harmonization[END_REF][START_REF] Chuan | A hybrid system for automatic generation of style-specific accompaniment[END_REF], through comparison of single predicted chords with ground-truth [START_REF] Cunha | An intelligent hybrid model for chord prediction[END_REF], or through subjective listening tests [START_REF] Simon | MySong: automatic accompaniment generation for vocal melodies[END_REF][START_REF] Lee | i-Ring: A system for humming transcription and chord generation[END_REF]. Because of the novelty of the proposed solutions, as well as that of the field of automatic harmonisation itself, many papers did not offer any evaluation [START_REF] Phon-Amnuaisuk | The four-part harmonisation problem: a comparison between genetic algorithms and a rule-based system[END_REF][START_REF] Ebcioglu | An expert system for chorale harmonization[END_REF][START_REF] Pachet | Musical harmonization with constraints: A survey[END_REF][START_REF] Gang | Tuning a neural network for harmonizing melodies in real-time[END_REF]. In this paper, we have chosen to perform two complementary evaluations: we first use the theoretic cross-entropy-based evaluation as in [START_REF] Allan | Harmonising chorales by probabilistic inference[END_REF][START_REF] Paiement | Probabilistic melodic harmonization[END_REF] as a convenient way to validate our interpolation-based approach and we then perform an objective evaluation of the generated chord sequences in the same manner as in [START_REF] Chuan | A comparison of statistical and rule-based models for stylespecific harmonization[END_REF]. The code of our algorithm is available online at http://versamus.inria.fr/software-and-data/harmonization. tbz2.

Out of the 2000 Wikifonia lead sheets, 100 lead sheets were used as a test set, 100 as a validation set, and the rest were used as a training set.

Cross-entropy

An efficient way of determining the modelling power of a model is to compute the normalized negative log-likelihood, or cross-entropy [START_REF] Jurafsky | Speech and Language Processing[END_REF]. For base-2 logarithms, it can be interpreted as the average number of bits (b) required to encode a single chord symbol (Shannon's optimal code length). So, naturally, the smaller the value, the higher the prediction power. The cross-entropy is calculated from the test chord, melody and tonality sequences C, M and T as

H(C) = - 1 T log 2 P(C|M, T, Λ) = - 1 T T t=1 log 2 P(C t |C 1:t-1 , X 1:t ) (14)
where Λ denotes the model parameters and T is the number of frames in the test set. The cross-entropy is upper bounded by the cross-entropy of the non- informative uniform chord prior, which in our case is equal to

H U = -log 2 1 N = 8.46 b, (15)
where N = 351 is the number of distinct chord labels. The cross-entropies obtained for different frame lengths (from 1/3 to 16 beats) and different generatively trained log-linear combinations of the bigram (B), tonality (T) and melody (M) models are plotted in the upper part of Fig. 9. We can observe that the predicting power of the melody (M) and tonality (T) models, although better than the prior alone, is poor for small frame lengths (about 5 bits/frame), but improves by about 0.5 bits/frame as the frame length increases. This is logical, as both the tonality and the chord are musical quantities that depend on a much wider context than a single melody note. However, using large frame lengths limits the temporal precision of the chord estimation, since chords can typically change on any beat (though typically on the down beat). We also observe the benefit of model interpolation: the combined melody and tonality (M+T) model is better than either of the two models alone, and the combined melody and tonality and bigram (M+T+B) model is better the each of these three models alone and than the M+T combination. For 2-beat frames, the latter improvement is equal to 0.37 bit/frame (11%). Still in the upper part of Fig. 9, we can see that the cross-entropy decreases monotonically with decreasing frame length for those models that include the bigram chord progression model, namely B and M+T+B. In fact, it would decrease asymptotically to zero for infinitesimal frame lengths, because predicting the next chord given the current one would be getting easier: one would simply have to predict the same chord and that prediction would be increasingly correct. Therefore we found it useful to complement this plot with a plot of cross-entropies normalised per beat instead of per frame in the bottom part of Fig. 9. Per-beat cross-entropies are proportional to the total amount of information in the entire dataset, given the model (because the number of beats in the dataset is fixed, as opposed to the number of frames). From that plot we can see that we actually get more informative bigram models for larger frame lengths. On the other hand, this evaluation measure is biased towards larger frame lengths since it is bounded by the per-beat cross-entropy of the uniform chord distribution H U /B = (log 2 N )/B, where B is the number of beats per frame, which decreases asymptotically towards zero as B increases.

In the light of the above discussion we conclude that the best frame length for harmonisation is in the range of few beats, where the cross-entropy is between 2 bits/frame for 1-beat frames and 4 bits/frame for 8-beat frames. Fig. 10 shows the reduction of cross-entropy achieved by log-linear interpolation with respect to linear interpolation (H lin -H loglin ). Although log-linear interpolation is more time consuming (due to the need for re-normalisation by Z), it offers significantly higher modelling power when several models are combined: up to 0.32 bit/frame (10%) for the M+T model and up to 0.26 bit/frame (6.5%) for the M+T+B model.

Harmonisation

In a second experiment we compare the chord sequence generated by the full discriminatively trained model M+T+B with the ground truth chord labels in the test lead sheet files. In this experiment, the timing of the ground truth chord sequence is preserved and does not depend on the chosen frame length. Chords were estimated via a Viterbi-like algorithm, which finds the most likely sequence of chords given the melody and tonality. The tonality was assumed to be known, because often several tonal interpretations are possible [START_REF] Aiello | Musical perceptions[END_REF] and we want the resulting chord sequences to be comparable to the ground truth. Similarly to the MIREX competition [START_REF] Downie | The music information retrieval evaluation exchange: Some observations and insights[END_REF][START_REF]Music Information Retrieval Evaluation eXchange[END_REF], the estimated chord types were subsequently clustered into a smaller number of triads (major, minor, augmented, diminished, suspended second or suspended fourth). We compare the chords either in terms of their root note only or in terms of their root note and their triad chord type using two alternative root accuracy measures: a binary one and a weighted one. The latter uses the following weights: 1 for correct root pitch class estimation, 0.5 for 5 or 7 semitone errors (perfect fourth or fifth), 0.3 for 3, 4, 8 or 9 semitone errors (minor or major third and minor or major sixth), 0.1 for 2 or 10 semitone errors (major second or minor seventh) and 0 for other errors.

For comparison, we have used the results of the state-of-the-art rule-based lead sheet-like harmonisation system of Temperley and Sleator, which is freely available on-line at [START_REF] Temperley | Harmonic Analyzer[END_REF]. The lead sheets were converted to the input format of that algorithm with a time precision ("BaseUnit") of a dotted sixty-fourth note and the metrical structure was generated based on the time signatures and the upbeat durations extracted from the MusicXML files. Note that this algorithm estimates only the chord roots, not the chord types.

The results for different frame lengths are plotted in Figs 11,12 and 13. For all evaluation metrics, the log-linearly interpolated models offer better accuracy than the linear ones and the best results are most often achieved for a frame length of 2 beats. For shorter frame lengths the melody provides less information, while for longer frame lengths the temporal resolution of the generated chord sequence becomes too coarse. For that frame length, the proposed algorithm with log-linear interpolation outperforms the reference algorithm by 5.5% absolute (17% relative) in terms of root note accuracy and by 4.1% absolute in terms of weighted root note accuracy.

Conclusion

In this paper we have presented a novel method of building versatile statistical models of chords for harmonisation by joining multiple simpler sub-models by means of linear or log-linear interpolation. To test this idea, we have trained and combined in this way three sub-models: the tonality, the melody and the chord bigram model. We have evaluated the resulting interpolated models in terms of their cross-entropy and observed that log-linear interpolation yields a model whose cross-entropy is lower than the best of the component models and also better than tha achieved by linear interpolation. We have then performed a series of harmonisation experiments, where we have observed that the proposed log-linearly interpolated model offers higher root chord accuracy than the reference rule-based harmoniser from [START_REF] Temperley | Modeling meter and harmony: A preferencerule approach[END_REF] by up to 5% absolute.

In future work, a larger number of more complex sub-models could be investigated for further improvement in terms of chord accuracy. Subjective listening tests could also be used to analyse the quality of the harmonisations in more details. Finally, the model interpolation methodology could be applied to other music information retrieval tasks that would potentially benefit from modelling several musical aspects simultaneously.

Figure 1 :

 1 Figure 1: Example harmonisation of the lead melody from Radiohead's Karma Police.

Figure 2 :

 2 Figure 2: A typical HMM for melody harmonisation (top) compared to the proposed model (bottom).

Figure 3 :

 3 Figure 3: Histogram of the key labels in the training dataset.

Figure 4 :

 4 Figure 4: Distribution of chords P(C t ) for 1-beat frames. On top: distribution of chord types; on the right: distribution of chord roots. Chord types are sorted by their occurrence frequency.

Figure 6 :

 6 Figure 6: Conditional distribution of chords P 1 (C t |T t ) for 1-beat frames, plotted for the tonality of T t = C-maj. The plot is accompanied with linear-scale bar plots of chord root and chord type distributions on the sides.

Figure 7 :

 7 Figure 7: Conditional distribution of chords P 3 (C t |M t ) for 1-beat frames, plot-ted for a single C note (top) and for a time frame containing three notes: C, E and G (bottom). Each plot is accompanied with linear-scale bar plots of chord root and chord type distributions on the sides.

Figure 8 :

 8 Figure 8: Values of the generatively trained interpolation coefficients for linear (top) and log-linear (bottom) interpolation.

Figure 10 :

 10 Figure 10: Cross-entropy reduction of log-linear over linear interpolation.

Figure 11 :

 11 Figure 11: Root note accuracy obtained for different model frame lengths.

Figure 12 :Figure 13 :

 1213 Figure 12: Weighted root note accuracy obtained for different model frame lengths.

  Figure 5: Conditional distribution of chords P 1 (C t |C t-1 ) (chord transition prob-ability) for 1-beat frames (top) and 16-beat frames (bottom), plotted for the previous chord C t-1 = G-major. The plot is accompanied with linear-scale bar plots of chord root and chord type distributions on the sides.

	Root pitch class
	major	dominant	minor	minor-seventh	major-seventh	major-sixth	dominant-ninth	dominant-seventh	none	minor-sixth	half-diminished	suspended-fourth	diminished-seventh	diminished	augmented-seventh	minor-ninth	dominant-13th	augmented	major-ninth	suspended-second	dominant-11th	minor-11th	minor-major	major-minor	augmented-ninth	power	major-13th
									1	0.2		0.3	0.4		0.5	0.6		0.7		0.8	0.9		1

Average cross-entropy per frame for log-linear interpolation

		10
					M+T+B
					M
					T
		8			B
					M+T
	Average cross-entropy [bits/frame]	4 6			uniform chord distribution prior only
		2
		0
		1/3	1/2	1	2	4	8	16
					Frame length [beats]
			Average cross-entropy per beat for log-linear interpolation
		10
					M+T+B
					M
					T
		8			B
					M+T
					uniform chord distribution
	Average cross-entropy [bits/beat]	4 6			prior only
		2
		0
		1/3	1/2	1	2	4	8	16
					Frame length [beats]
	Figure 9: Cross-entropies calculated for the test dataset, normalised per time
	frame (top) and per beat (bottom). M stands for the melody model, T for the
	tonality model and B for the bigram chord progression model. Models were
	combined using log-linear interpolation.

Per-

frame entropy reduction of log-linear over linear interpolation

		0.6
			M+T+B
	Average cross-entropy [bits/frame]	0.2 0.4	M T B M+T
		0.0
		1/3	1/2	1	2	4	8	16
					Frame length [beats]

RR n°8110

Publisher Inria Domaine de Voluceau -Rocquencourt BP 105 -78153 Le Chesnay Cedex inria.fr ISSN 0249-6399

Acknowledgment

This work was supported by INRIA under the Associate Team Program VER-SAMUS (http://versamus.inria.fr/).