Fabien Durand

HD0L ω-EQUIVALENCE AND PERIODICITY PROBLEMS IN THE PRIMITIVE CASE (TO THE MEMORY OF G. RAUZY)

In this paper I would like to witness the mathematical inventiveness of G. Rauzy through a personnal exchanges I had with him. The objects that will emerge will be used to treat the decidability of the HD0L ω-equivalence and periodicity problems in the primitive case.

1. Introduction 1.1. Some words about Gérard Rauzy. I started my PhD in 1992 under the direction of G. Rauzy. The topic of my PhD was Substitutions and Bratteli diagrams. G. Rauzy was a worldwide recognized expert on substitutions (and many other topics) but Bratteli diagram representations of substitution subshifts was not one of these subjects. Thus I had an unformal second advisor, B. Host. Later, as a common agreement, he became my advisor. Some days after I begun my PhD he told me that he noticed something amazing about substitutions. He wrote on the blackboard the Fibonacci fixed point on the alphabet {0, 1} starting with 0. He considered all occurrences of the letter 0 and coded the words between two consecutive such occurrences. It is not surprising (see [Morse and Hedlund 1940]) that only two words appear : 01 and 0. If you respectively code them by 0 and 1 then you again obtain the same sequence. This is not surprising too because this sequence is the fixed point of the morphism defined by 0 → 01 and 1 → 0. Then he wrote a long prefix of the fixed point of the Morse sequence that we call x. Applying the same coding process (see below) to x it clearly appeared that the new sequence called the derived sequence of x, and denoted by D(x), is not the same for the obvious reason that 3 different letters occur. He repeated the process to this new sequence. Let's do it. The result was a sequence on a 4 letter alphabet. Once again he repeated the process. Amazingly we obtained the same sequence (the same long prefix). Repeating the process will clearly produce the same sequence. The process is periodic.

x

= 011 01 0 011 0 01 011 01 0 01 011 0 011 01 0 . . . D(x) = 0 1 2 0 2 1 0 1 2 1 0 2 0 1 2 . . .

Then he showed me another example where the period was 2. He was wondering if such a period always exists for fixed points of primitive substitutions. He knew that when this process is eventually periodic then the sequence is primitive substitutive and he explained me the proof on the blackboard. This proof was very useful for my PhD and for the present paper : this is the construction exposed in Proposition 17. It became a second subject in my PhD that I succeeded to solve in [Durand 1998a] and that provided me many arguments to obtain an extension of Cobham's theorem for primitive substitutions ( [Durand 1998b]).

When he showed me this process I found curious to look at it. Why was it interesting ? In fact, at this time, I did not know much about the work of G. Rauzy and almost nothing about interval exchange transformations (iets) and what he did with. I realized quickly that this process was in fact the induction process defined on dynamical systems by the Poincaré first return map : the derived sequence D(x) generates a subshift that is conjugate to the induced dynamical system (X, T ) on the cylinder set [0], where (X, T ) is the subshift generated by x. Thus, this corresponds to the induction process he defined [Rauzy 1979] to study iets. He induced on the left most interval and showed that the induced tranformation is again an iet with the same permutation. This is also the case for primitive substitutions : derived sequences of purely primitive substitutive sequences are purely primitive substitutive sequences ( [Durand 1998a]). His goal was to tackle the Keane's conjecture saying that almost all iets are uniquely ergodic. This was solved independently in [Masur 1982] and [Veech 1982]. Later M. Boshernitzan and C. R. Carroll [Boshernitzan and Carroll 1997] proved that, for an iet defined on a quadratic field, "every consistent method of induction is eventually periodic" (up to rescaling). This is what I obtained in [Durand 1998a] in the context of primitive substitutions and that was expected by G. Rauzy. Thus I imagine the intuition of G. Rauzy was certainly leaded by iet considerations. I am very grateful to Gérard Rauzy for the wonderful subject he gave me and for the opportunity he offered me to do research.

1.2. The HD0L ω-equivalence problem. In this paper we propose to apply what the author developped in [Durand 1998a] and [Durand 1998b] (to answer the question of G. Rauzy) to solve the following problems in the primitive case. To describe these problems we need some notations. Let σ : A * → A * and τ : B * → B * be morphisms, and, u and v be words such that lim n→∞ σ n (u) and lim n→∞ τ n (v) converge to respectively x ∈ A N and y ∈ B N . We also need two morphisms φ :

A * → C * and ψ : B * → D * . 1.2.1. The HD0L ω-equivalence problem.
This problem is open for more than 30 years and ask whether it is decidable to know if two morphic sequences are equal. More precisely, Is the equality "φ(x) = ψ(y)" decidable?

In this paper we answer positively to this problem for primitive morphisms. The equality problem "x = y" (also called the D0L ω-equivalence problem) was solved in 1984 by K. Culik II and T. Harju [START_REF] Culik | [END_REF]]. An other proof was given by J. Honkala in [Honkala 2007] (see also [Honkala 2009a]).

1.2.2. More on the equality of sequences generated by morphisms. Once the equality "φ(x) = ψ(y)" or "x = y" is satisfied, it is natural to look for relations between the involved morphisms.

In [START_REF] Séébold | On the conjugation of standard morphisms[END_REF]] it is shown than when two invertible substitutions on a two letter alphabet share the same fixed point, then these substitutions are powers of some other morphism (for different proofs see also [Berthé, Frettlöh and Sirvent:preprint] and [Rao and Wen 2010]). For more general results on this topic see [Krieger 2008, Diekert and[START_REF] Diekert | [END_REF]].

We will give a proof of a weaker result but that can be extended to bigger alphabets : if two primitive substitutions share the same fixed point and satisfy some conditions on their return words then they have a common power.

1.2.3. The decidability of the periodicity of HD0L infinite sequences. This is also an open problem for more than 30 years. It ask :

Is it decidable to know whether φ(x) is ultimately periodic ?

J.-J. Pansiot proved in [START_REF] Pansiot | Decidability of periodicity for infinite words[END_REF]] that the periodicity of x is decidable. We answer positively to this problem in the primitive case and in a forthcoming paper the author should complete this work to close the problem.

In the sequel the vocabulary of D0L and HD0L systems will be replace by the vocabulary of substitutions.

1.3. Some definitions.

1.3.1. Words and sequences. An alphabet A is a finite set of elements called letters.

A word over A is an element of the free monoid generated by A, denoted by A * . Let x = x 0 x 1 • • • x n-1 (with x i ∈ A, 0 ≤ i ≤ n -1) be a word, its length is n and is denoted by |x|. The empty word is denoted by ǫ, |ǫ| = 0. The set of nonempty words over A is denoted by A + . The elements of A N are called sequences.

If x = x 0 x 1 • • • is a sequence (with x i ∈ A, i ∈ N) and I = [k, l] an interval of N we set x I = x k x k+1 • • •
x l and we say that x I is a factor of x. If k = 0, we say that x I is a prefix of x. The set of factors of length n of x is written L n (x) and the set of factors of x, or the language of x, is denoted L(x). The occurrences in x of a word u are the integers i such that x [i,i+|u|-1] = u. If u has an occurrence in x, we also say that u appears in x. When x is a word, we use the same terminology with similar definitions. The sequence x is ultimately periodic if there exist a word u and a non-empty word v such that x = uv ω , where

v ω = vvv • • • . In this case |v| is called a length period of x. It is periodic if u is the empty word. A word u is recurrent in x if it appears in
x infinitely many times. A sequence x is uniformly recurrent if every factor u of x appears infinitely often in x and the greatest difference of two successive occurrences of u is bounded.

Morphisms and matrices.

Let A and B be two alphabets. Let σ be a morphism from A * to B * . When σ(A) = B, we say σ is a coding. Thus, codings are onto. If σ(A) is included in B + , it induces by concatenation a map from A N to B N . These two maps are also called σ. To the morphism σ is naturally associated the matrix M σ = (m i,j ) i∈B,j∈A where m i,j is the number of occurrences of i in the word σ(j).

Substitutions and substitutive sequences.

A substitution is a morphism σ : A * → A * . If there exists a letter a ∈ A and a word u ∈ A + such that σ(a) = au and moreover, if lim n→+∞ |σ n (a)| = +∞, then σ is said to be prolongable on a.

Since for all n ∈ N, σ n (a) is a prefix of σ n+1 (a) and because |σ n (a)| tends to infinity with n, the sequence (σ n (aaa • • • )) n≥0 converges (for the usual product topology on A N ) to a sequence denoted by σ ω (a). The morphism σ being continuous for the product topology, σ ω (a) is a fixed point of σ: σ(σ ω (a)) = σ ω (a). A sequence obtained in this way by iterating a prolongable substitution is said to be purely substitutive (w.r.t. σ). If x ∈ A N is purely substitutive and φ : A * → B * is a morphism then the sequence y = φ(x) is said to be a morphic sequence (w.r.t. σ).

When φ is a coding, we say y is substitutive (w.r.t. σ).

Whenever the matrix associated to σ is primitive (i.e., when it has a power with strictly positive coefficients) we say that σ is a primitive substitution. We say σ is erasing if there exists b ∈ A such that σ(b) is the empty word.

Substitutions and return words

The main results of the paper are consequences of results in [Durand 1998a] (see also [Durand 2000]) taking care to be explicit on some bounds (Section 2.1). Since we are dealing with decidability problems, it is important to pay attention that these bounds are algorithmically computable because they are involved in the decidability problems we are dealing with.

2.1. Useful results on return words of morphic sequences. Let x be a uniformly recurrent sequence on the alphabet A and u a non-empty prefix of x. We call return word to u every word

x [i,j-1] = x i x i+1 • • • x j-1
, where i and j are two successive occurrences of u in x. The reader can check that a word v is a return word to u of x if and only if vu belongs to L(x) = {x [i,j] ; 0 ≤ i ≤ j}, u is a prefix of vu and u has exactly two occurrences in vu. The set of return words to u is finite, because x is uniformly recurrent, and is denoted by R x,u . The sequence x can be written as a concatenation

x = m 0 m 1 m 2 • • • , m i ∈ R x,u , i ∈ N,
of return words to u, and this decomposition is unique. We enumerate the elements of R x,u in the order of their first appearance in (m n ) n∈N . This defines a bijective map Θ

x,u : R x,u → R x,u ⊂ A * where R x,u = {1, • • • , #(R x,u )}.
It defines a morphism and the set Θ x,u (R *

x,u ) consists of all concatenations of return words to u. We should point out that the order used to define this morphism will be of particular importance in the sequel. In fact, you can change the enumeration and this will change the morphism. What is important in the sequel is that once you fixed an algorithm to define Θ x,u then you should use the same for all x and u. We denote by D u (x) the unique sequence on the alphabet R x,u characterized by

Θ x,u (D u (x)) = x.
We call it the derived sequence of x on u. It is clearly uniformly recurrent. The following proposition points out the basic properties of return words.

Proposition 1 ( [Durand 1998a]). Let x be a uniformly recurrent sequence and u a non-empty prefix of x.

(1) The set R x,u is a code, i.e., Θ x,u : R *

x,u → Θ x,u (R * x,u ) is one-to-one.

(2) If u and v are two prefixes of x such that u is a prefix of v then each return word to v belongs to Θ x,u (R * x,u ), i.e., it is a concatenation of return words to u.

(3) Let v be a non-empty prefix of D u (x) and w

= Θ x,u (v)u. Then • w is a prefix x, • Θ x,u Θ Du(x),v = Θ x,w and • D v (D u (x)) = D w (x).
Corollary 2. Let x be a uniformly recurrent sequence and u a non-empty prefix of x. If u is a prefix of v and v is a prefix of x then there is a unique morphism

Θ satisfying Θ x,u Θ = Θ x,v .
The next lemma asserts that the derived sequences of purely primitive substitutive sequences are purely primitive substitutive sequences.

Proposition 3 ( [Durand 1998a]). Let x = σ ω (a) where σ : A * → A * is a primitive substitution, and, u be a non-empty prefix of x. The derived sequence D u (x) is the unique fixed point of the unique substitution σ u satisfying

Θ x,u σ u = σΘ x,u .
The substitution σ u is called the return substitution on u.

Remark that it is easy to compute σ u and Θ x,u . Take w the first return word to u having an occurrence in x : Θ x,u (1) = w. Then, compute σ(w). It is necessarily a (computable) concatenation of return words to u and a prefix of x :

σ(w) = Θ x,u (v 1 ) for some v 1 ∈ R * x,u . This defines σ u (1) = v 1 .
If no return words different from w appears then we stop here and x is periodic : x = www . . . . If not, we do the same with Θ x,u (2) and we continue until no new return word appears. This will stop in finite time because x is ultimately recurrent. Remark that this procedure also defines the map Θ x,u . We proved the following lemma.

Lemma 4. Let x = σ ω (a) where σ : A * → A * is a primitive substitution. Let u be a prefix of x. Then, σ u and Θ x,u are algorithmically computable.

Bounds will be given with the help of the next 3 lemmas. For any morphism σ : A * → B * , we set |σ| = max a∈A |σ|. In the sequel R σ will denote the maximal difference between two successive occurrences of a word of length 2 in any fixed point of the primitive substitution σ. In order to bound R σ we need the following lemma.

Lemma 5. [Horn and Johnson 1990] If M is a primitive n × n matrix then there exists k ≤ (n -1)n n such that M k has strictly positive entries. Lemma 6. If σ is a primitive substitution defined on an alphabet with d letters, then R σ is bounded by 2|σ| (d-1)d d .

Proof. The proof is left to the reader.

Lemma 7. For all primitive substitution defined on an alphabet with d letters, then for all n

|σ n | ≤ Q σ min a∈A |σ n (a)|.
where

Q σ = max max 0≤n≤(d-1)d d +1 |σ n | mina |σ n (a)| , |σ| .
Proof. The proof is left to the reader. Now we can give upper and lower bounds for the length of the return words and for the cardinality of the sets of return words.

Theorem 8 ( [Durand 1998a, Durand, Host andSkau 1999]). Suppose x = σ ω (a), where σ is a primitive substitution, is non-periodic. For all non-empty prefixes u of x,

(1) for all words v belonging to R x,u , 1 Kσ |u| ≤ |v| ≤ K σ |u|, and

(2) #(R x,u ) ≤ 4K 3 σ , where K σ = Q σ R σ |σ|.
Proof. The upper bound in (1) is easily deduced from [Durand 1998a]. The lower bound can also be deduced from [Durand 1998a] but it is easier to obtain from [Durand, Host and Skau 1999]. The last bound can be deduced from (1).

Note that the right inequality of Statement ( 1) is sufficient to imply the two other inequalities (see [Durand, Host and Skau 1999]). From this theorem and the definition of σ u we deduce the following corollary.

Corollary 9. For all primitive substitution σ and all u such that u is a prefix of σ(u) we have

|σ u | ≤ |σ|K 2
σ . The result we stated above were used to prove the two following results that will be involved in the decidability problems. The first will be important in the purely substitutive (or D0L) context and the second in the morphic (or HD0L) case.

Proposition 10 ( [Durand 1998a]). Let σ be a primitive substitution. The set of the return substitutions of σ is finite and bounded by 4K 3 σ |σ|K 2 σ +1 .

Proposition 11. Let x = σ ω (a), where σ is a primitive substitution, and φ : A * → B * be a coding. For all prefix u of x there exists a unique morphism λ u satisfying

φΘ x,u = Θ φ(x),φ(u) λ u . Moreover, (1) λ u (D u (x)) = D φ(u) (φ(x)), (2) #{λ u ; u prefix of x} ≤ (K σ + 1) K 2 σ .
Proposition 12. With the settings of Proposition 11, let u be a prefix of x. Then, Θ φ(x),φ(u) and λ u can be algorithmically defined.

HD0L ω-equivalence problem

3.1. Decidability results for purely substitutive sequences in a primitive context. In this section we propose a strategy to decide whether two purely substitutive sequences are equal when the underlying substitutions are primitive.

The following lemma is easy to establish (see [Durand, Host and Skau 1999]). In the sequel the constants K σ and K τ are those given by Theorem 8.

Lemma 13. Let x be a uniformly recurrent sequence such that there exists K satisfying : for all u ∈ L(x) and all w ∈ R x,u , |w| ≤ K|u|. Then, for all n, all words of length (K + 1)n of L(x) have an occurrence of all words of length n.

Theorem 14. Let σ : A * → A * and τ : B * → B * be primitive substitutions,

x = σ ω (a) ∈ A N and y = σ ω (b) ∈ B N . Let K = 1 + (4K σ ) 3 |σ|K 2 σ +1 (4K τ ) 3 |τ |K 2 τ +1 , u n = x[0, (K σ + 1) n ], v n = y[0, (K σ + 1
) n ] and T n = (σ un , τ vn ). Then, the following are equivalent :

(1) x = y;

(2) there exist n and m less than K such that Corollary 15. The D0L ω-equivalence problem is decidable in the primitive case.

u n = v n , u m = v m , Θ x,
Proof. This is a direct consequence of Theorem 14. Let us keep the notation of this theorem. The bounds K σ , K τ and K are algorithmically computable. So are the u n and v n , n ≤ K. We gave the algorithm to compute return substitutions after their definition, hence T n is algorithmically computable and the D0L ω-equivalence problem is decidable.

Decidability results for morphic sequences in the primitive context.

The strategy to tackle the decidability of the HD0L ω-equivalence problem in the primitive case is the same as for the D0L case.

Theorem 16. Let σ : A * → A * and τ : B * → B * be two primitive substitutions,

x = σ ω (a) ∈ A N and y = σ ω (b) ∈ B N . Let φ : A * → C * and ψ : B * → D * be two codings. Let K = 1 + (4K σ ) 3 |σ|K 2 σ +1 (K σ + 1) K 2 σ (4K τ ) 3 |τ |K 2 τ +1 (K τ + 1) K 2 τ , u n = φ(x[0, (K σ + 1) n ]), v n = ψ(y[0, (K σ + 1) n ]) and T n = (λ un , σ un , λ vn , τ vn ).
Then, the following are equivalent :

(1) x = y;

(2) there exist n and m less than K such that

u n = v n , u m = v m , Θ x,un = Θ y,un , Θ x,um = Θ y,um and T n = T m .
Proof. The proof proceed as the proof of Theorem 14.

The following proposition will allow us to reduce the HD0L ω-equivalence problem to substitutive sequences. This result can also be found in [Cobham 1968[START_REF] Pansiot | Hiérarchie et fermeture de certaines classes de tag-systèmes[END_REF][START_REF] Cassaigne | Quelques propriétés des mots substitutifs[END_REF], Durand 1998a, Honkala 2009b]. Below I present the proof G. Rauzy showed me.

Proposition 17. Let x be a morphic sequence w.r.t. a primitive substitution. Then, x is substitutive w.r.t. a primitive substitution.

Proof. Let y be a purely substitutive sequence generated by the primitive substitution σ : A * → A * and ρ : A * → B * be a morphism such that x = ρ(y). Suppose y = σ ω (a). From Lemma 5, there exists k ≤ (m -1) m , where m = #A, such that all images of σ k have an occurrence of all letters of A. Thus, the morphism φ = ρσ k does not have the empty word as an image of letter. Moreover, x = φ(y).

Let D = {(c, k); c ∈ C, 0 ≤ k ≤ |φ(c)| -1} and ψ : C * → D * the morphism defined by: ψ(c) = (c, 0) . . . (c, |φ(c)| -1).
There is an integer

n (≤ max a∈A |φ(a)|) such that |σ n (c)| ≥ |φ(c)| for all c in C.
Let τ be the morphism from D to D * defined by:

τ ((c, k)) = ψ(σ n (c) [k,k] ) if 0 ≤ k < |φ(c)| -1, and τ ((c, |φ(c)| -1)) = ψ(σ n (c) [|φ(c)|-1,|σ n (c)|-1] ) otherwise.
For all c in C we have

τ (ψ(c)) = τ ((c, 0) • • • (c, |φ(c)| -1)) = ψ(σ n (c) [0,0] ) • • • ψ(σ n (c) [|σ n (c)|-1,|σ n (c)|-1] ) = ψ(σ n (c)),
hence τ (ψ(y)) = ψ(σ n (y)) = ψ(y). Hence ψ(y) is a fixed point of the substitution τ which begins with some letter (e, 0) and τ ψ = ψσ n . Let χ be the coding from

D * to B * defined by χ((c, k)) = φ(c) [k,k] for all (c, k) in D. For all c in C we obtain χ(ψ(c)) = χ((c, 0) • • • (c, |φ(c)| -1)) = φ(c),
and then χ(ψ(y)) = φ(y). Moreover, the relation τ ψ = ψσ n implies that for all k ∈ N we have τ k ψ = ψσ kn . Hence if σ is primitive then τ is primitive too.

Corollary 18. The HD0L-ω-equivalence problem is decidable in the primitive case.

Proof. It is a consequence of Proposition 17 and Theorem 16.

More on substitutions sharing the same fixed point

A substitution on a finite alphabet is invertible if when extended to a free group, it is an automorphism. Let x be a sequence on the finite alphabet A and suppose all letters of A have an occurrence in x. The stabilizer of x, denoted by Stab(x), is the set of morphisms σ : A * → A * such that σ(x) = x. We say that x is rigid whenever the monoid Stab(x) is generated by a single element (see [Krieger 2008] for more informations on this topic). It has been shown in [START_REF] Séébold | On the conjugation of standard morphisms[END_REF]] that all invertible substitutions generate rigid sequences (Theorem 19). There was in fact a gap in the proof. In [Rao and Wen 2010] the authors gave a geometrical proof (see also [Berthé, Frettlöh and Sirvent:preprint] for a different proof). The gap in the proof in [START_REF] Séébold | On the conjugation of standard morphisms[END_REF]] was recently corrected in [Richomme and Séébold:preprint]. The statement of this result is the following.

Theorem 19. Let σ and τ be primitive invertible substitutions defined on the alphabet {a, b}. They have powers sharing the same fixed point if and only if there exist a substitution ζ and positive integers i and j such that σ = ζ i and ζ j = τ .

As pointed out in [Krieger 2008], the invertibility assumption is crucial in this theorem. The following example due to Séébold (see [Krieger 2008]) gives a counterexample :

σ : a → ab, b → baabba, τ : a → abbaab, b → ba.
Indeed, σ ω (a) = τ ω (a) but σ and τ are not powers of the same substitution. For a 3-letter alphabet, it has been shown in [Wen and Zhang 1999] (see also [Richomme 2003]) that the monoid of invertible substitutions is not finitely generated. This leaded D. Krieger [Krieger 2008] to find a non rigid purely substitutive sequence x w.r.t. an invertible substitution. Below we propose a weaker version of Theorem 19 (Corollary 22) that applies to any kind of alphabets. Before we need a lemma, a proposition and the following notation : for any word w, V (w) is the vector whose entry of index a is the number of a in w. It is sometimes called the Parikh vector of w or the abelianized of w.

Lemma 20. Let x = σ ω (a) and y = τ ω (a) where σ and τ are two primitive substitutions (on an alphabet with d letters). Suppose x and y have a common prefix u satisfying :

(1) σ(w) = τ (w) for all w ∈ R x,u ;

(2) {V (w); w ∈ R x,u } generates Z d .

Then τ = σ.

Proof. The proof is left to the reader.

Proposition 21. [Durand 1998b] Let τ and σ be two primitive substitutions. Then, they share the same non-periodic fixed point σ ω (a) = τ ω (a) if and only if σ ω (a) and τ ω (a) share a common prefix u such that there exist two integers i and j verifying

τ i u = σ j u . (4.1)
Remark that Equality (4.1) means that τ i and σ j coincide on the set R x,u . Thus, using Lemma 20 we obtain the announced corollary.

Corollary 22. Let τ and σ be two primitive substitutions (on an alphabet with d letters) sharing the same non-periodic fixed point x = σ ω (a) = τ ω (a). Let i, j be the positive integers and u the prefix of x given in Proposition 21. Suppose the set {V (w); w ∈ R x,u } generates Z d . Then, there exist two integers i and j such that

τ i = σ j .
Applied to invertible substitutions on a 2-letter alphabet this corollary implies that when two invertible primitive substitutions share the same fixed point then they have non trivial common powers because it can be checked that for all u the set {V (w); w ∈ R x,u } generates Z 2 . For the counterexample of Seebold recalled above this means that all u satisfying Equality (4.1) is such that all vectors in {V (w); w ∈ R x,u } are Q-colinear.

Proposition 21 provides a way to study rigity problems : with a weaker assumption than the one in Corollary 22 but with specific families of substitutions it is clear that the conclusion of this corollary would also hold. To finish with this section we would like to mention the following complete substitutive extension of Cobham's theorem obtained in [Durand 2011]. This gives a strong restriction (without the primitivity assumption) to generate the same morphic sequence.

Theorem 23. Let α and β be two multiplicatively independent Perron numbers. Let A be a finite alphabet and x be a sequence of A N . Then, x is both α-substitutive and β-substitutive if and only if x is ultimately periodic, where "x is α-substitutive" means that x is substitutive w.r.t. a substitution whose dominant eigenvalue of its incidence matrix is α. This result also holds true for morphic sequences whenever the morphisms (instead of the codings) are non-erasing. For purely substitutive substitutions it (the necessary condition) was known from [Durand 2002].

Periodicity of HD0L infinite sequences

In the sequel σ : A * → A * is a primitive morphism, φ : A * → B * is a morphism, y = σ ω (a) is a sequence of A N and x = φ(y). We have to find an algorithm telling if x is ultimately periodic or not. Remark that as the morphism is primitive, the sequence x is uniformly recurrent, and, ultimate periodicity of x is equivalent to its periodicity. The following lemma is easy to prove. Lemma 24. For any periodic sequence z with length period L, if u is a word of L(z) whose length is greater than L then #R u = 1.

We recall (Proposition 1) that for all n there exists a prefix u n of y such that D n (x) = λ n (D un (y)) and u 1 = x 0 where λ n is the unique morphism from R * un to R * φ(un) defined by φΘ x,un = Θ φ(x),φ(un) λ n . We set σ n = σ un .

Proposition 25. Suppose φ is a coding and the incidence matrix of σ has strictly positive entries. Let K = 4K 3 σ |σ|K 2 σ +1 (K σ + 1) K 2 σ . Then, x is periodic if and only if there exists i less than K such that #R φ(ui) = 1.

Proof. If #R φ(ui) = 1, then it is clear that x is periodic. Now suppose x is periodic. From the choice of K and the pigeon hole principle, there exist i and j, i < j, such that (λ i , σ i ) = (σ j , λ j ). Thus, R ui = R uj and R φ(ui) = R φ(uj ) . Moreover, remark that we have D i (x) = D j (x). Thus, D i (x) = D i+k(j-i) (x) for all k ≥ 1. From Proposition 1, lim k→+∞ |u i+k(j-i) | = +∞. We conclude with Lemma 24.

Theorem 26. The HD0L periodicity problem for infinite sequences is decidable in the context of primitive morphisms.

Proof. First, let us make some reductions of the problem. From Lemma 5, we can suppose the incidence matrix of σ has strictly positive entries. Suppose φ is erasing. It cannot erase all the letters. Hence ψ = φσ is a non erasing morphism and remark that x = ψ(y). Consequently, we can suppose φ non-erasing. To end with the reductions, from Proposition 17 we can suppose φ is a coding. Now, with the notations introduced in this section, it suffices to compute R φ(un) from n = 1 to n = K. If for some i less than K, #R φ(un) = 1 then x is periodic, otherwise x is not periodic. It remains to explain how to compute R φ(un) . We recall the maps λ un : R y,un → R x,φ(un) are algorithmically definable (see Proposition 12). Hence it can be easily check that #R φ(un) = 1.

  un = Θ y,vn , Θ x,um = Θ y,vm and T n = T m .Proof. The necessary condition is obtained applying Proposition 10 and the pigeon hole principle. Let us show it is sufficient : suppose there exist n, m satisfying n < m ≤ K such thatu n = v n , u m = v m , Θ x,un = Θ y,vn , Θ x,um = Θ y,vm and T n = T m .From Proposition 3 we deduce that D un (x) = D um (x) and D un (y) = D um (y). Moreover, from Corollary 2 there exists a unique morphism Θ :R * → R * such that Θ x,un Θ = Θ x,um where R = R x,un = R x,um = R y,un = R y,um . Then, Θ x,un Θ(D um (x)) = Θ x,um (D um (x)) = x and Θ x,un Θ(D um (y)) = Θ x,um (D um (y)) = y.Thus, from the definition of D un (x) and D un (y) we deduce :Θ(D um (x)) = D un (x) = D um (x) Θ(D um (y)) = D un (y) = D um (y)and D un (x) and D un (y) are fixed points of Θ. But remark that from Lemma 13 and the choice of the u n , in all images of Θ there are all the letters of R and all images of Θ start with the same letter. Hence Θ is a primitive substitution and Θ has a unique fixed point. Consequently, D un (x) = D un (y) and x = Θ x,un (D un (x)) = Θ n (D un (y)) = y.


come

Mixed effect analyses are increasingly employed for the analysis of longitudinal efficacy or safety categorical data measured in clinical trials [START_REF] Ito | Exposure-response analysis for spontaneously reported dizziness in pregabalin-treated patient with generalized anxiety disorder[END_REF][START_REF] Mandemaand | Population pharmacodynamic model for ketorolac analgesia[END_REF][START_REF] Zingmark | Population pharmacokinetics of clomethiazole and its effect on the natural course of sedation in acute stroke patients[END_REF][START_REF] Zingmark | Modelling a spontaneously reported side effect by use of a Markov mixed-effects model[END_REF][START_REF] Sheiner | A new approach to the analysis of analgesic drug trials, illustrated with bromfenac data[END_REF]. For this purpose, proportionalodds model are frequently used. Other models, such as differential-odds model has also been proposed [START_REF] Kjellsson | Comparison of proportional and differential odds models for mixed-effects analysis of categorical data[END_REF]. Results from these analyses, i.e. models along with parameter estimates, are often further utilized for simulation of novel scenarios with respect to new dosing schedules or new patient populations. It is also advocated that these models should be used as an essential part of drug development programs. Therefore it is critical that these parameter estimates are unbiased and reliable.

We focus here on maximum likelihood estimation. However, as in all nonlinear mixed models, the integral of the likelihood function cannot be explicitly solved and various approximations are employed to approximate the true likelihood [START_REF] Verbeke | Mixed models for the analysis of categorical repeated measures[END_REF]. The most commonly used approximation is Laplace, available in the software NONMEM. Bias in parameter estimates for these models with NONMEM VI and SAS v.8 has been studied in detail and it has been shown that the use of Laplace approximation may result in severely biased estimates, especially when response categories are non-evenly distributed within the studied population [START_REF] Jonsson | Estimating bias in population parameters for some models for repeated measures ordinal data using NONMEM and NLMIXED[END_REF]. This is common when analyzing clinical data and therefore the importance of using exact evaluations of the likelihood integral is further accentuated. SAS implements a more exact evaluation of the likelihood using Adaptive Gaussian Quadrature . However this approach is time consuming, can be applied to models with limited number of random effects, and is not flexible for the analysis of pharmacokinetic/pharmacodynamic data after various repeated dosage regimen [START_REF] Fitzmaurice | Longitudinal Data Analysis[END_REF] In recent years, there have been several approaches/algorithms developed for the analysis of continuous data which are able to find maximum likelihood estimates without a need to compute the likelihood numerically [START_REF] Bauer | A survey of population analysis methods and software for complex pharmacokinetic and pharmacodynamic models with examples[END_REF][START_REF] Kuhn | Maximum likelihood estimation in nonlinear mixed effects models[END_REF] . A stochastic approximation version of EM algorithm linked to a Markov Chain Monte Carlo procedure has been suggested for maximum likelihood estimation within the non-linear mixed effects framework [START_REF] Kuhn | Maximum likelihood estimation in nonlinear mixed effects models[END_REF][START_REF] Lavielleand | Estimation of population pharmacokinetic parameters of saquinavir in HIV patients with the MONOLIX software[END_REF]. This procedure has been demonstrated to possess excellent statistical convergence properties as well as the ability to provide an estimator close to the maximum likelihood estimate in only a few iterations [START_REF] Lavielleand | Estimation of population pharmacokinetic parameters of saquinavir in HIV patients with the MONOLIX software[END_REF][START_REF] Savicand | Performance in population models for count data, part II: a new SAEM algorithm[END_REF]. In addition to the estimation of the maximum likelihood parameters, the SAEM algorithm also provides the user with the estimate of the Fisher Information Matrix, used to assess parameter estimate uncertainty. However, there have not been studies reported with respect to application of the stochastic algorithms to the analysis of the ordered categorical data

The aims of this study were (i) to extend the SAEM algorithm for estimation of parameters in categorical data mixed models, (ii) to evaluate its performance both for parameter and standard errors estimation via Monte Carlo simulation study and, (iii) to illustrate the performance of the algorithm on a real data example where both continuous PK and discrete PD data are simultaneously analyzed.

A proportional odds model with individual-specific random effect has been employed throughout the exercise to study properties of the new SAEM algorithm. This model has been used in the area of PKPD modeling [START_REF] Mandemaand | Population pharmacodynamic model for ketorolac analgesia[END_REF][START_REF] Sheiner | A new approach to the analysis of analgesic drug trials, illustrated with bromfenac data[END_REF][START_REF] Ezzetand | A random effects model for ordinal responses from a crossover trial[END_REF] and has been presented in detail elsewhere [START_REF] Mandemaand | Population pharmacodynamic model for ketorolac analgesia[END_REF][START_REF] Sheiner | A new approach to the analysis of analgesic drug trials, illustrated with bromfenac data[END_REF][START_REF] Jonsson | Estimating bias in population parameters for some models for repeated measures ordinal data using NONMEM and NLMIXED[END_REF][START_REF] Ezzetand | A random effects model for ordinal responses from a crossover trial[END_REF].

Methods:

The proportional odds model with random intercept

We assume that the response is an ordered categorical variable which takes its values in (0,1,…,M). Let y ij be the j th observation in the i th individual, i = 1, …, N. In the proportional odds model with random intercept, the cumulative probability of y ij being larger or equal to m (m=1, …, M), can be defined by the following logistic regression model logit 𝑃 𝑦 𝑖𝑗 ≥ 𝑚 = 𝛼 1 + … + 𝛼 𝑚 + (𝛽, 𝑥 𝑖𝑗 ) + 𝜂 𝑖 Equation 1

where logit denotes the logit function, 𝛼 1 + … + 𝛼 𝑚 specifies the baseline for category m (m=1, …, M); h is the function defining predictors or covariate effect, β is a vector of fixed effects which is the same across all categories, x ij is the predictor vector (e.g. time, dose, concentrations) for observation j of individual i and η i is the random effect of individual i. It is assumed that the random effects are normally distributed with mean 0 and variance ω 2 .

Implementation of the SAEM algorithm for categorical data models

The SAEM algorithm described in [START_REF] Verbeke | Mixed models for the analysis of categorical repeated measures[END_REF] for continuous data models has been extended to the ordered categorical data models in a similar manner as it has been done for the count data models [START_REF] Savicand | Performance in population models for count data, part II: a new SAEM algorithm[END_REF]. Let 𝜇 = 𝜶 𝟏 , 𝜶 𝟐 , … , 𝜶 𝑴 , 𝜷 𝟏 , 𝜷 𝟐 … 𝜷 𝑳 be the vector of fixed effects of the model and  be the variance-covariance matrix of the random effects  i (in our example, i is scalar and reduces to the variance ω 2 of  i ). Then, SAEM is an iterative procedure where at iteration k, a new set of random effects

 (k) =( i (k)
) is drawn with the conditional distribution 𝑝(𝜂| 𝑦 ; 𝜇 𝑘  , Ω 𝑘 ). Then, the new population parameters 𝜇 𝑘+1  , Ω 𝑘+1 are obtained by maximizing 𝑄 𝑘+1 𝜇, Ω defined as follows:

𝑄 𝑘+1 𝜇, Ω = 𝑄 𝑘 𝜇, Ω + 𝛾 𝑘 𝑙 𝑦, 𝜂 during the next 100 iterations.

An MCMC algorithm was used for the simulation step (see [START_REF] Verbeke | Mixed models for the analysis of categorical repeated measures[END_REF][START_REF] Jonsson | Estimating bias in population parameters for some models for repeated measures ordinal data using NONMEM and NLMIXED[END_REF] for more details).

Estimation of the Fisher information matrix

Let  ( ) be the set of population parameters to be estimated, and let 𝜽 be the The procedure consists in computing first 𝜃 with SAEM then applying the Louis formula with 𝜃 = 𝜃 which requires the computation of the conditional expectation and conditional variance defined in equation 4. These quantities are estimated by Monte-Carlo: 300 iterations of MCMC were performed for the numerical experiments. All extensions for SAEM algorithm described here have been implemented in software MONOLIX 3.1.

Simulation settings

The performance of the SAEM algorithm was evaluated via Monte Carlo simulation-. To Scenario D-E included a specific baseline, placebo and drug model through two additional  parameters (Eq. 6)

𝑙𝑜𝑔𝑖𝑡 𝑃(𝑦 𝑖𝑗 ≥ 𝑚) = 𝑎 1 + … + 𝑎 𝑚 + 𝛽 1 𝑐 𝑖𝑗 + 𝛽 2 𝑑 𝑖𝑗 + 𝜂 1 ; 1 ≤ 𝑚 ≤ 3 Equation 6
The placebo model was implemented as a step function ( 𝑐 𝑖𝑗 = 0 𝑖𝑓 𝑗 = 1 and Typical parameter values were chosen so as to mimic desired distribution of responses.

𝑐 𝑖𝑗 = 1 if 𝑗 = 2,3, 4
The studied variance range was 0.5 -40.

For each scenario, one hundred datasets each containing 1000 individuals were simulated with MATLAB. All estimation procedures were performed using MONOLIX 3.1.

Overview of studied scenarios is shown in Table I. For more details on the simulation design used, reader is kindly asked to refer to the original publication of Jönsson et al [START_REF] Jonsson | Estimating bias in population parameters for some models for repeated measures ordinal data using NONMEM and NLMIXED[END_REF]. For simplicity in the notations, all these formula are vectorial formula which holds for each component of Outcomes of all Monte-Carlo simulation studies exploring both, the parameter estimation procedure and estimation of Fisher information matrix, were presented as box-plots of relative estimation errors (REE) where bias and imprecision of the method, as defined by equation 8b and 8c, can easily be visualized.

Evaluation of the SAEM algorithm and the standard error estimates

CPU times needed for estimation of (i) population parameters, (ii) Empirical Bayes

Estimates (EBEs) which are individual random effects and (iii) standard error estimates, were also measured to assess the efficiency of the algorithm and the runtime for the analysis.

Illustration on a real data

The well-known real PKPD dataset of warfarin was used to evaluate novel SAEM algorithm and its ability to simultaneously analyse continuous and categorical data. The data were collected in 33 patients after a single dose of warfarin for 140h post dose. In total 251 pharmacokinetic (PK) observations and 232 pharmacodynamic (PD) observations (corresponds to inhibition of prothrombin complex synthesis -PCA (%))

were available [START_REF] O'reillyand | Studies on coumarin anticoagulant drugs. Initiation of warfarin therapy without a loading dose[END_REF][START_REF] O'reilly | Studies on the Coumarin Anticoagulant Drugs: The Pharmacodynamics of Warfarin in Man[END_REF]. Original PD variable was continuous variable expressed in percentages (0 -100 %); however for our purpose, we categorized the PCA variable into three ordered categories: 0 (if PCA is more than 50%), 1 (if PCA is between 33% and 50%) and 2 (if PCA was less than 33%). Of note, categorization of the continuous variable is done for illustration purpose only and it is not recommended to be done in Results:

Simulation study

Overall, the estimation procedure with the SAEM algorithm for mixed categorical data models, showed satisfactory performance with low bias and high precision.

Convergence was 100% for both parameter and standard error estimation. For parameter estimation, the absolute value of relative bias was less than 7.9% and 8.13%

for fixed effects and the random effect variances and RRMSE was less than 27% and 30% for fixed effects and the random effect variances over all tested scenarios. For standard error estimation, the absolute value of relative bias was less than 3.4% and 5.8% for fixed effects and random effect variances and RRMSE was less than 2.3% and 5.6% for fixed effects and random effect variances. The random effect variances, shown to be severely biased when estimated with Laplace method implemented in NONMEM (8) (8), were precisely estimated with SAEM, exhibiting relative bias ranging from 0.03% -8.13% across all studied scenarios. Detailed results for each scenario are listed below.

The distribution of REE for all scenarios and all parameter and standard error estimates are shown in Figure 2A-E

Illustration of a real data

With respect to the warfarin real data example, both parameter and standard error estimation was successful. Estimation procedure was completed in less than 2 minutes, for the model containing 8 typical parameters, 6 variances and 2 residual error

Discussion

The new SAEM algorithm has been developed, implemented and evaluated for application to categorical data models in the non-linear mixed effects framework. Five different scenarios using proportional odds model were evaluated, including those with non-even distribution of response categories. The algorithm was also implemented for computation of Fisher Information matrix in order to assess the uncertainty estimate.

The SAEM algorithm performed well under all tested model scenarios resulting in accurate and precise estimation of all parameters. Variances of scenarios with non-even distribution of response categories were accurately and precisely estimated, which was not reported previously in analysis with LAPLACE method (8). Shrinkage in EBE leads to linearization around zero for the random effects, close to a FO method, which is known to be biased [START_REF] Karlssonand | Diagnosing model diagnostics[END_REF][START_REF] Savicand | Importance of shrinkage in empirical bayes estimates for diagnostics: problems and solutions[END_REF]. Additionally, random effects enter models in a non-linear fashion; therefore these are most likely to suffer from the poor integral approximation, which was indeed observed in the previous work with severely biased variances [START_REF] Jonsson | Estimating bias in population parameters for some models for repeated measures ordinal data using NONMEM and NLMIXED[END_REF]. Gaussian quadrature method, as implemented in SAS, performed better than LAPLACE due to better numerical approximation of the likelihood integral , therefore EBE shrinkage influence was less pronounced compared to the LAPLACE approximation (8). Similar pattern was also observed when performance of these estimation methods was evaluated for count data [START_REF] Plan | Performance in population models for count data, part I: maximum likelihood approximations[END_REF]. Of note, it has been reported that this Gaussian Quadrture may become unstable and time consuming for more complex type of problems (7, 21).

The SAEM algorithm does not involve any likelihood estimation conditioned on EBEs

or any approximation of the model in computation of the likelihood integral and therefore does not suffer from any related biases. SAEM simulates large number of individual parameters using not only conditional modes, but also conditional variances at the current iteration. These conditional variances are large; therefore EBE shrinkage is not a problem under these circumstances.

Of note, small significant bias was observed in Scenario D for estimation of β 1 parameter, which is magnitude of treatment effect. This bias is most likely related to the small number of observations per subject as it disappeared when the number of observations per subject was increased.

The SAEM algorithm provides estimation of both the likelihood and Fisher information matrix, without linearization of the model. This is a favorable property of the algorithm, which leads to accurate and unbiased parameter and standard error estimates. The importance of unbiased standard error estimates has seldom been the topic of discussion. Standard errors are utilized in different aspects of pharmacometrics -they are an important aspect of prospective simulations, determination of the optimal study design, Wald test and exploration of competing study design scenarios .The SAEM algorithm appeared to satisfy requisite precision and unbiased estimates of parameter uncertainty.

In the previous analysis reported by Jönsson et al [START_REF] Jonsson | Estimating bias in population parameters for some models for repeated measures ordinal data using NONMEM and NLMIXED[END_REF] authors concluded that CPU time was not too burdensome and estimations were generally fast for methods investigated. The SAEM algorithm is easily applied for simultaneous modeling of continuous and discrete data and the most common application of this feature is in development of the PKPD models, with discrete PD variable. This case was also illustrated in our example with warfarin data. The advantages of simultaneous over sequential PKPD analysis has been demonstrated previously [START_REF] Zhang | Simultaneous vs. sequential analysis for population PK/PD data I: best-case performance[END_REF][START_REF] Zhang | Simultaneous vs. sequential analysis for population PK/PD data II: robustness of methods[END_REF], however to our knowledge such an analysis when PD variable is discrete has never been reported in the literature, even though simultaneous modeling of continuous and discrete data is possible with NONMEM VI.

The reason for that is that LAPLACE algorithm often becomes unstable whenever the model structure is more complex. The new SAEM algorithm as implemented in MONOLIX offer simple model coding and fast and stable estimation procedure.

The SAEM algorithm, which forms the core of MONOLIX is a freeware available at http://www.monolix.org and is based on thoroughly evaluated and documented thorough statistical theory. Monolix is an ongoing project implementing new statistical developments in a dynamic environment. The new version of MONOLIX program includes the extension of the algorithm for the analysis ordered categorical data as well as for count data [START_REF] Savicand | Performance in population models for count data, part II: a new SAEM algorithm[END_REF].

Conclusions

In conclusion, SAEM algorithm has been extended for the analysis of ordered categorical data. The parameters and standard errors are precisely and accurately estimated. The estimation procedure is stable and fast. Algorithm is easily extended for simultaneous modelling of continuous and discrete data.

  Implementation and Evaluation of the SAEM algorithm for longitudinal ordered categorical data with an illustration in pharmacokinetics-pharmacodynamics " Authors: Radojka M. Savic 1 , France Mentré 1 and Marc Lavielle 2 Institution: (1) UMR 738, INSERM -Université Paris Diderot, Paris, France (2) INRIA Saclay & Department of Mathematics, University Paris 11, Orsay , France email: Radojka.Savic@inserm.fr tel: +33 1 57 27 75 38 fax: +33 1 57 27 75 21 KEY WORDS: SAEM, CATEGORICAL DATA, MIXED MODELS, MONOLIX, PROPORTIONAL ODDS MODEL Introduction:

  allow a fair comparison with other algorithms, we used identical scenarios as presented previously in the paper of Jönsson et al where authors explored performance of Laplace and Adaptive Gaussian quadrature algorithms (8). Overall, five different scenarios (A-E) were used. In all scenarios response was a four level categorical variable that takes its values in {0,1, 2,3}. Scenarios A-C describe a baseline model 𝑙𝑜𝑔𝑖𝑡 𝑃(𝑦 𝑖𝑗 ≥ 𝑚) = 𝑎 1 + … + 𝑎 𝑚 + 𝜂 𝑖 ; 1 ≤ 𝑚 ≤ 3 Equation 5 with three different distributions of response categories: even (scenario A), moderately skewed (scenario B) and skewed (scenario C).

  ,), while the drug model was implemented as a linear function of the dose (𝑑 𝑖𝑗 = 0 𝑖𝑓 𝑗 = 1and 𝑑 𝑖𝑗 = 0,7.5,15,30 𝑖𝑓 𝑗 = 1,2,3,4 respectively. The distribution of response categories was even (scenario D) and skewed (scenario E).

For𝑹𝑩

  each scenario, the SAEM algorithm was used with the K=100 simulated datasets for computing the K parameter estimates, 𝜽 𝒌 , 𝒌 = 𝟏, … 𝑲 . The Fisher information matrix was also estimated for each data set, and its inverse was used to compute the K standard error estimates, 𝒔𝒆 𝒌 , 𝒌 = 𝟏, … 𝑲 . The empirical standard errors se* (i.e. the RMSE) were computed by equation 7:   stands for the true parameter value. To assess statistical properties of the proposed estimators, for each parameter, relative estimation errors 𝑹𝑬𝑬 𝜽 𝒌 , 𝒌 = 𝟏, … 𝑲 were computed as shown in equation 8a, where 𝒙 𝒌 = 𝜽 𝒌 . Similarly, for each estimated parameter standard error, relative estimation error 𝑹𝑬𝑬 𝒔𝒆 𝒌 , 𝒌 = 𝟏, … 𝑲 was computed, as shown in Equation 8a, where 𝒙 𝒌 = 𝒔𝒆 𝒌 . Each REE is expressed as a percentage (%). From the REEs, relative bias (RB), and relative root mean square errors (RRMSE) were computed for each parameter in each scenario as shown in Equations 8b-c. 𝑹𝑬𝑬 𝒙 𝒌 = 𝒙 𝒌 -𝒙 * 𝜽 * × 𝟏𝟎𝟎 Equation 8a 𝟏𝟎𝟎 𝒂𝒏𝒅 𝒙 = 𝜽 𝒐𝒓 𝒔𝒆 Equation 8c

the real analysis.

  The cut-offs chosen, are close to international normalized ration (INR) values commonly used in clinical practice to target optimal warfarin therapy. Low INR values (< 2) are associated with high risk of having a cloth (corresponding to category 0), high INR values (>3) with high risk of bleeding (corresponding to category 2), while targeted value of INR, corresponding to optimal therapy is in between 2 and 3 (corresponding to category 1).. The raw PKPD data are shown in Figure 1. The PK model fitted was one compartment model with first order absorption and a lag time. Effect compartment model was used to mimic an effect delay. Proportional odds model with random intercept was used to fit ordered categorical response. The drug model was a linear function of warfarin concentration.

Figure 3 .

 3 Figure 3. The output of MONOLIX run representing parameter estimates and respective

  The explanation for previously observed biases with LAPLACE was related to the poor approximation of the likelihood integral. Similarly to FOCE (first order conditional estimation approximation), LAPLACE approximation likelihood estimation computes and it involves linearization of the likelihood function by means of estimating EBEs at each iteration step. Whenever the estimated EBE distribution does not reflect the true random effect distribution, the method is expected to perform poorly. Reason for deviations of EBE distribution from the true random effect distribution is due to shrinkage phenomenon. Whenever data are sparse, which may be due to design, variability or the model itself, individual random effects will shrink toward zero and empirical variance of EBE will be smaller than the corresponding estimated Ω.

This was similarly observed

  with the new SAEM algorithm, with the median time for parameter estimation being less than half a minute. This is somewhat slower than reported times with LAPLACE (9.87 -17.1 s) and in the lower range of the reported times with Gaussian Quadrature (5.92 -165 s), for different GQ methods for scenario D and E). Of note, LAPLACE and GQ runs were performed on the computer with slightly faster processor (Pentium 2.8 GHz vs Intel 2.40 GHz for SAEM). All studied models converged successfully (100%), for both parameter estimation and standard error estimation with average CPU time being measured in seconds.

Figure legend:Figure 1 .

 1 Figure legend:

Figure

  Figure 2A-E. Distribution of relative estimation error (REE) for all parameters (left panel) and

Figure 3 .

 3 Figure 3. Implementation of the simultaneous analysis of continuous and discrete data in

Figure 4 .

 4 Figure 4. MONOLIX output for real data example. Parameter estimates are shown along with

Figure 5 .

 5 Figure 5. Probability change over time for each warfarin response category, based on the

  Figure 2A

Figure 4

 4 Figure 4

Figure 5 .

 5 Figure 5.

maximum likelihood estimate of  computed with SAEM. The Fisher Information matrix is defined as -

  𝝏 𝜽 𝟐 𝒍(𝒚; 𝜽 ) where 𝒍(𝒚; 𝜽 )

is the log-likelihood of the observations, computed with 𝜽 = 𝜽 . Several numerical experiments have shown that linearization of the model for estimating the Fisher information matrix (as implemented in MONOLIX 2.4) is satisfactory in case of continuous data (15). In this case, the linearization of the structural model allows transformation of the non- linear model into a Gaussian model, for which one the Fisher information matrix can be computed in a closed form. However

, this approach cannot be applied for discrete data models. As alternative we propose to compute a stochastic approximation of the Fisher Information matrix using the Louis formula (see

[START_REF] Kuhn | Maximum likelihood estimation in nonlinear mixed effects models[END_REF]

for more details): 𝜕 𝜃 2 𝑙 𝑦; 𝜃 = 𝐸 𝜕 𝜃 2 𝑙 𝛾, 𝜂; 𝜃 ;𝛾| 𝜃 + 𝑉𝑎𝑟(𝜕 𝜃 𝑙 𝛾, 𝜂; 𝜃 ;𝛾| 𝜃) Equation 4

for parameter, EBEs and standard error estimation are given in Table IV, for all studied scenarios.

  . The numerical results showing accuracy and precision for parameter estimation, measured as relative bias and relative root mean square error, are shown in Table II. The numerical results showing accuracy and precision for relative standard error estimation, measured as bias and root mean square error, are shown in

	The average CPU (Central Processing Unit) time per run over all scenarios was 29.6 s for
	parameter estimation, and 6.5 s for standard error estimation, with Matlab/C++
	implementation of the algorithm, when ran on laptop DELL D830 2.40GHz configuration.
	Median CPU times
	Table III indicating low bias (<5.78%) and high precision (RRMSE<7.42)

Table legend Table I .

 legendI Original study design and simulation settings. Distribution of response categories for originally simulated data sets and true parameter values used in simulations are presented.

Table II .

 II Relative bias and relative root mean square error (in %) for parameter estimates for all studied scenarios. These

results correspond to the visual ones shown in the left panel of Figures 2a-e.Table III . Relative bias and relative root mean square error (in %) for standard error estimates for all studied scenarios. These results correspond to the visual ones shown in the right panel of Figures 2a-e. Table IV. Median CPU time for parameter, EBE and standard error estimations for all studied scenarios.

 III

Acknowledgments

Radojka Savic was financially supported by a Postdoc grant from the Swedish Academy of Pharmaceutical Sciences (Apotekarsocieteten). We thank MONOLIX team, Hector Mesa and Kaelig Chatel for their help with implementation of the algorithm in the MONOLIX software. We also thank two anonymous reviewers for their valuable comments on the manuscript.


Vincent Breteau
email: vincent.breteau@enpc.fr

Fabien Leurent

Job decentralization and transportation use in a monocentric city *

Our objective in this paper is twofold: first, we want to give a theoretical founding to empirical findings of several works that emphasize the fact that while distance traveled increases with household location distance from the city center, transportation time tends to decrease, thus offering a strong incentive to sprawl. Second, we want to analyze the impact of job dispersal on city size, overall distance traveled and transportation cost, along with other urban variables, and spatial equity. We therefore develop an extended monocentric model of city taking into account employment dispersal and varying unit commuting costs. Using this model, we show that under specific conditions including employment dispersal and high marginal transportation cost around city center, the distance traveled by households from home to workplace increases with their distance from the city center, while private transportation costs they endure decrease. Then, based on a surplus analysis, we show that city size moderately increases with the level of employment dispersal, while overall home-to-work distance traveled decreases, suggesting that job decentralization might entail savings in social costs of transportation. However, our findings show that such dispersal could entail spatial inequity: the households living near the city center could suffer a welfare loss.

Introduction

Urban sprawl remains one of the most acute problems in urban areas and is one of the main stakes of urban economics today. Sprawl induces excessive space consumption and, concerning more specifically transportation, has consequences on energy consumption1 and congestion. It is well acknowledged in a vast literature, both theoretical and empirical, that urban sprawl and transportation are deeply linked (e.g. [START_REF] Fujita | Urban Economic Theory: Land Use and City Size[END_REF][START_REF] Kahn | La qualité de la vie et la productivité dans les villes étalées par opposition aux villes denses aux États-Unis[END_REF][START_REF] Glaeser | The greenness of cities[END_REF], even though their exact relationship in a given city remains doubtful. [START_REF] Melia | The paradox of intensification[END_REF] show indeed that urban intensification does not necessarily entail welfare improvement. It seems that beyond density, urban form is a largely determining factor of transportation use [START_REF] Aguiléra | urbaines et migrations alternantes -Les enseignements d'une des aires urbaines de Lille, Lyon et Marseille[END_REF][START_REF] Charron | La relation entre la forme urbaine et la distance de navettage : les apports du concept de « possibilité de navettage[END_REF]. As part of the equation, job decentralization is also often considered as a driver or a byproduct [START_REF] Bruegman | Sprawl: A Compact History[END_REF][START_REF] Buchfield | Causes of sprawl: A portrait from space[END_REF] of urban sprawl. More generally speaking, whereas the compact city is often seen as the green city, several works suggest that the cost of compacity is not insignificant (see for instance Castel, 2005, concerning France).

As already mentioned, theoretical studies have already tackled the link between urban sprawl and transportation. The seminal work of [START_REF] Wheaton | A comparative static analysis of urban spatial structure[END_REF] has emphasized, in a monocentric setting, the role of transportation cost decrease or income increase in urban sprawl. However, many other variables may play a role: preference for green environment, congestion, housing price, urban blight, etc. Among these, job localization, and more precisely job decentralization, plays an ambiguous role. Some authors consider it as one of the drivers of urban sprawl, pushing households at the urban periphery (e.g. [START_REF] Buchfield | Causes of sprawl: A portrait from space[END_REF]. But other authors emphasized the fact that jobs may follow the sprawl of households instead of drive it (e.g. [START_REF] Bruegman | Sprawl: A Compact History[END_REF]). Yet, few theoretical work have focused on this variable. The vast majority of theoretical urban models based on the standard monocentric model [START_REF] Alonso | Location and Land Use; Toward a General Theory of Land Rent[END_REF] assume indeed that the CBD size is zero or that transportation is costless in the CBD, which is actually equivalent. This is a quite reductive assumption, as job dispersion has been empirically shown (e.g. [START_REF] Genre-Grandpierre | Des « réseaux lents contre la dépendance automobile? concept et implications en milieu urbain. L'Espace Géographique[END_REF][START_REF] Mcmillen | Suburban subcenters and employment density in metropolitan chicago[END_REF][START_REF] Glaeser | Decentralized employment and the transformation of the American city[END_REF]. Some authors such as [START_REF] White | Firm suburbanisation and urban subcenters[END_REF][START_REF] White | Location choice and commuting behavior in cities with decentralized employment[END_REF][START_REF] White | Urban areas with decentralized employment: Theory and empirical work[END_REF] have explicitly introduced the locations of jobs in a monocentric setting. These works focused however on a long-term equilibrium where wages depend on the locations of the firms. While empirical works have shown the existence of a wage gradient (e.g. [START_REF] Timothy | Intra-urban wage variations, employment location, and commuting times[END_REF], [START_REF] Mills | Urban Economics[END_REF] have emphasized that the wage adjustment process entails the fact that the rent function form is unaffected by the localization of jobs, or differently put that job localization has no effect on urban sprawl, which is a doubtful conclusion.

In [START_REF] Breteau | Housing and commuting in an extended monocentric city[END_REF], henceforth BL10, we developed an extended monocentric model in which jobs are distributed exogenously in an employment area where the transportation cost is nonzero. In addition, we also assume households to be homogeneous in terms of utility and gross income and each household to have a fixed place of employment that influences its residential location. Our model is therefore one of medium-term equilibrium, where households are the only participants in the housing market and do not compete with firms. We have proved in BL10 the existence and uniqueness of the equilibrium, and derived the main comparative statics of the model.

In this paper, we use the model to provide answers to the following questions: does the decentralization of jobs constitute a potential driver of urban sprawl? Can job decentralization entail a transport-sober urban development? We thus give a founding to empirical findings showing that while distance traveled increases with household location distance from the city center, transportation time tends to decrease, offering a strong incentive to sprawl [START_REF] Fouchier | Les densités urbaines et le développement durable -Le cas de l'Île-de-France et des villes nouvelles[END_REF][START_REF] Genre-Grandpierre | Des « réseaux lents contre la dépendance automobile? concept et implications en milieu urbain. L'Espace Géographique[END_REF]. We also show that whereas job decentralization has a moderate negative impact on city size, the impact is positive and strong on overall distance traveled and transportation cost.

The paper is organized as follows: the next section describes the main assumptions of the model and recalls the main results of BL10. Section 3 presents our results concerning the first question addressed and section 4 discusses our results concerning the second question. Section 5 concludes.

The model and its assumptions

Firms and workplaces

We consider a city in which firms, and therefore jobs, are distributed around a centre in a disk of radius ρ f : we call it the employment area, composed of a CBD of radius r 0 exclusively reserved to jobs, and a mixed area ranging from r 0 to ρ f . Jobs are assumed to be distributed following a radial density f (ρ) that is non-zero on an interval of [0, ρ f ], yielding a cumulative distribution function F (ρ) = ρ 0 f (r)dr. The function F is therefore strictly increasing. The number of jobs is fixed at N , therefore F (ρ f ) = N .

Households

As regards households, we assume: (H1) that the households are located in a ring around the CBD. (H2) That the households are homogeneous in their preferences. (H3) That each household has only one working member employed by a firm in the employment area. (H4) That each household receives an income Y from this job, independent of residential and job location2 . Finally, (H5) each household has a fixed workplace. The choice of a residential location in r by a household whose workplace is ρ, leaves them a net income of I = Y -T (ρ, r), where T is the transportation cost. Households utility U depends on the size of its dwelling, s, and on the quantity z of a composite consumption good treated as numéraire. The function U is assumed to be increasing and continually differentiable into each of its variables. In general terms, each household is deemed to be a rational decision maker, which seeks to maximize its utility subject to its budget constraint.

Spatial structure of the residential area

We are interested here in the land market within the residential zone only: the unit price of land in r, or land rent, is denoted R(r). The opportunity cost of land, corresponding to an alternative use (e.g. agricultural), is denoted R A . The economic program of the household working at ρ is expressed as follows:

max r,z,s U (z, s) s.t. z + R(r)s ≤ Y -T (ρ, r).
(1)

The density of the households in r is denoted h(r). This is, with the land rent, the main endogenous variable in our model. For reasons of symmetry with the distribution of jobs, we introduce the cumulative household distribution function, H(r) = r r 0 h(r)dr. Conservation of the relative orders of jobs and households 3 allows us to define a function ρ → r ω (ρ) yielding, for a given workplace, the corresponding equilibrium location of the household, that is such that: H(r ω (ρ)) = F (ρ). Regarding the transportation cost, it is assumed that T (ρ, r) is a decreasing function in ρ and increasing in r. Finally, the amount of space available at r is denoted L(r).

Equilibrium characterizing system

On the demand side, we define the indirect utility function of a household as follows:

V (R, I ρ ) = max s,z U (s, z) z + Rs ≤ I ρ . ( 2
)
We then note W (r, ρ) = V (R(r), Y -T (ρ, r)), so that the economic program of the household given by ( 1) may be rewritten:

max r W (r, ρ). ( 3
)
On the supply side, as no land is left vacant and as land is allocated to the highest bider, we have, for all r ≤ r f :

h(r)s(r) = L(r), (4) R(r) ≥ R A .
(5)

Condition (3) yields:

∂W (r, ρ)/∂r = 0 at point ρ = ρ ω (r), (6)
that is:

R (r) ∂V ∂R = ∂T (ρ ω (r), r) ∂r ∂V ∂I . ( 7
)
Using Roy's identity, we obtain the Muth condition:

R (r) = - ∂T (ρ ω (r), r) ∂r /ŝ(R, I ρ,r ), (8)
with I ρ,r = Y -T (ρ ω (r), r), the net revenue, ρ ω (r) = F -1 •H(r) the reciprocal of r ω and ŝ the Marshallian demand for space. Condition (4) then yields:

H (r) = L(r)/ŝ(R, I ρ,r ), (9)
for the equilibrium lot size s(r) corresponds to ŝ(R, I ρ,r ), the bid-max lot size.

The urban equilibrium is thus characterized by the following differential system: R (r) = -∂T (ρω(r),r)

∂r /ŝ(R, I ρ,r ) H (r) = L(r)/ŝ(R, I ρ,r ) . ( 10
)
Using this system we were able to show existence and uniqueness of the equilibrium, defining for that purpose a terminal rent function4 , similar to the boundary rent function in [START_REF] Fujita | Urban Economic Theory: Land Use and City Size[END_REF]: see BL10 for details.

A specific model

In order to achieve analytical resolution when it is needed, we use specific forms for the exogenous variables of the model. Concerning the utility function, the literature has extensively considered the Cobb-Doublas form (or log-linear form, see [START_REF] Fujita | Urban Economic Theory: Land Use and City Size[END_REF]:

U (z, s) = U 0 z α s β , ( 11
)
with α, β > 0 and α + β = 1. We assume that residential land capacity and job distribution are uniform, that is L(r) = λ > 0 for all r > r 0 , f (ρ) = N/ρ f for ρ ∈ [0, ρ f ] and f (ρ) = 0 for ρ > ρ f , with N the number of households in the city. Finally, we assume that the generalized cost of transport may be written as follows:

T (ρ, r) = a 0 + ar -a ρ, ( 12
)
where a and a are the respective unit transport cost in the employment area and outside, and a 0 = ā0 + (a -a)ρ f with ā0 the fixed cost of transport, reflecting, among others, car ownership and parking costs.

In this specific framework, the system (10) becomes:

R (r) = -a λ H (r) H (r) = 1 β λR(r) Y -a 0 -ar+a ρω(r) , ( 13
)
which yields, after resolution, a closed form for r ω (ρ):

∀ρ ∈ [0, ρ f ], r ω (ρ) = A + B 1 - ρ ρ -C 1 - ρ ρ β , ( 14
)
where:

A = Y -a 0 a + ρ f a a λR A aN + 1 , B = β α ρ f a a λR A aN + 1 , C = A + B -r 0 , and ρ = ρ f λR A aN + 1 .

Decentralization of jobs and congestion as drivers of urban sprawl

As already mentioned in the introduction section, some authors (especially [START_REF] Fouchier | Les densités urbaines et le développement durable -Le cas de l'Île-de-France et des villes nouvelles[END_REF] have emphasized the fact that distance traveled by households and the time they spend traveling do not follow the same evolution with distance from the CBD. While traveled distance generally increases with distance to the city center, which constitutes a major source of adverse effects of urban sprawl, the time spent in transport tends to decrease, thus creating a strong incentive to sprawl.

The standard monocentric model is unable to account for this disjunction between distance traveled and time spent in transport. Within the framework of our model however, it is possible to analyze this phenomenon. For this, we study the traveled distance for home-to-work trips, denoted D(r) = rρ ω (r) and the transport cost5 , denoted T (r) = T (ρ ω (r), r). The respective evolutions of D and T with r are given by:

d D(r) dr = 1 -ρ ω (r), ( 15
) d T (r) dr = ρ ω (r) ∂ ∂ρ T (ρ ω (r), r) + ∂ ∂r T (ρ ω (r), r) = ∂T ∂ρ ρ ω (r) + ∂T /∂r ∂T /∂ρ . ( 16
)
We can assume that ρ ω < 1 for all r, which reflects observation that households are more dispersed than jobs, even only slightly [START_REF] Wheaton | Commuting, congestion and dispersal in cities with mixed land-use[END_REF]. That gives D (r) > 0. We also have ∂T /∂ρ ≤ 0 and ∂T /∂r ≥ 0. The sign of T (r) depends on value of ∂T /∂r ∂T /∂ρ .

We define two typical situations: the first one, called congested, is defined by |∂T /∂ρ| ≥ ∂T /∂r, meaning that the unit transport cost is higher in the center of the city than at the periphery. Thus, in this situation, ∂T /∂r ∂T /∂ρ ≥ -1. The second typical situation, called massified, is defined by |∂T /∂ρ| ≤ ∂T /∂r, meaning that the transportation system is very efficient at the center of the city allowing the unit transport cost to be lower there than at the periphery. Thus, in this situation, ∂T ∂T /∂ρ ≤ -1. Each of these typical situations for the transportation system may then lead to a typical configuration for the city. The massified situation leads to a configuration we called quasi-monocentric for distance traveled and time spent both increase with distance from the city center, as predicted by the standard monocentric model. On the contrary, the congested situation leads to a configuration we called deconcentrated, in which distance traveled increases with distance from the city center, while time spent decreases.

With the specific model described in subsection 2.5, the sufficient condition for the quasi-monocentric configuration is:

ρ f λR A aN + 1 ≤ β Y -a 0 -ar 0 min 1 a , 1 a , ( 17
)
which may be obtained for a not too high, a a and a small ρ f . As to the deconcentrated configuration, for a not too high to ensure

λR A λR A +aN α -β > 0, the sufficient condition is: αβ Y -a 0 -ar 0 a λR A λR A + aN α -β -1 ≤ ρ f λR A aN + 1 ≤ β Y -a 0 -ar 0 a ,
(18) which may be obtained for a a and a quite large ρ f . Fig. 1 illustrates this deconcentrated configuration for the following values of the parameters: a = 0.5 euro/km, a = 3 euros/km and ρ f = 5 km.

This shows that, when taking into account a spatial distribution of jobs, congestion in the center of the city may lead to a decrease in time spent in transport with distance from the city center, while distance traveled increases. This creates for households a potential incentive to sprawl, as they may thus benefit from a smaller transport cost. In the same time, the distance they travel increases, with the negative effects this may have on energy consumption, GES and local pollutant emissions. In terms of transport policy, this result suggests that policies consisting of capacity reduction in central areas without an efficient alternative transportation mode (rapid transit) potentially lead to more sprawl.

Transport-sober urban development

We now focus on the second question we raised in the introduction of this paper: is decentralization of jobs able to induce a transport-sober urban development? To answer this question, we compare two contrasted scenarios with a reference one. These scenarii differ by the level of decentralization of jobs: low, medium (for the reference) and high, the total number of jobs being fixed. Using our notations, it means that ρ f 1 < ρ f 0 < ρ f 2 , where index 0 stands for the reference scenario, index 1 for the low decentralization scenario, and index 2 for the high decentralization one.

We simultaneously consider two situations concerning the transport cost in the employment area: uncongested, meaning for instance that the use of several modes of transport allows for a low level of congestion in the city center, and congested, if car is the main mode of transport.

Surplus-neutral scenarios

As the location of jobs has an impact on the utility level of households, we want a means to compare comparable scenarios. When all households achieve the same utility level at equilibrium, as in De [START_REF] De Palma | Congestion pricing and long term urban form: Application to Île-de-France[END_REF] for instance, it may be obtained by adjusting income by a population tax, following the Herbert-Stevens model [START_REF] Fujita | Urban Economic Theory: Land Use and City Size[END_REF]. In our case, however, all households do not achieve the same utility level. thus compare scenarios which overall surplus of switching from one situation to the other is zero. We note σ(n) the surplus that a household indexed n may benefit when switching from a scenario i to a scenario j:

σ(n) = Y j -Y i + CV I(n) , ( 19
)
where CV I is the compensating variation of income. With a Cobb-Douglas utility function, CVI may be written as:

CV I(n) = T ρ j (n), r j (n) -Y i + Y j -T ρ j (n), r j (n) V i V j , ( 20
)
where V i and V j are the utility levels of the household in both scenarios. Injecting ( 20) in ( 19), we have the expression of surplus for household n:

σ(n) = Y j -T (ρ j (n), r j (n)) V j -V i V j . ( 21
)
Using the expression of the indirect utility function, and replacing T (ρ, r) by its expression, we obtain:

σ(n) = (Y j -Y i )+ (ρ f j -ρ f i ) a α λR A + a(N -n) aN λR A + aN λR A + a(N -n) α -1 . (22)
For the whole city population, we have:

S i→j = N 0 σ(n)dn = N (Y j -Y i )+ (ρ f j -ρ f i ) a α N 0 λR A + a(N -n) aN λR A + aN λR A + a(N -n) α -1 dn. (23)
The explicit calculation of the integral yields:

N 0 λR A + a(N -n) aN λR A + aN λR A + a(N -n) α -1 dn = η(α) -η(0) aN ≥ 0, (24)
where function η is defined by:

→ η( ) = - (λR A + aN ) a(2 -) (λR A + aN ) 2--(λR A ) 2-. ( 25
)
Finally, as we want to compare the scenarios on the same basis, we choose the income levels Y 1 and Y 2 for scenarios 1 and 2 so that S 0→1,2 = 0:

Y 1,2 = Y 0 + ρ f 0 -ρ f 1,2 N a α η(α) -η(0) aN (26)

Parameters choice

In order to compute numerical we to the values for the parameters of the model. Hence, we consider a city with 1, 700, 000 households, and a CBD of radius 1 km. This corresponds roughly to the Greater Paris region characteristics. The Cobb-Douglas utility function parameters are set to α = 0.72 and β = 0.28, meaning that about one-third of the income net of transportation cost is devoted to housing6 . The land capacity is set to λ = 6 km, the agricultural rent to 8 euros/m 2 /month, and the gross income in the reference scenario to Y 0 = 2, 500 euros/month. Concerning the transportation cost, the unit cost is set to a = 1 euro/km in the purely residential area, to a = 1 euro/km for the uncongested situation and to a = 2 euros/km in the congested situation, in the employment area. The 'fixed' part of the transportation cost, a 0 = ā0 + ρ f (a -a), is assumed independent of ρ f , through an adjustment of ā0 . As ā0 mainly accounts for parking costs (monetary and search cost) or, in the uncongested case, for park-and-ride cost or fixed transit cost, we simply assume that, in the congested case, these costs decrease with the decentralization level, while in the uncongested case, they keep their reference values.

The radius of the employment area is set to ρ f 0 = 5 km for the reference scenario, ρ f 1 = 2 km for the low decentralization scenario and ρ f 2 = 8 km for the high decentralization scenario. The gross income established by equation ( 26) and the values of ā0 are given in Tab. 1.

Results from numerical simulations

The main results of the simulations summarized Tab. 2, in which r f is the city radius, D the average distance traveled by households for home-to-work trips, T the average transportation cost they face. These results suggest that employment decentralization has a quite small impact on the city radius: the elasticity of city radius with respect to employment area radius is 0.01 when the city center is congested. However, in the same time, the distance traveled for home-to-work trips and the transportation cost strongly decrease: the elasticities are respectively -0.4 and -0.43. To put it another way, even though the model only deals with home-to-work trips, the findings suggest that compacification [START_REF] Melia | The paradox of intensification[END_REF] use the word intensification) may not have an important impact in terms of land use but an important negative impact in terms of distance traveled, and thus energy consumption and pollution, especially when the level of congestion is high.

Concerning the land rent and the lot sizes, we find a very small impact of the level of dencentralization on land rent, mainly because of the income adjustment, which allows households to keep the level of their bid almost unchanged. The lot size (and symmetrically the density) is impacted differently in the center and at the periphery, as illustrated by Fig. 2.

Spatial equity

In our model, the location of jobs has an impact on the distribution of utility levels through space. The workplace constitutes indeed an indirect income. To measure and compare the utility levels in monetary terms, we use the equivalent variation of income, or rather the equivalent income: the income that a reference household, for instance the household working at the very center of the CBD, i.e. in ρ = 0 and living in r 0 , should receive in order to achieve the utility level ρ of another household. This equivalent income is defined as:

Ỹρ=0 (ρ) = T ω (ρ = 0) + E(R ρ=0 , u ρ ), (27)
where E(R, u) is the expenditure function, giving the net income that allows a household facing a land rent of R to achieve utility level u.

With our Cobb-Douglas utility function, this expression becomes:

Ỹρ=0 (ρ) = Y + a α ρ f λR A + aN aN   1 -    λR A + aN 1 -ρ ρ f λR A + aN    α    . ( 28
)
Fig. 3 illustrates the result for the three scenarios. With a high level of employment decentralization, the equivalent income curves is steeper, meaning that the level of inequality is higher than with a low level of decentralization. Moreover, in the high decentralization scenario, households working close to the city center achieve a lower utility level than in the reference scenario, while the diagnostic is reverse for households working far from the center. This result suggests that the households working in the CBD and the the households working at the periphery of the employment area would probably not support the same policies toward jobs localization.

Conclusion

This paper examined the impact of job decentralization on urban sprawl through a twofold analysis. First we found, with a theoretical framework based on the monocentric urban model, that job decentralization associated with transport congestion may yield a strong incentive to sprawl for households: households living far from the city center benefit from a low transportation cost to their workplace, while the distance they travel is long, entailing energy consumption, GES and pollution emissions.

Second, through a comparison of three scenarios, we found that, other things being equal, a more centralized employment area leads to a moderately smaller city, the elasticity of the city radius with respect to employment area radius being 0.01, but a large variation in the distance traveled, the elasticity being -0.4. This suggests, along with authors such that [START_REF] Melia | The paradox of intensification[END_REF], that urban intensification is hardly able to reduce land consumption but may on the contrary induce long traveled distances for home-to-work trips.

Finally, within this theoretical framework, our findings qualify the widely held view that the compact city is the most efficient urban form. Our results suggest that a city where employment is partly decentralized may be efficient in terms of transport, without excessive land consumption.

Figure 1 :

 1 Figure 1: Transportation cost and home-to-work traveled distance in a deconcentrated city

Figure 2 :

 2 Figure 2: Lot size in the three scenarios.

Figure 3 :

 3 Figure 3: Equivalent income in the three scenarios.

Table 1 :

 1 Gross income and 'fixed fixed' part of transportation cost for the different situations and scenarios.

	Congested Scenario	Income (eu/m) ā0 (eu/trip)
		Reference	2, 500
	No	Low decent.	2, 526.5	3
		High decent.	2, 473.5
		Reference	2, 500	3
	Yes	Low decent.	2, 553	6
		High decent.	2, 447	0

Table 2 :

 2 Main results of the numerical simulations.

	Congested Scenario	r f (km) D (km/trip) T (eu/trip)
		Reference	18.24	6.51	9.51
	No	Low decent.	18.19	8.02	11.02
		High decent.	18.28	4.99	7.99
		Reference	17.92	6.30	11.80
	Yes	Low decent.	17.83	7.83	14.83
		High decent.	18.02	4.77	8.77

And hence on GES emissions and local pollution.

i.e. employers are indifferent to the residential location of their employees: this is a reductive hypothesis for the statistical distribution of incomes, since empirical research has shown the existence of an income gradient within cities[START_REF] Eberts | An empirical investigation of intraurban wage gradients[END_REF][START_REF] Mcmillen | Work location, residence location, and the intraurban wage gradient[END_REF][START_REF] Timothy | Intra-urban wage variations, employment location, and commuting times[END_REF]. Nevertheless,[START_REF] Glaeser | Decentralized employment and the transformation of the American city[END_REF] have shown that employment deconcentration tends to generate homogenization of incomes.

This property is shown in BL10 with a proof based on the Rule 2.3 inFujita (1989, p.28), a smaller ρ entailing a steeper bid rent for the household, and thus a residential location closer to the CBD.

This function R gives, for a given value R 0 of the rent at the beginning of the residential area, the value of the rent at the distance within which all households are located. The urban equilibrium is then characterized by R * 0 such that R(R * 0 ) = R A .

This cost may be seen as a generalized transport cost, in which the time cost constitutes the major part.

And about one-quarter of the gross income, as the transport budget is 14 % of the gross income in France[START_REF] Arthaut | Le budget transports des ménages depuis 40 ans -la domination de l'automobile s'est accrue[END_REF]. This figure of roughly 25 % is in line with[START_REF] Omalek | Le logement: une dépense importante pour les ménages franciliens modestes[END_REF].

* Support for this research was provided, in part by the Chaire ParisTech éco-conception, and in part by the French Ministry of Environment, Transportation and Housing. We remain fully responsible for the research and conclusions contained herein.


K Radzicki

S Bonelli

Impulse Response Function Analysis model application to the thermical seepage monitoring in the earth dams

L'érosion interne, liée avec la présence de fuites, est un danger bien connu pour la sécurité des barrages en terre et des digues. Actuellement, l'analyse thermométrique du remblais d'un ouvrage hydraulique est considérée comme la méthode la plus efficace et prometteuse pour identifier ces processus. Le papier présente une méthode appelée Analyse Thermique de la Fonction de Réponse Impulsionnelle (IRFTA) ainsi qu'un exemple de son application à l'analyse des mesures du température réalisés dans un barrage en terre. Le modèle IRFTA permet de décrire précisément le couplage transport de l'eau et de la chaleur dans le corps du barrage et en conséquence, il permet d'identifier des fuites mêmes minimales ainsi que leur variation.

Mots clefs : barrage en terre, fuite, surveillance thermique, reponse impulsionelle 1 Introduction

Nowadays one of the most effective and promising method which is used for leakages identification on the dams and dikes is the thermal analysis method [START_REF] Artieres | Active and Passive Defences against Internal Erosion in Assessment of the Risk Internal Erosion of Water Retaining Structures : Dams, Dykes and Levees[END_REF], [START_REF] Radzicki | Analyse retard des mesures de températures dans les digues avec application à la détection de fuites[END_REF], [START_REF] Artieres | Monitoring of Earthdams and stability with Fibre-Optics based monitoring system[END_REF]. Heat transport in the body of the earth hydraulic structure is described by energy Equation [START_REF] Artieres | Active and Passive Defences against Internal Erosion in Assessment of the Risk Internal Erosion of Water Retaining Structures : Dams, Dykes and Levees[END_REF]. The second and the third term of this equation describe respectively the conductive and advective heat transport processes, where advective process is defined as the transport of heat with the mass of flowing water. [START_REF] Johansson | Seepage monitoring in Embankment Dams[END_REF] 2

0 f T T T C v C t x x λ ∂ ∂ ∂ - - = ∂ ∂ ∂ (1)
where: T -temperature C -volumetric heat capacity of porous domain C f -volumetric heat capacity of water λ -thermal conductivity of porous domain Temperatures of the air and the water in the reservoir are the principal thermal loadings for the dam, supposing that other heat sources like geothermal and frozen processes, radiation and wind influence are neglected [START_REF] Johansson | Seepage monitoring in Embankment Dams[END_REF]. For the null water velocity there is only conductive, slow heat transport from the dam surfaces inward the dam. With rising of the water velocity, temperature from the reservoir is moved quicker with the masse of water. It results in the temperature field perturbation. Similarly, there are also thermal significant differences in the dike's body temperature between the zones of low and of faster seepage. Finally, an analysis of the temperature distribution in the dam's body allows for leakage identification. Moreover, temperature measurements can be realized with the fiber optic cable at every meter of its length 2

[1], [START_REF] Radzicki | Analyse retard des mesures de températures dans les digues avec application à la détection de fuites[END_REF], [START_REF] Johansson | Application of fiber-optics systems in embankment dams for temperature, strain and pressure measurements-some comparison and experiences[END_REF], [START_REF] Artieres | Monitoring of Earthdams and stability with Fibre-Optics based monitoring system[END_REF]. In consequence, this technology called distributed temperature sensing (DTS) gives a possibility of continuous monitoring of the structure in space. However, a correct analysis of the measured temperature is possible only with application of the suitable heat transport models. Particularly, it refers to temperature measurements performed in the downstream toe of dam, which zone is very attractive for dams thermal leakage monitoring. Firstly, in the most cases leakage path that cross the body of dam comes to this zone. Secondly, installation of the fiber optic cable in this zone is cheap and easy and can be realized on the existing dams. Nevertheless, a model to be applied to analyze thermal data from the downstream toe of dam musts take into consideration both external thermal loadings (air temperature or/and water temperature) as well the soil saturation degrees changes.

In this paper we present Impulse Response Function Thermal Analysis (IRFTA) model [START_REF] Artieres | Active and Passive Defences against Internal Erosion in Assessment of the Risk Internal Erosion of Water Retaining Structures : Dams, Dykes and Levees[END_REF], [START_REF] Radzicki | Analyse retard des mesures de températures dans les digues avec application à la détection de fuites[END_REF], [START_REF] Radzicki | Thermal seepage monitoring in the earth dams with impulse response function analysis model[END_REF] that allows for physical parameters identification of heat and water transport in the dams including leakage identification and its intensity estimation. This model is specially predicted to be used for analysis of the temperature data measured inside of the downstream toe of earth dam. After theoretical background, IRFTA model application results for dike of canal Oraison are described.

Background of IRFTA model

Energy Equation ( 1) describes a parabolic-linear problem. It means that heat transport (diffusional-advective) can be assumed to behave in the linear manner if the thermal porous medium properties and the water velocity are constant.

In consequence we can use the Green's function methodology to build a suitable model of the relevant problem. Using this approach the loading (input signal) a(t) and the system's response (output signal) y(t) are connected by the impulse response function of the system h(t) as follows:

( ) 0 ( ) ( ) ( ) ( ) t y t h a t h t a d τ τ τ = * = - ∫ (2)
where * is the mathematical convolution operator.

In other words, impulse response function describes how input signal (in form of Dirac delta) is modified by the porous zone of the dam. In our model we used an approximation of the impulse response function in the form of the two-parameters (α,η) exponential decay.

( ) ( , )

h t R α η ≈ (3)
The role of the parameters is explained by harmonic analysis. Under slowly varying loading conditions, η representing time-lag which quantifies the time elapsing between the onset of the loading and the response of the system in the point of measurements and α is the signal damping factor. Finally the IRFTA model has the following form :

( ) ( ) ( ) ( ) ( , ) , * , * C w w air air T x t R x t t R x t t θ θ θ = + + (4)
where: θ C -constant R w , R air -impulse response function approximation respectively for the water temperature and the air temperature loadings θ w , θ air -water temperature and air temperature loadings on dam surface Measured temperature T(x,t) is formed by superposition of responses of dam for the water temperature and the air temperature loadings, that are represented respectively in the second and the third term of Equation [START_REF] Radzicki | Thermal seepage monitoring in the earth dams with impulse response function analysis model[END_REF]. Finally the IRFTA model has four parameters. Two of them α w and η w inform about transformation of the thermal signal from the upstream face (water temperature loading). Downstream thermal signal (air temperature loading) modification is described by next two parameters α air and η air .

In particular conditions, IRFTA model can be applied in reduced form. If temperature sensor is located directly in the saturated zone of seepage and the air temperature influence is neglected we can use the following model :

( ) ( ) ( , ) , * C w w T x t h x t t θ θ = + (5)
Contrary, if there is non water temperature influence on the fiber optics temperature, it is possible to use model :

( ) ( ) ( , ) , * C air air T x t h x t t θ θ = + (6) 3
3 Analysis of dike's canal temperature with IRFTA model

Description of the site

An example of the application of the IRFTA model described in this paper focuses on thermal analysis of seepage process in about 27 metres height dike of the Oraison canal. Cross-section of this canal is presented at Figure 1. Bottom and slopes of the canal are covered by protection elements made from reinforced concrete slabs being simultaneously an impermeable layer. Fiber optic cable is situated at the land-side toe of the dike at the distance of 1000m. Next, it changes its location over a few tens of meters increasing at the top of the berm (Figure 2). For its entire length it is located at the depth of 0,8m below the soil surface.

Results of the analysis

Preliminary analysis showed that temperatures at the most part of the measurement points is influenced significantly by both air and water temperature. In consequence, a full 4-parameter IRFTA model defined with Equation (4) has been applied for modelisation. Reproduction of the data by model was excellent. For all measurements points a coefficient of determination R2 was higher then 0,99 (for the values of function 1-R2 consequently is lower than 0,01). Analysis of the IRFTA model parameters values and their variations allowed to identify several hydro-thermal zones of the dike in relation to different seepage intensity. Due to limit place in the article we describe only three, chosen zones. First zone is localized between 700 and 815 m of the fiber optic cable. Results of the analysis for this zone are presented at the Figure 3. For the air temperature influence on the fiber optic temperature, we observe that values of α air parameter varies from 0,8 to 0,7. It means respectively only from 20% to 30% signal damping and in consequence significant influence of the air temperature. Calculated time-lag of the response of the fiber temperature for the air temperature loading (η air ) was about 20 to 30 days. In the Figure 4 we see that time-lag between the air temperature maximum value and the fiber optic temperature maximum value is similar. The same we observe for the minimum values. This confirms a correct estimation of time-lag values by IRFTA model.

4

FIG. 3 -Values of the IRFTA models parameters for temperature analysis between 700 and 815m distance of the fiber optic cable.

FIG. 4 -Time-lag of the response of the fiber optic temperature for air temperature influence.

On the other side, the influence of water temperature is very week. Damping of water temperature signal equals from 80% to 90% (α w varies from 0,2 to 0,1) and the time of response (η w ) is from 14 to 27 days. These values of IRFTA models parameters exclude important seepage process in this zone of the dike.

The second presented zone is localized on the distance from 880 m to 950 m of the fiber optic. Results of the modelisation are presented at the Figure 5. Comparing the values of the parameters between previous zone and this one we see that for signal transformation of the water temperature influence, damping is little lower particularly for the distance from 880m to 925m where α w equals from 0,2 to 0,25 (damping varies from 80% to 75%). Simultaneously transport of the heat is faster for the water temperature signal. η w varies between 9 and 20 days. Simultaneously, also the air temperature signal transport is stronger (α air equals from 0,76 to 0,91) and faster (η air equals from 11 to 20 days) than as the case of the previous zone. These changes in air temperature signal transformation parameters have clear physical explication. Due to more significant seepage process, humidity zone around the seepage zone is larger and degree of saturation of the soil in this zone is higher. It results in changes in values of thermal parameters of porous medium as volumetric heat capacity and thermal conductivity. Finally it causes acceleration of conductive heat transport process, which is stronger and faster, also from land-side of the dike (air temperature influence) inward the dam. However, even in spite of the fact of more significant seepage presence in the discussed zone, values of the IRFTA model parameters exclude existence of the leakage problem.

5

FIG. 5 -Values of IRFTA models parameters for temperature analysis between 880 and 950m distance of the fiber optic cable.

In the third presented zone fiber optic cable (from 990 to 1100 m) is situated close to the crest of the berm. In consequence, this cable is outside of seepage influence. Results of the analysis for this zone are presented at the Figure 6. We observe only significant air temperature influence (α air varies from 0,8 to 0,96). For water temperature influence, α w values are closed to 0 and time-lag η w seeks to very large values of some thousands days. Physically it means a null water temperature influence. In some singular points because of convergence difficulty in real data reproduction these values are different. However α w even there is lower then 0,1.

FIG. 6 -Values of IRFTA models parameters for temperature analysis between 990 and 1100 m distance of the fiber optic cable.

Conclusion

In the paper we showed the application of IRFTA model to the thermal analysis of hydraulic field of the dike of the canal. We see that with this method, very small seepage process influence related only to degree of humidity variation can be easy detected. Moreover, physical definition of the model parameters allows for estimation of the seepage filtration intensity. Analysis of the model's parameters values and their variation gives also a possibility of clear physical interpretation of the observed thermal-hydraulic processes, not possible to be performed with only statistical models. It is particularly important for the earth dams and dikes of the canals behaviours assessments. In consequence, thermal monitoring and analysis of temperature with IRFTA model increase of hydraulic 6 structure security level. It minimizes also the cost of eventual reparation work of the dam due to fact that erosion process (linked with leakage process) is early and precisely defined. IRFTA model application is particularly predicted for analysis of the temperature measured with fiber optic cable. This technology allows for continuous monitoring of the leakages process along the dike. Installation of the fiber optics cable in the downstream (land-side) toe of earth dam or dike is easy and cheap and can be recommended as efficient monitoring system of existing structures.

FIG. 1 -FIG. 2 -

 12 FIG. 1 -Cross-section of Oraison canal.