Internal stabilization of the Oseen-Stokes equations by Stratonovich noise

Viorel Barbu (Iaşi, Romania)

Introduction and statement of the problem

This work is concerned with internal stabilization via Stratonovich noise feedback controller of Oseen-Stokes system (1)

∂X ∂t -ν∆X + (f • ∇)X + (X • ∇)g = ∇p in (0, ∞) × O, ∇ • X = 0 in (0, ∞) × O, X(0, ξ) = x(t), ξ ∈ O, X = 0 on (0, ∞) × ∂O.
Here, O is an open and bounded subset of R d , d = 2, 3, with smooth boundary ∂O and f, g ∈ C 2 (O; R d ) are given functions. In the special case g ≡ 0, system (1) describes the dynamic of a fluid Stokes flow with partial inclusion of convection acceleration (f •∇)X (X is the velocity field). The same equation describes the disturbance flow induced by a moving body with velocity f through the fluid. Should we mention also that in the special case f ≡ g ≡ X e , where X e is the equilibrium (steady-state) solution of the Navier-Stokes equation ( 2)

∂X ∂t -ν∆X + (X • ∇)X = ∇p + f e , ∇ • X = 0, X| ∂O = 0,
and f e ∈ C(O; R d ), system [START_REF] Arnold | Stochastic Differential Equations; Theory and Applications[END_REF] is the linearization of (2) around X e . In this way, the stabilization of (1) can be interpreted as the first order stabilization procedure of steady-state Navier-Stokes flows.

Our aim here is to design a stochastic feedback controller of the form

(3) u = 1 O 0 M k=1 R k (X) • βk , R k ∈ L((L 2 (O)) d ),
which stabilizes in probability system (1) and has the support in an open subdomain O 0 of O.

Here, {β k } M k=1 is a system of mutually independent Brownian motions on a probability space {Ω, F, P} with filtration {F t } t>0 while the corresponding closed loop system (4)

dX -ν∆X dt + (f • ∇)X dt + (X • ∇)g dt = 1 O 0 M k=1 R k (X) • dβ k + ∇p dt X(0) = x
is taken in Stratonovich sense (see, e.g., [START_REF] Arnold | Stochastic Differential Equations; Theory and Applications[END_REF]) and this is the significance of the symbol R k (X) • βk in the expression of the noise controller (3). We have denoted by

1 O 0 the characteristic function of the open set O 0 ⊂ O.
In [START_REF] Barbu | Stabilization of Navier-Stokes Flows[END_REF], [START_REF] Barbu | The internal stabilization by noise of the linearized Navier-Stokes equation[END_REF], [START_REF] Barbu | Stabilization of a plane periodic channel flow by noise wall normal controllers[END_REF], [START_REF] Barbu | Internal stabilization of the Navier-Stokes equation[END_REF], the author has designed similar stabilizable Ito noise controllers for equation [START_REF] Arnold | Stochastic Differential Equations; Theory and Applications[END_REF] and related Navier-Stokes equations. However, it should be said that, with respect to Ito noise controllers, the Stratonovich feedback controller (3) has the advantage to be stable with respect to smooth changes βε k of the noise βk and this fact is crucial not only from the conceptual point of view, but also for numerical simulations and practical implementation into system (1) of the random stabilizable feedback controller

u ε (t) = 1 O 0 M k=1 R k (X(t)) βε k (t),
where βε k is a smooth approximation of β ε k . As regards the literature on stabilization of linear differential systems by Stratonovich noise, the pioneering works [START_REF] Arnold | Stabilization by noise revisited[END_REF], [START_REF] Arnold | Stabilization of linear systems by noise[END_REF] should be primarily cited. For linear PDEs, this procedure was developed in [START_REF] Caraballo | Comparison of the long time-behavior of linear Ito and Stratonovich partial differential equations[END_REF], [START_REF] Caraballo | Stabilization of linear PDEs by Stratonovich noise[END_REF] which are related to this work. For general results on internal stabilization of Navier-Stokes systems with deterministic feedback controllers, we refer to [START_REF] Barbu | Internal stabilization of Navier-Stokes equations with finitedimensional controllers[END_REF].

The noise stabilizing feedback controller

Consider the standard space of free divergence vectors

H = {y ∈ (L 2 (O)) d ; ∇ • y = 0 in O, y • n = 0 on
∂O} and denote by A 0 : D(A 0 ) ⊂ H → H the realization of the Oseen-Stokes operator in this space, that is, (5)

A 0 y = P (-ν∆y + (f • ∇)y + (y • ∇)g), y ∈ D(A 0 ),
where

D(A 0 ) = H ∩ (H 2 (O)) d ∩ (H 1 0 (O)) d .
Here, P is the Leray projector on H and H 2 (O), H 1 0 (O) are standard Sobolev spaces on O. In the following, it will be more convenient to represent equation (1) in the complex Hilbert space H = H + iH by extending A 0 to A :

D(A) ⊂ H → H via standard procedure, A(x + iy) = A 0 x + iA 0 y.
The operator A has a countable set of eigenvalues {λ j } ∞ j=1 (eventually complex) with the corresponding eigenvectors ϕ j . Denote by A * the adjoint of A with eigenvalues λ j and eigenvectors ϕ * j . Each eigenvalue is repeated in the following according to its algebraic multiplicity m j . We denote by N the number of eigenvalues λ j for which [START_REF] Barbu | The internal stabilization by noise of the linearized Navier-Stokes equation[END_REF] Re

λ j > 0 for j ≤ N, λ 1 + λ 2 + • • • + λ N > 0.
(In the above sequence, each λ j is taken together its conjugate λ j and, clearly, there is such a natural number N .)

Set X u =lin span{ϕ j } N j=1 and denote by X s the algebraic complement of X u in X . It is well known that X u and X s are both invariant for A and, if we set

A u = A 0 | X u , A s = A 0 | X s , we have for their spectra σ(A u ) = {λ j } N j=1 , σ(A s ) = {λ j } ∞ j=N +1 and, since -A s is the generator of an analytic C-semigroup e -A s t in X s , we have (7) e -A s t L(H) ≤ C exp(-Re λ N +1 t), t ≥ 0,
(see, e.g., [START_REF] Barbu | Stabilization of Navier-Stokes Flows[END_REF], p. 14). In the following, we shall assume that (i) All the eigenvalues {λ j } N j=1 are semisimple.

This means that the algebraic multiplicity of each λ j , j = 1, ..., N , coincides with its geometric multiplicity or, in other words, the finitedimensional operator (matrix) A u is diagonalizable. As we will see later on, this assumption is not essentially necessary but it simplifies however the construction of the stabilizing controller because it reduces the unstable part of the system to a diagonal finite-dimensional differential system. In particular, it follows by (i) that we can choose the dual systems {ϕ j } and {ϕ * j } in such a way that (8) We note that the uncontrolled Oseen-Stokes system (1) can be rewritten in the space H as [START_REF] Caraballo | Comparison of the long time-behavior of linear Ito and Stratonovich partial differential equations[END_REF] dX dt

ϕ i , ϕ * j = δ ij , i, j = 1, ...,
+ AX = 0, t ≥ 0, X(0) = x,
and setting X u = N j=1 y j ϕ j , X s = (I -P N )X, where P N is the algebraic projector on X u , we have

dX u dt + A u X u = 0, X u (0) = P N x, ( 10
)
dX s dt + A s X s = 0, X s (0) = (I -P N )x. (11) We set A u = { Aϕ j , ϕ * k } N j,k=1 = diag λ j N j=1
and so, by ( 8), we may rewrite [START_REF] Caraballo | Stabilization of linear PDEs by Stratonovich noise[END_REF] where γ 0 > 0. The matrix C k is explicitly constructed in [START_REF] Arnold | Stabilization of linear systems by noise[END_REF] and it will be used below to construct a stabilizable feedback controller of the form (3). Namely, we set in (3)

in terms of y = y j = X u , ϕ * j N j=1 as (12) dy dt + A u y = 0, y(0) = P N x, ϕ * j N j=1 . Since Tr(-A u ) = -λ 1 -λ 2 -• • • -λ N < 0, it
(15) R k (X) = N i,j=1 C k ij X, ϕ * j φ i , k = 1, ..., M, where C k ij N i,j=1 = C k , ( 16
) φ i = N =1 α i ϕ * , i = 1, ..., N,
and α i are chosen in such a way that (17

) N =1 α i γ j = δ ij , i, j = 1, ..., N.
Here, γ j = O 0 ϕ * ϕ * j dξ and since, by the unique continuation property (see [START_REF] Barbu | Stabilization of Navier-Stokes Flows[END_REF], p. 157), the eigenfunction system {ϕ * j } is linearly independent on O 0 , we infer that the matrix γ j N ,j=1 is not singular and, therefore, there is a unique system {α i } which satisfies (17). Then, by ( 16), we see that

(18) 1 O 0 φ i , ϕ * j = δ ij , i, j = 1, ..., N.

As mentioned earlier, O 0 is an arbitrary open subset of O.

Theorem 1 is the main result.

Theorem 1 The solution X to the closed loop system (4), where R k are defined by (15), is exponentially stable in probability, that is,

(19) |X(t)| H ≤ Ce -γt |x| H , ∀t ≥ 0, P-a.s.,
where γ > 0.

The proof of Theorem 1

The idea of the proof, already used in stabilization theory of infinitedimensional systems with finite-dimensional unstable subspaces, is to stabilize the finite-dimensional system (10) by a feedback controller of the form (3) and to reconstruct consequently the system via the infinite-dimensional stable complement (11).

Namely, taking into account ( 10), (11), we write the closed loop system (4), that is,

(20) dX + AX dt = P 1 O 0 M k=1 R k (X) • dβ k as (21) dX u + A u X u dt = P N M k=1 N i,j=1 C k ij X u , ϕ * j P (, φ i ) • dβ k , ( 22
) dX s + A s X s dt = (I -P N ) M k=1 N i,j=1 C k ij X u , ϕ * j P (φ i ) • dβ k ,
where X u + X s = X. Taking into account (18) and that X u = N j=1 y j ϕ j , we may rewrite (21) as ( 23)

dy + λ y dt = M k=1 C k j y j • dβ k , = 1, ..., N,
and so, by (19), we have that (We have denoted by C several positive constants independent of t.)

As regards the existence of a solution X s to (22), this is standard and follows from the general theory of linear infinite-dimensional stochastic equations. Now, in order to estimate X s , it is convenient to replace (22) by its Ito formulation (see, e.g., [START_REF] Caraballo | Stabilization of linear PDEs by Stratonovich noise[END_REF])

(25) dX s +A s X s dt= 1 2 M k=1 [P (1 O 0 R k ] 2 X u dt+ M k=1 P (1 O 0 R k (X u ))dβ k .
Taking into account that, by ( 7), e -A s t is an exponentially stable semigroup on X s , without loss of generality we may assume that Re A s x, x ≥ γ|x| 2 H . (Otherwise, proceeding as in [START_REF] Barbu | The internal stabilization by noise of the linearized Navier-Stokes equation[END_REF], we replace the scalar product x, y by Qx, y where Q is the solution to the Lyapunov equation A s Q + QA * s = γI.) Then, applying Ito's formula in (25), we see that (26)

d|X s (t)| 2 H +2γ|X s (t)| 2 H dt= X s (t), M k=1 [P (1 O 0 R k ] 2 X u (t) dt + M k=1 |P (1 O R k (X u ))| 2 dt + 2 M k=1 P (1 O R k (X u ))X s dβ k .
Taking into account that, by (24), |X u (t)| H ≤ e -γt |X u (0)|, we infer by (26) that

E|X s (t)| 2 H ≤ C|X s (0)| 2 e -γt , ∀t ≥ 0. Since M (t) = t 0 M k=1 P (1 O 0 R k (X u )X s dβ k is a local martingale, it
follows by (26) and Lemma 3.1 in [START_REF] Barbu | Internal stabilization of Navier-Stokes equations with finitedimensional controllers[END_REF] that

|X s (t)| H ≤ C|X s (0)|e -γt , ∀t ≥ 0, P-a.s.,
which completes the proof.

The design of a real valued noise controller

A nice feature of the feedback controller given by Theorem 1 is its simple structure. Moreover, its computation does not rise any special problem because the stabilizing matrices C k for the diagonal system (23) can be explicitly expressed.

If all λ j , j = 1, ..., N , are real, then the feedback controller {R k } given by Theorem 1 as well as the closed loop system (4) are real, too. However, if some λ j are complex, then the above feedback controller does not stabilize system (1), but its complex realization in the space H = H + iH. In this case, though the stabilizing feedback controller (3) is quite simple, its implementation into real systems rises some delicate problems because the feedback controller is in implicit form as function of Re X, Im X. In order to circumvent these inconvenience, we will derive below from the above construction a real valued noise feedback controller which has a stabilizing effect on the Oseen-Stokes system [START_REF] Arnold | Stochastic Differential Equations; Theory and Applications[END_REF].

To this end, we denote by {ψ j } N j=1 the orthogonalized system (via Schmidt procedure) {Re ϕ j , Im ϕ j , j = 1, ..., N }. Taking into account that each ϕ j in the above system arises together with its conjugate, the dimension of the system {ψ j } is again N . We also set

X re u = lin span{ψ j } N j=1
and note that X re u = {Re y; y ∈ X u } and H = X re u ⊕ X re s , where X re s = {Re y; y ∈ X s }. Moreover, the operator A 0 leaves invariant spaces X re u , X re s and if we set

A re u = A 0 | X re u , A re s = A 0 | X re s , we have, of course, (27) e -A re s t L(H) ≤ Ce -γ 0 t , ∀t > 0,
for some γ 0 > 0. It is also easily seen that Tr[-A re u ] = Tr[-A u ] < 0. Then, we define as in (15) the operators

(28) R k (x) = N i,j=1 C k ij X, ϕ j φ i , k = 1, ...., M,
where

C k = C k ij N
i,j=1 is the matrix system which stabilizes in probability, via Theorem 2 in [START_REF] Arnold | Stabilization of linear systems by noise[END_REF], the finite-dimensional system

dy + A re u y dt = M k=1 C k y • dβ k .
Here, A re u = A re u ψ i , ψ j N i,j=1 and φ i is given by

φ i = N =1 α i ψ , i = 1, ...., N,
where α i are chosen in such a way that

N =1 α i γ j = δ ij , γ j = O 0 ψ ψ j dξ.
As in the previous case, the system {ψ j } N j=1 is still linearly independent on O 0 which implies that det γ j = 0 and the above system has a unique solution α i .

Then, arrived at this point, the proof of Theorem 1 applies neatly to conclude that the feedback noise controller (3) with R k instead of R k is stabilizing in the probability system [START_REF] Caraballo | Comparison of the long time-behavior of linear Ito and Stratonovich partial differential equations[END_REF].

Namely, we have the following stabilization result. Remark 3 It should be emphasized that there is a close connection between the unique continuation property of eigenfunctions ϕ j (or ϕ * j ) of the Oseen-Stokes operator and the above construction of a stabilizing finite-dimensional feedback controller (3). In fact, as seen above, the design of u in the form (15) or ( 28) is essentially based on this sharp property through existence in the algebraic system (17).

Remark 4 One might speculate that the noise feedback controller

(29) u = 1 O 0 M k=1 R k (X -X e ) • dβ k
inserted in the right hand side of Navier-Stokes system (2) stabilizes exponentially in probability the equilibrium solution X e to (2). In general, this might not be true, but it happens for Ito noise of the form (29) for any x in a sufficiently small neighborhood of X e (see [START_REF] Barbu | Internal stabilization of the Navier-Stokes equation[END_REF]) and one might expect that the fixed point argument used there for the equivalent random system is still applicable in the present case. We expect to give details in a later work.

Conclusions

We have designed in this paper a Stratonovich stochastic feedback controller which exponentially stabilizes in probability a general Oseen-Stokes system from fluid dynamics. The controller has the support in an arbitrary open subset O 0 of the velocity field domain O ⊂ R d , d = 2, 3, and has a finite-dimensional linear structure which involves the dual eigenfunctions corresponding to unstable eigenvalues of the systems. The stabilization effect is independent of the Reynold number 1/ν through the dimension N of the stabilizing controller (3) might depend on ν.

  N. (Here, and everywhere in the following, •, • stands for the scalar product inH and H. By | • | H and | • | H we denote the corresponding norms.)

  follows by Theorem 2 in[START_REF] Arnold | Stabilization of linear systems by noise[END_REF] that there is a sequence of skew-symmetric matrices {C k } M k=1 , where M = N -1 such that the solution y to the Stratonovich stochastic system (13)dy + A u y dt = M k=1 C k y • dβ k , t ≥ 0,has the property (14) |y(t)| ≤ C|y(0)|e -γ 0 t , P-a.s., ∀t > 0,

  (24) |y (t)| ≤ Ce -γt |y (0)|, P-a.s., ∀t > 0, = 1, ..., N.

  Theorem 2 The solution X to the closed loop system (4) with R k = R k defined by (28) satisfies (19) in the real norm |X(t)| H .


Maria Do

Rosário De Pinho

NECESSARY CONDITIONS FOR CONSTRAINED OPTIMAL CONTROL PROBLEMS

This is possible if one assumes g and b to be continuously. differentiable, for example.  Extensions to state constrained problems and problems with both state and mixed constraints.

 Determination of NC for general Free Time problems

  Necessary Conditions for State Constrained Problems 2. Refining Necessary Conditions  weak assumptions and/or stronger conclusions  easily translated in computational schemes  special interest on state constrained problems 3. Necessary Conditions for various non standard problems  Hybrid optimal control, Free time problems,  infinite horizon problems, nonregular mixed constrained problems 5. Relations between sets of Necessary Conditions 6. Applications vs Theory CHALLENGES ON NC FOR COCP NONSMOOTH PROBLEMS Necessary conditions under "minimal assumptions" Based on Clarke's monograpgh 2005.If G is pseudo Lipchitz, then it satisfies the bounded slope condition.NEW DEVELOPMENTSThis result is a main tool in the derivation of refinements of necessary conditions for various problems.For example:(New)Nonsmooth Maximum Principle For standard optimal control problems-for convex-linear case it is a sufficient and necessary condition  For general mixed constrained problems including problems with DAE  For problems with state constraints ? Necessary conditions for general free time problems to have some MULTIPLIERS in the Euler inclusion instead of the Normal cone.



  For mixed (state-control) constraintsIdentify other situations where multipliers can stated in terms of multipliers and not the normal cone to the constrained set.

SPECIAL CASEs

In many problems we have state-control constraints with ´´two'' controls:

Then the bounded slope condition is imposed only on S(t).

PROBLEMS WITH DIFFERENTIAL AND ALGEBRAIC EQUATIONS

Assuming some smoothness…

Regularity Assumption Drawback:

Question: Any alternative regularity assumption for the bounded slope condition and w/o Lipschitz continuity wrt u?

Answer: in some cases, YES. Roughly speaking we get NC assuming, for that the following holds but we also need to assume that COVERING CONDITION


Mohamedou Ould Haye

Anne Philippe

Marginal density estimation for linear processes with cyclical long memory

Keywords: Confidence band, empirical process, limit theorem, mean integrated squared error

 to the stationary processes for which the singularities of the spectral density are not limited to the origin. We show that the convergence rates and the limiting distribution may be different in this context.

1 Introduction [START_REF] Hosking | Fractional differencing[END_REF] introduced long memory processes with quasi periodic behaviour. This fact corresponds, for stationary processes, to spectral densities which exhibit singularities at non zero frequencies. Many authors have contributed to the construction of fractional models with singularities/poles outside the origin, see for instance, [START_REF] Gray | On generalized fractional processes: a correction[END_REF][START_REF] Gray | On generalized fractional processes[END_REF]; [START_REF] Hassler | mis)specification of long memory in seasonal time series[END_REF]; [START_REF] Viano | Long-range dependence and mixing for discrete time fractional processes[END_REF]; [START_REF] Leipus | Modelling long-memory time series with finite or infinite variance: a general approach[END_REF]; [START_REF] Bisognin | Properties of seasonal long memory processes[END_REF].

We can distinguish between two types of long memory: one regular and the other cyclical according to whether the spectral density has a pole at the origin or outside the origin. From a statistical point of view, the estimators of the long memory parameter have been adapted to yield some estimates if cyclical effects are assumed. In a parametric context, the √ n-consistency of the maximum likelihood estimate or the Whittle estimate has been proved (see [START_REF] Hosoya | A limit theory for long-range dependence and statistical inference on related models[END_REF]; [START_REF] Giraitis | Gaussian estimation of parametric spectral density with unknown pole[END_REF] when the pole is unknown). Semi parametric estimates can be more or less easily adapted to the cyclical case (see [START_REF] Hidalgo | Estimation of the location and exponent of the spectral singularity of a long memory process[END_REF]; Arteche andRobinson (2000, 1999); [START_REF] Hsu | Semiparametric estimation for seasonal long-memory time series using generalized exponential models[END_REF]; [START_REF] Reisen | Estimating seasonal long-memory processes: a Monte Carlo study[END_REF]; [START_REF] Whitcher | Wavelet-based estimation for seasonal long-memory processes[END_REF]).

When we consider empirical process related statistics, the situation is more delicate. The normalisation and the limit distribution can be different according to whether the memory is regular or cyclical. An important literature is devoted to the convergence of the empirical process, see for instance [START_REF] Ho | On the asymptotic expansion of the empirical process of longmemory moving averages[END_REF]; [START_REF] Giraitis | Central limit theorem for the empirical process of a linear sequence with long memory[END_REF] in regular case and Ould Haye (2002) Ould [START_REF] Ould Haye | A noncentral limit theorem for the empirical process of linear sequences with seasonal long memory[END_REF] in cyclical case.

In this paper we give some convergence results on the kernel estimator of the marginal density f . Let (X 1 , • • • , X n ) be an observed sample from f , the kernel estimator of f is defined by fn (x) = 1 nm n n j=1 K

x -X j m n .

(1.1)

where the bandwidth m n is a sequence such that m n → +∞ and nm n → 0 as n → ∞, and K is a kernel function. Consider the following infinite moving average process,

X t = t j=-∞ b t-j ξ j , t ≥ 1 (1.2)
where

• the sequence (b k ) k has the form b k = k -(α+1)/2 j∈J a j cos kλ j + o(1) , k → ∞ (1.3)
where α ∈ (0, 1) and λ j = 0 for all j ∈ J, a finite non empty subset of N.

• (ξ n ) n is a sequence of independent and identically distributed random variables with zero mean and finite variance Eξ 2 0 = σ 2 < ∞.

From [START_REF] Giraitis | A generalized fractionally differencing approach in long-memory modeling[END_REF], the covariance function r of (X t ) defined by (1.2) and (1.3) has the form r(h) = h -α j∈J a j cos hλ j + o(1) .

(1.4) as h tends to infinity.

A large class of linear processes satisfying these conditions is obtained by filtering a white noise (ξ i )

X t = G(B)ξ t with G(z) = g(z) m j=-m 1 -e iλj z (αj -1)/2 , m ≥ 1, (1.5)
where B is the backshift operator and where g is an analytic function on {|z| < 1}, continuous on {|z| ≤ 1} and g(z) = 0 if |z| = 1, and where 0 < α j ≤ 1, α j = α -j , λ -j = -λ j , j = 0, . . . , m, and

0 = λ 0 < λ 1 < . . . < λ m < π.
Taking α = min{α j , j = 0, . . . , m}, and J = {j ≥ 0 :

α j = α}, if α < α 0 /2 then the condition (1.3) is satisfied.
Note that the condition on the coefficient α ensures that ∞ h=1 |r(h)| = ∞, thus the process has a long-memory. However this condition is not enough to characterize the cyclical long memory.

1. When α < α 0 /2. h j=1 r(j) = o   h j=1 r(j) 2   as h → ∞. Therefore the process (X 2
t ) has also a long memory, which is more persistent than (X t ) (see Remark 1 for the exact expressions). This fact characterises cyclical long memory, and the asymptotic behavior of many statistics (see below for the empirical process) can be drastically different. We focus on this case in this paper.

2. When α > α 0 /2, the cyclical behavior is less persistent than the regular long memory (singularity at frequency zero). The presence of singularities outside zero do not modify the convergence results obtained in the regular case.

3. When α = α 0 /2, both (X t ) and (X 2 t ) will contribute to the limiting distribution, which will be a combinaison of the two previous cases.

We consider the empirical process associated with the process (X n ) n≥1 defined by

F n (x) = 1 n n j=1 1 {Xj ≤x} .
Ould [START_REF] Ould Haye | A noncentral limit theorem for the empirical process of linear sequences with seasonal long memory[END_REF] proved the following results for the linear process (X n ) defined in (1.5). Assume that Eξ 4 0 < ∞, the cumulative distribution function of ξ 0 is 5 times differentiable with continuous bounded and integrable derivatives on R. Denote

d n = n 1-α , and D = (2 -2α)(1 -2α) 4Γ(α) cos(απ/2) .
If α < α 0 /2, then, as n tends to infinity, we have

d -1 n [nt] F [nt] (x) -F (x) =⇒ F ′′ (x) 2 R(t), (1.6)
where R is a linear combination of independent Rosenblatt processes with the same parameter α

R(t) = R α,Λ (t) = D -1 j∈J c j R
(1)

j (t) + R (2) j (t) , (1.7)
where Λ = {λ j , j ∈ J}, and where

• c 0 = h 0 /2, c j = h j if j = 0 and h j = g(e iλj ) k =j 1 -e i(λ k -λj ) (α-1)/2 ,
• R (i) j (t), i = 1, 2 and j ∈ J are Rosenblatt processes with parameter 1 -α, independent except for j = 0, R

(1)

0 (t) = R (2) 0 (t).
The paper is organized as follows. In Section 2, we establish a limit theorem for the kernel estimate. This extends one of [START_REF] Ho | On the asymptotic expansion of the empirical process of longmemory moving averages[END_REF]'s results, in particular we show the contribution and the effect of the singularities of the spectral density outside the origin to the convergence rate and the limiting distribution. Then we apply our limit theorem to construct confidence bands for the density function.

Similarly to [START_REF] Hall | Convergence rates in density estimation for data from infinite-order moving average processes[END_REF]; [START_REF] Mielniczuk | On the asymptotic mean integrated squared error of a kernel density estimator for dependent data[END_REF], we provide in Section 3, the asymptotic behavior of the mean integrated squared error, and we show that the equivalence, one had in regular long memory' can be modified when the singularities of the spectral density are not limited to the origin.

Asymptotic distribution of the kernel estimator

Hereafter, we assume that the kernel K is a continuous function with compact support and K(x)dx = 1. Concerning the bandwidth m n , we assume that m n → 0 and nm n → ∞, as n tends to infinity. The equality fn

(x) -E fn (x) = 1 m n R K x -u m n d F n (u) -F (u) (2.8)
clearly shows the relationship between the estimate fn (x) and the empirical process F n (x). The process fn (x) is sometimes called the empirical density process.

For every integer n ≥ 1, we define the following statistics

Y n,1 = n k=1 X k , Y n,2 = n k=1 s<r b r b s ξ k-s ξ k-r , (2.9) and S n,2 (x) = n F n (x) -F (x)) + F ′ (x)Y n,1 - 1 2 F ′′ (x)Y n,2 .
(2.10)

Remark 1 For linear processes defined in (1.5), the following equivalences as n tends to infinity, have been proved by Ould Haye and Philippe ( 2003)

Var(Y n,2 ) ∼ 1 4 Var n j=1 (X 2 j -E(X 2 1 )) ∼ Cn 2-2α . (2.11)
and

Var(Y n,1 ) = Var n j=1 X j ∼ Cn 2-α0 .
(2.12) Therefore (2.12) and (2.11) imply that the convergence rate obtained in Proposition 2.1 is smaller than the convergence rate of Xn .

Let us define the class of Parzen kernels of order s.

Definition 2.1 A kernel function K is said to be a Parzen kernel of order s ≥ 2 if it satisfies the following conditions

1. R K(u)du = 1, 2. for every 1 ≤ j ≤ s -1, R u j K(u)du = 0, 3. R |u s ||K(u)|du < ∞.
Bretagnolle and Huber (1979) proved the existence of such kernels, for which, an explicit construction can be found in [START_REF] Gasser | Kernel estimation of regression functions[END_REF].

Proposition 2.1 Consider a process (X n ) defined in (1.5). Assume that α < α 0 /2, Eξ 4 0 < ∞, the cumulative distribution function of ξ 0 is 5 times differentiable with continuous bounded and integrable derivatives on R. Let K be a Parzen kernel of order 4 having bounded total variation. Assume that the bandwidth has the form

m n = n -δ , where α 4 < δ < α 2 .
Then, as n tends to infinity

n α sup x∈R | fn (x) -f (x)| d -→ sup x∈R f ′′ (x) 2 |R α,Λ |. (2.13) where R α,Λ = R α,Λ (1). Moreover, n α ( fn (x) -f (x)) C b (R) =⇒ - f ′′ (x) 2 R α,Λ , (2.14)
where

C b (R)
=⇒ denotes the convergence in C b (R), the space of continuous bounded functions.

Proof: The difference between fn and f can be expressed as

fn (x) -f (x) = fn (x) -E fn (x) + E fn (x) -f (x) = 1 m n K(u)d F n (x -m n u) -F (x -m n u) + f (x -m n u) -f (x) K(u)du.
We first replace F n -F by its expression in (2.10). Then we apply the integration by parts formula on the first integral. For the second, we apply the Taylor-Lagrange formula. There exists a real number u * such that |u * -x| < |m n u| and

fn (x) -f (x) = -1 nm n S n,2 (x -m n u)dK(u) + Y n,1 n f ′ (x -m n u)K(u)du - Y n,2 n f ′′ (x) K(u)du + Y n,2 n m n f (3) (u * )uK(u)+ + -m n uf ′ (x) + m 2 n u 2 2 f ′′ (x) - m 3 n u 3 6 f (3) (x) + m 4 n u 4 24 f (4) (u * ) K(u)du =:a n (x) + b n (x) + c n (x) + d n (x) + e n (x).
Now, a proof similar to that of Theorem 2.2 in [START_REF] Ho | On the asymptotic expansion of the empirical process of longmemory moving averages[END_REF] allows us to write for 2δ

< α n α+δ-1 sup x∈R |S n,2 (x)| a.s.
-→ 0, as n → ∞.

(2.15)

And, thus we have

n α sup x∈R |a n (x)| P -→ 0, as n → ∞ (2.16)
where P -→ denotes the convergence in probability. For the sequences b n (x), d n (x), e n (x), we get the same convergence in probability as in (2.16) by bounding the variances. To obtain the bounds, we start from the variances of Y n,1 and Y n,2 defined in (2.12) and (2.11), and we use the fact that K is a Parzen kernel and that f is 4 times differentiable and bounded derivatives. We get, as n tends to infinity,

Var(n α sup x∈R |b n (x)|) ≤ n 2α-2 Var Y n,1 sup x∈R |f ′ (x)| |K(u)|du = Cn 2α-2 n 2-α0 = Cn 2α-α0 -→ 0, Var(n α sup x∈R |d n (x)|) ≤ n 2α-2 Var Y n,2 m n sup x∈R |f (3) (x)| |uK(u)|du = Cn 2α-2 n 2-2α n -δ -→ 0, n α sup x∈R |e n (x)| ≤ sup x∈R |f (4) (x)| n α-4δ 24 u 4 |K(u)|du = O(n α-4δ ) -→ 0,
These four convergences in probability imply that both sequences

n α sup x∈R | fn (x) -f (x)| and n α sup x∈R |f ′′ (x)| Y n,2 n = n α sup x∈R |c n (x)|
have the same limit as n tends to infinity. According to Lemma 2.1 in Ould Haye and Philippe (2003), this common limit is equal to

sup x∈R f ′′ (x) 2 |R α,Λ |.
Hence (2.13) is proved. According to (2.11), we notice that the rate n -α given in (2.13) is the convergence rate of n -1 n j=1 (X 2 j -E(X 2 1 )). Similarly, as n tends to infinity, the finite-dimensional distributions of

n α ( fn (x) -f (x)) and -n α f ′′ (x) Y n,2 n = n α c n (x)
converge simultaneously to the finite-dimensional distributions of -(f ′′ (x)/2)R α,Λ . This concludes the proof of (2.14) because (2.13) implies the tightness of n α ( fn (x) -f (x)).

Remark 2 We clearly see that the choice of the class of Parzen kernels allows the bias E fn (x)f (x) to become negligible. If K is not a Parzen kernel, the contribution of the bias e n (x) is not negligible with respect to b n (x). Therefore, (2.13) is false for a standard kernel unless we replace fn (x) -f (x) by fn (x) -E fn (x) in (2.13).

Remark 3 The result (2.13) in Proposition 2.1 can be applied to obtain a goodness of fit test on the marginal density.

Remark 4 The result (2.13) in Proposition 2.1 provides confidence bands for f which depend on the derivative f ′′ . In general, f ′′ is not available, and thus the confidence band cannot be calculated. Then f ′′ can be replaced by its kernel estimate given by

f ′′ n (x) = 1 nm 3 n n j=1 K ′′ x -X j m n .
(note that it is necessary to assume that the kernel function K is twice differentiable.)

Proposition 2.2 Under the same hypotheses as in Proposition 2.1 and if the kernel function K is twice differentiable and its derivative K ′′ is continuous, then for each interval [a, b] on which f ′′ is positive, we have

2n α sup x∈[a,b] fn (x) -f (x) f ′′ n (x) d -→ |R α,Λ |.
(2.17)

In other words, as n tends to infinity, for every t > 0, we have

P fn (x) - t f ′′ n (x) 2n α ≤ f (x) ≤ fn (x) + t f ′′ n (x) 2n α , a ≤ x ≤ b → P |R α,Λ | < t .
(2.18)

In Proposition 2.3, we give a consistent estimate of the quantiles of process R α,Λ . Using (2.18), this allows us to obtain asymptotic confidence band for the density f (x) which is valid for every

x ∈ [a, b]. Proof : Let φ be the function defined on C b (R) by φ(g) = sup x∈[a,b] g(x) f ′′ (x)
Since φ is continuous, (2.14) ensures the following convergence :

2n α sup x∈[a,b] fn (x) -f (x) f ′′ (x) d -→ |R α,Λ |, as n → ∞. (2.19)
Now, we prove that the difference

Y n (x) := n α fn (x) -f (x) f ′′ (x) - fn (x) -f (x) f ′′ n (x) satisfies sup x∈[a,b] |Y n (x)| P -→ 0, as n → ∞.
This convergence is obtained as follows. We rewrite Y n (x) as

|Y n (x)| = n α fn (x) -f (x) f ′′ (x) f ′′ n (x) -f ′′ (x) f ′′ n (x)
.

and by (2.19), it is enough to prove that

sup x∈R f ′′ n (x) -f ′′ (x) f ′′ n (x) P -→ 0, as n → ∞.
(2.20)

The difference between f ′′ n and f ′′ can be written as

f ′′ n (x) -f ′′ (x) = -1 nm 3 n S n,2 (x -m n u)dK ′′ (u) + Y n,1 n f (3) (x -m n u)K(u)du- - Y n,2 n f (4) (x -m n u)K(u)du + f ′′ (x -hu) -f ′′ (x) K(u)du.
by replacing f with f ′′ and fn with f ′′ n and following the same lines as the proof of Proposition 2.1. Then, we get sup

x∈R | f ′′ n (x) -f ′′ (x)| = O n -(2δ∧(1-3δ)) .
Since 0 < δ < 1/4, we have sup

x∈R | f ′′ n (x) -f ′′ (x)| P -→ 0, moreover, the derivative f ′′ satisfies inf x∈[a,b] |f ′′ (x)| > 0.
Thus, we get (2.20). This concludes the proof.

Proposition 2.3 Fix β ∈ (0, 1). Let c(α, Λ, β) be the quantile of order β of the process R α,Λ defined in (1.7). If (α n , Λ n ) are consistent (in probability) estimators of (α, Λ). then

c(α n , Λ n , β) P → c(α, Λ, β) (2.21)
Remark 5 In the references given in the introduction, the parametric and semi parametric methods provide estimators of (α, Λ) which satisfy the condition required in Proposition 2.3.

Proof : We want to show (2.21) which will be obtained if we show that the application (γ, θ) → c(γ, θ, β) is continuous, as (α n , Λ n ) P → (α, Λ). To prove this continuity we prove that the mappings g, h below are continuous,

((0, 1) × [0, π] |J| , |.|) g → (C b (R), . ) h → ((0, 1), |.|),
where . is the uniform metric, and in the following decomposition F γ,θ is the distribution function

of R γ,θ . (γ, θ) → [g(γ, θ) = F γ,θ ] → [h(F γ,θ ) = c(γ, θ, β)].
Continuity of g can be proved as follows. Consider a deterministic sequence (γ n , θ n ) such that (γ n , θ n ) → (γ, θ) as n → ∞. Then to prove that F γn,θn → F γ,θ uniformly it will be enough to show that R γn,θn =⇒ R γ,θ . To obtain the latter weak convergence it will suffice to show that every sequence of Rosenblatt variables (R γn ) with parameter γ n converges weakly to a Rosenblatt variable R γ with parameter γ, as R γn,θn is a linear combination of independent Rosenblatt variables R γn with the coefficients c j /D that are continuous functions of γ n , θ n . We have from [START_REF] Major | Multiple Wiener-Ito Integrals[END_REF])

R γn = R 2 e i(x+y) -1 i(x + y) W n (dx, dy)
where W n (dx, dy) = |x| (γn-1)/2 |y| (γn-1)/2 W (dx, dy)

with W (dx, dy) being the standard Gaussian random measure, and since

|x| (γn-1)/2 |y| (γn-1)/2 → |x| (γ-1)/2 |y| (γ-1)/2
then we have the required convergence. Now to prove the continuity of h it is enough to note that the quantile function is continuous (with respect to the uniform metric) over the class of monotonic continuous distribution functions, i.e. if

F n -F → 0 then h(F n , β)) → h(F, β).
Of course here we do have F γn,θn -F γ,θ → 0, as we just established that R γn,θm =⇒ R γ,θ .

3 Asymptotic mean integrated squared error (MISE)

The mean integrated squared error (MISE) of the estimate fn is defined by

R E fn (x) -f (x) 2 dx.
For a wide class of linear processes including the processes with short and regular long memories, [START_REF] Hall | Convergence rates in density estimation for data from infinite-order moving average processes[END_REF] and [START_REF] Mielniczuk | On the asymptotic mean integrated squared error of a kernel density estimator for dependent data[END_REF] studied the asymptotic behavior of the MISE. In particular, they established the following equivalence, when n tends to infinity,

R E fn (x) -f (x) 2 dx ∼ R E 0 fn (x) -f (x) 2 dx + Var(X n ) R f ′ (x) 2 dx (3.22)
where E 0 denotes the expectation with respect to the distribution of n independent random variables distributed from the density f . In particular, the equivalence (3.22) shows that the convergence rate of the MISE cannot be faster than the convergence of Var(X n ). In other words, the convergence rate of the kernel density estimates is bounded from above by the convergence rate of the empirical mean. This is the optimal rate. Hereafter, we assume that the distribution of the innovation (ξ k ) satisfies

[Z] There exist δ > 0 and C < ∞ such that the characteristic function of ξ 0 satisfies

|Ee iuξ0 | ≤ C(1 + |u|) -δ (3.23)
Theorem 3.1 Let (X n ) be a linear process defined in (1.2) and (1.5) such that the distribution of ξ 0 satisfies [Z] and Eξ 4 0 < ∞. Assume that α < 1 3 ∧ α0 2 and the kernel K is a bounded symmetric density function. Then the MISE satisfies, as n tends to infinity,

M ISE( fn ) ∼ R E 0 fn (x) -f (x) 2 dx + 1 4 Var 1 n n j=1 (X 2 j -E(X 2 1 )) R f ′′ (x) 2 dx (3.24)
where E 0 denotes the expectation with respect to the distribution of n independent random variables distributed from the density f . [START_REF] Ould Haye | A noncentral limit theorem for the empirical process of linear sequences with seasonal long memory[END_REF]). Equation (3.24) shows that this term is a ceiling rate of MISE independently of the choice of the kernel and bandwidth.

Remark 6 The variance Var

1 n n j=1 (X 2 j -E(X 2 1 )) is also equivalent to 4 Var( 1 n Y n,2 ) (see Ould

Proof.

Notation : for an arbitrary function g, we denote by ĝ its Fourier transform.

The proof consists in adapting the proof of [START_REF] Mielniczuk | On the asymptotic mean integrated squared error of a kernel density estimator for dependent data[END_REF] to the cyclical case. Using [START_REF] Hall | Convergence rates in density estimation for data from infinite-order moving average processes[END_REF] decomposition of the MISE, we have

MISE( fn ) = R E 0 fn (x) -f (x) 2 dx+ + 1 nπ n-1 j=1 (1 -j/n) | K(m n t)| 2 Re (E(e it(X1-Xj+1) ) -| f (t)| 2 dt (3.25) := MISE 0 + W n
Let f j be the joint density of (X 1 , X j+1 ). We extend the expansion of f j , obtained by [START_REF] Giraitis | Asymptotic normality of regression estimators with long memory errors[END_REF] for the first order, to the second order as follows: there exists a function ℓ j : R 2 → R such that

f j (x, y) = f (x)f (y) + r(j)f ′ (x)f ′ (y) + 1 2 r(j) 2 f ′′ (x)f ′′ (y) + ℓ j (x, y) ∀(x, y) ∈ R 2 (3.26)
where r is given in (1.4).

We have E(e it(X1-Xj+1) ) = e it(x-y) f (x)f (y) dx dy + r(j) e it(x-y) f ′ (x)f ′ (y) dx dy+ + 1 2 r(j) 2 e it(x-y) f ′′ (x)f ′′ (y) dx dy + e it(x-y) ℓ j (x, y) dx dy

= | f (t)| 2 + r(j)| f ′ (t)| 2 + 1 2 r(j) 2 | f ′′ (t)| 2 + ℓ j (t, -t).
(3.27)

Similarly to [START_REF] Mielniczuk | On the asymptotic mean integrated squared error of a kernel density estimator for dependent data[END_REF], W n in (3.25) can be written as

W n = 2 n n-1 j=1 (1 -j/n)r(j) |K mn ⋆ f ′ | 2 (t) dt + 1 n n-1 j=1 (1 -j/n)r(j) 2 |K mn ⋆ f ′′ | 2 (t) dt+ + 1 nπ n-1 j=1 (1 -j/n) | K(m n t)| 2 Re ℓ j (t, -t) dt where K mn (x) = m -1 n K(xm -1 n )
, and where f ⋆ g is the convolution of f and g. Moreover we have, for k = 1, 2,

|K mn ⋆ f (k) | 2 (t) dt = f (k) (t) 2 dt + o(1), n → ∞.
We obtain

W n = 2 n n-1 j=1 (1 - j n )r(j) f ′ (t) 2 dt + o(1) + 1 n n-1 j=1 (1 - j n )r(j) 2 f ′′ (t) 2 dt + o(1) + + 1 nπ n-1 j=1 (1 - j n ) | K(m n t)| 2 Re ℓ j (t, -t) dt. (3.28)
According to [START_REF] Giraitis | A limit theorem for polynomials of a linear process with long-range dependence[END_REF], we have

Var 1 n n j=1 (X 2 j -E(X 2 1 )) = 2 n 2 1≤i,j≤n r 2 (i -j) + O(n -1 ), = 2 n 2 (nr(0) + 2 n-1 j=1 (n -j)r(j) 2 ) + O(n -1 ) = 4 n n-1 j=1 (1 -j/n)r(j) 2 + O(n -1 ) := γ(n) (3.29)
Moreover, using the form of r given in (1.4) and the fact that α < 1/3, we get

γ(n) = 4 n n-1 j=1 (1 -j/n)j -2α ( k∈J a k cos jλ k + o(1) ) 2 + O(n -1 ) = 2 n n-1 j=1 (1 -j/n)j -2α k∈J a 2 k + O(n -1 ) = 2 n n 1-2α 1 1 -2α - 1 2 -2α k∈J a 2 k + O(n -1 ) = n -2α 1 (1 -2α)(1 -α) k∈J a 2 k + O(n -1 ) ∼ Cn -2α (3.30)
As α < 1 3 ∧ α0 2 and using (2.12), we get

2 n n-1 j=1 (1 -j/n)r(j) = 1 n 2 Var(Y n,1 ) -r(0)n -1 = O(n -α0 ) + O(n -1 ) = o(n -2α
W n = 1 4 Var 1 n n j=1 (X 2 j -E(X 2 1 )) f ′′ (t) 2 dt + o(n -2α )+ + 1 nπ n-1 j=1 (1 -j/n) | K(m n t)| 2 Re ℓ j (t, -t) dt.
Since r(j) 2 behaves asymptotically as j -2α , and

1 nπ n-1 j=1 (1 -j/n) | K(m n t)| 2 Re ℓ j (t, -t) dt ≤ 1 nπ n-1 j=1 (1 -j/n) | ℓ j (t, -t)| dt (3.32)
the proof is completed using the following lemma proven below.

Lemma 3.1 Under the same assumption of Theorem 3.1, .33) for ǫ an arbitrary positive number smaller than 1 -3α 10 .

| ℓ j (t, -t)| dt = O(j -2α-ǫ ). ( 3
Proof of Lemma 3.1 By definition of ℓ j in (3.26), we have

ℓ j (x, y) = f j (x, y) -f (x) f (y)(1 -xyr(j) + 1 2 x 2 y 2 r(j) 2 ) We split the integral R | ℓ j (t, -t)| dt = |t|>j ǫ | ℓ j (t, -t)| dt + |t|<j ǫ | ℓ j (t, -t)| dt (3.34)
where ǫ is an arbitrary positive number smaller than 1 -3α 10 .

Under assumption (3.23), [START_REF] Giraitis | Asymptotic normality of regression estimators with long memory errors[END_REF] proved for the regular long memory that for arbitrary k

| f j (x 1 , x 2 )| ≤ c(k)(1 + |x|) -k for all x = (x 1 , x 2 ) ∈ R 2 and | f (x)| ≤ c(k)(1 + |x|) -k for all x ∈ R.
Their proof can be adapted to the cyclical case i.e. when the coefficients (b j ) j∈N satisfies (1.3). Using their notation, it suffices to construct a finite set J 1 such that for all j ∈ J 1 : |b -j | > 2|b t-j | + c 1 where c 1 does not depend on t. Since (|b j |) j∈Z is not summable, there exists a subsequence (j u ) u∈Z such that b -ju = 0. We can take J 1 a subset of {j u : u ∈ Z} with [δ|J 1 |] = k + 3. Indeed, for j ∈ J 1 , we have |b -j | > C(J 1 )|j -(α+1)/2 |, and for t large enough there exists c1

|j|

-(1+α)/2 > 2/C(J 1 )|t -j| -(1+α)/2 + c1 .
Therefore, there exists c 1 such that for all j ∈ J 1 ,

|b -j | > 2|b t-j | + c 1 .
For all k ′ , the first integral in (3.34) satisfies

|t|>j ǫ | ℓ j (t, -t)| dt ≤ j -ǫk ′ |t|>j ǫ |t| k ′ | ℓ j (t, -t)| dt = O(j -ǫk ′ ).
Therefore we can take any arbitrary k ′ such that k ′ > (2α + ǫ)/ǫ. For the second integral in (3.34), it is enough to show that sup

|u|<j ǫ | ℓ j (u)| = O(j -2α-2ǫ ). (3.35)
The proof is quite similar to that of equation (2.20) in Giraitis et al [START_REF] Giraitis | Asymptotic normality of regression estimators with long memory errors[END_REF] adding the terms of order two in the expansion.

We write the difference f j (x, y) -f (x) f (y) from products of the characteristic function φ of ξ 1 .

f j (x, y) -f (x) f (y) = I1 I1 I1 φ(xb -i + yb t-i ) - I1 I1 I1 φ(xb -i )φ(yb t-i ) := a 1 a 2 a 3 -a ′ 1 a ′ 2 a ′ 3 = (a ′ 1 -a 1 )a 2 a 3 + (a ′ 2 -a 2 )a ′ 1 a 3 + (a ′ 3 -a 3 )a ′ 1 a ′ 2
where I 1 = {|i| < j 2ǫ }, I 3 = {|t -i| < j 2ǫ } and I 3 = Z -(I 1 ∪ I 2 ). We will deduce (3.35) from |a i | < 1, |a ′ i | < 1 and the following facts, for all u < t ǫ a i -a ′ i = O(j -2α-2ǫ ), i = 1, 2 (3.36) a 3 -a ′ 3 = a ′ 3 (-xyr(j) + 1 2 x 2 y 2 r(j) 2 ) + O(j -2α-2ǫ ).

(3.37)

Similarly to [START_REF] Giraitis | Asymptotic normality of regression estimators with long memory errors[END_REF], we prove (3.36) with i = 1 (or similarly for i = 2) as follows

|a 1 -a ′ 1 | ≤ |i|≤j 2ǫ |φ(xb -i + yb j-i ) -φ(xb -i )φ(yb j-i )| ≤ |i|≤j 2ǫ |xb j-i |
As |i| ≤ j 2ǫ and x ≤ j ǫ , we have

|xb j-i | ≤ Cj ǫ j -(1+α)/2 = j -2α-4ǫ O(1)
since ǫ < 1 -3α 10 . Therefore |a 1 -a ′ 1 | = j -2α-2ǫ O(1). To prove (3.37), we follow the same calculations as [START_REF] Giraitis | Asymptotic normality of regression estimators with long memory errors[END_REF] x 2 y 2 r(j) 2 + O(j -2α-2ǫ )).

  page 325. Since |xb -i | + |yb j-i | = o(1), we write a 3 -a ′ 3 of the forma 3 -a ′ 3 = a ′ 3 (e Qj (x,y) -1) = a ′ 3 (Q j (x, y) + 1 2 Q j (x, y) 2 + o(Q j (x, y) 2 ))whereQ j (x, y) = i∈I3 Ψ(xb -i , yb j-i ) = -xy i∈I3 b -i b j-i + O( i∈I3 (xb -i ) 2 |yb j-i | + |xb -i ||yb j-i | 2 ) := -xy i∈I3 b -i b j-i + R n and Ψ(x, y) = log(φ(x + y)) -log(φ(x)) -log(φ(y))and we show thatQ j (x, y) = -xyr(j) + O( I1∪I2 |x||y||b -i ||b j-i | + i x 2 |y||b -i | 2 |b j-i |) = Q j (x, y) 2 = x 2 y 2 r(j) 2 + x 2 y 2 ( i∈I1∪I2 b -i b j-i ) 2 -2x 2 y 2 -xy i∈Z b -i b j-i i∈I1∪I2 b -i b j-i + R 2 n -2R n xy i∈I3 b -i b j-i For |x| < j ǫ et |y| < j ǫ we have I1∪I2 |x||y||b -i ||b j-i | = j 2ǫ-(1+α)/2 O(1) = j -2α-2ǫ O(1) since ǫ < (1 -3α)/8 and i x 2 |y||b -i | 2 |b j-i | = j -α/2-1/2+3ǫ O(1) = j -2α-2ǫ O(1)since ǫ < (1 -3α)/10. These asymptotic behaviors ensure that for |x| < j ǫ et |y| < j ǫ we have a 3 -a ′ 3 = a ′ 3 (xyr(j) + 1 2

Acknowledgement

The authors would like to thank the anonymous referee for their helpful comments and suggestions, that improved the presentation of the paper.


Alexandros Bousios

David Waxman

Stephen R Pearce

Absence of close-facing retrotransposons: A comparison of molecular data and theory

Retrotransposons occur in extremely large numbers in many eukaryotic genomes. However, little is known of the factors which affect the distribution of close proximity retrotransposons. In this work we investigate the frequency of close facing retrotransposons in a plant species with extremely high numbers of retrotransposons. Molecular observations are compared with predictions of a mathematical model that assumes a uniform probability of retrotransposon insertion into the genome. The mathematical model plays the role of a null hypothesis. We find that compared with the predictions of the model, there is a statistically significant deficit of identical copies of facing retroelements that are close to one another. This suggests that an efficient mechanism exists that removes or limits close facing retroelements.

Retroelements, which include retroviruses and retrotransposons, are selfish genetic elements that are present in all eukaryotic genomes [START_REF] Coffin | Retroviruses[END_REF]. Apart from causing major illnesses, such as HIV/AIDS, retroelements comprise a substantial and changing component of eukaryotic genomes, and have a profound effect on the composition and evolution of all eukaryotic species [START_REF] Feschotte | Plant transposable elements: Where genetics meets genomics[END_REF][START_REF] Kazazian | Mobile elements: drivers of genome evolution[END_REF][START_REF] Wicker | Genome-wide comparative analysis of copia retrotransposons in Triticeae, rice and Arabidopsis reveals conserved ancient evolutionary lineages and distinct dynamics of individual copia families[END_REF].

In this work we focus solely on the LTR (long terminal repeat) retrotransposons [START_REF] Kumar | Plant retrotransposons[END_REF][START_REF] Vitte | LTR retrotransposons and flowering plant genome size: emergence of the increase/decrease model[END_REF]. These mobile genetic elements are a subset of retroelements that are considered to be defective retroviruses as they lack components which enable them to escape the cell and be infective. Despite this, they retain the retrovirus replicative mode of transposition, which enables them to reinfect the nuclear genome in which they reside [START_REF] Willhelm | Reverse transcription of retroviruses and LTR retrotransposons[END_REF]. Successive rounds of insertion therefore lead to an increase in retrotransposon copy number, with extremely large numbers of retrotransposons being found in many eukaryotic genomes (Pearce et al., 1996a;[START_REF] Sanmiguel | Nested retrotransposons in the intergenic regions of the maize genome[END_REF][START_REF] Meyers | Abundance, distribution and transcriptional activity of repetitive elements in the maize genome[END_REF].

Although retrotransposons are dispersed throughout the genome, a direct implication of a high copy number of such elements, or a particular preference for insertion into particular regions of the genome, is that the number of closely spaced elements will be high. In order to investigate the distribution of closely spaced elements, we adopted a molecular technique, IRAP (Inter-Retrotransposon Amplified Polymorphism), which is able to detect two copies of the same retrotransposon when they are facing each other and reside within the range of PCR [START_REF] Kalendar | IRAP and REMAP: Two new retrotransposon-based DNA fingerprinting techniques[END_REF][START_REF] Kalendar | IRAP and REMAP for retrotransposon-based genotyping and fingerprinting[END_REF]. We have formulated a simple mathematical model to provide concrete predictions of the molecular data.

A c c e p t e d m a n u s c r i p t 3

The system adopted as the source of empirical data was required to have high copy numbers of transpositionally active retroelements. We chose the plant, Agave tequilana, for this purpose, as this organism has abundant and dynamic retrotransposon populations with very high copy numbers [START_REF] Bousios | Isolation and characterization of Ty1-copia retrotransposon sequences in the blue agave (Agave tequilana Weber var. azul) and their development as SSAP markers for phylogenetic analysis[END_REF]. This organism is a good subject for this study, as the total retroelement load is likely to be greater than many species studied including rice [START_REF] Mccarthy | Long terminal repeat retrotransposons of Oryza sativa[END_REF] and Arabidopsis thaliana [START_REF] Pereira | Insertion bias and purifying selection of retrotransposons in the Arabidopsis thaliana genome[END_REF]). Within A. tequilana we restricted analysis to two examples of the major classes of retrotransposons: Ty3-gypsy and Ty1-copia elements which both occur within this organism in very large numbers [START_REF] Bousios | Isolation and characterization of Ty1-copia retrotransposon sequences in the blue agave (Agave tequilana Weber var. azul) and their development as SSAP markers for phylogenetic analysis[END_REF]. This last aspect makes these elements in this organism a source of highly reliable statistics, and hence ideally suited to the present study.

In the work described here only an outward facing primer from the 5'LTR was used, so only close proximity elements in a 5' facing orientation will produce an IRAP product and be detected (figure 1).

Figure 1 The number and size of close-proximity elements observed in the IRAP data reveal the number of adjacent, 5' facing retrotransposons and the distances between them. This allows us to test the hypothesis that retrotransposons have randomly inserted into the genome with no preference to be located near other elements of the same type, and no preference to insert in a particular genomic location. The hypothesis is introduced into this work in the form of a mathematical model which explicitly assumes a uniform probability of insertion anywhere within the genome.

Such a probability of insertion does not lead to regularly spaced retrotransposons, since statistical fluctuations will automatically lead to variation of the distances between adjacent retrotransposons. The molecular data differs significantly from the predictions of the model (which assumes a uniform probability of insertion) and indicates that close facing retrotransposon insertions are extremely rare or perhaps even absent in the genome. This may occur by either a universal mechanism to prevent close facing retrotransposon insertions or an efficient mechanism to remove retrotransposons with this orientation.

Results

Determining retrotransposon copy number

In order to compare observations and the predictions of the mathematical model introduced here, we require an estimate of the retrotransposon copy number, which is denoted n. We used a molecular technique to estimate n, namely quantitative Southern hybridization (Pearce et al. 1996a); see methods and Appendix A for details of this determination. This estimate is accurate within a factor of 2 so the actual copy numbers have a value, n, that lies within the range to .

Modelling retrotransposon insertions

In Furthermore, when / nr L is small ( / 1 nr L << ) as is relevant to this study, the distribution ( ) f r is approximately constant:

( ) / f r n L ≈ (1)
and hence behaves as a uniform distribution.

The retrotransposons that are actually detected by PCR are those that lie within a finite range of one another; a distance we denote D. The value of D adopted in this work is bp (we use the abbreviation bp for base pairs) which is the size limit of the products of PCR (data not shown). The statistical nature of the dynamics of retrotransposons means that the number of IRAP bands observed in a particular genome is a random variable, which we denote B. When copy number n is large ( 1 n >> ) and the ratio / nD L is small ( / 1 nD L << ), as it is for our data sets, the number of IRAP bands is predicted from the model to have a mean or expected value

(written E[B]) of 2 [ ] /(4 ) E B n D L ≈ (2)
(see Appendix C for details). From analysis of the experimental data we arrived at the estimates in Table 1 for the number of observed IRAP bands, B obs . We note that the overall number of observed IRAP bands, B obs , for both types of retrotransposon, is close to the lower theoretical prediction for the expected number of bands, E[B]. This finding is suggestive that the overall number of observed IRAP bands is consistent with the estimated copy number although it is at the lower end of the range expected from a uniform probability of insertion. This is not surprising since a "real" population of retrotransposons (in contrast to an idealised theoretical population) will have a degree of sequence variability in their terminal sequences (through mutation) which will prevent amplification by PCR. Although this sequence variability will reduce the observed number of IRAP bands (B obs ) it will not influence the overall form of the size distribution of IRAP products.

In figure 2 we have plotted the theoretically predicted size distribution of IRAP fragments within the A. tequilana genome and compared this directly with the size profiles observed using IRAP for Teq1 and Teq17 retrotransposons.

Figure 2

Size distribution of IRAP products

The frequencies of IRAP products in the 0-1500bp size range that were observed for the two A.

tequilana retrotransposons Teq1 and Teq17 were compared with the size distribution predicted in the 0-1500bp size range from random insertions of an element with similar copy numbers (denoted U in Figure 2). The IRAP experiments for both types of retroelements show a slightly greater number of fragments in the middle of the size range (800bp) with lower levels than those predicted by the model at both the higher and lower ends of the size range. We make the assumption that the number of fragments in the middle of the size range (800bp) as being indicative of a "true" uniform distribution of fragment sizes. The mismatch of the observed frequency of 800bp

fragments and the frequency of the theoretical distribution, labelled U, in Figure 2, can be taken as a measure of the accuracy of the estimate of copy number, n. The reduction in observed numbers of upper size fragments (>800bp) is not unexpected, given that longer PCR products are more difficult to amplify under the experimental conditions. The completely unanticipated observation is that IRAP products in the lower size range of the experiment are effectively absent. These 'effectively absent' fragments lie within the optimal size range for PCR amplification. This constitutes strong evidence that facing, closely spaced adjacent retroelements (i.e. those separated by less than 300bp) are an extreme rarity in this genome. The deficit of retroelements with small separations is statistically significant within the framework of the null hypothesis (i.e., under the assumption of a uniform probability of insertion). Obtaining the observed numbers of closely adjacent retroelements, given a uniform probability of insertion, occurs with probability .

Discussion

Eukaryotic genome structure is complex and dynamic and no part of these genomes is as dynamic as the repetitive component [START_REF] Morgante | Plant genome organization and diversity: the year of the junk![END_REF]. The distribution and copy number of

A c c e p t e d m a n u s c r i p t

8 retrotransposons observed in eukaryotes is the outcome of evolution, i.e., a complex history of amplification, mutation and deletion, resulting in profiles of retrotransposon insertion which are unique to particular lineages [START_REF] Vitte | Analysis of retrotransposon structural diversity uncovers properties and propensities in angiosperm genome evolution[END_REF]. The high number of repetitive sequences in eukaryotic genomes indicates that these make a significant contribution to genome size. In humans an expansion in genomic size was an ancient event, with transposition declining over the past 40 million years [START_REF] Lander | Initial sequencing and analysis of the human genome[END_REF]. By contrast, in many plants a genomic expansion has occurred in the last few million years [START_REF] Sanmiguel | The Paleontology of intergene retrotransposons of maize[END_REF][START_REF] Bennetzen | Transposable element contributions to plant gene and genome evolution[END_REF][START_REF] Vicient | Copia-like retrotransposons in the rice genome: few and assorted[END_REF].

For species with larger genomes, gene densities appear to be especially non-random, with the general picture emerging that genomes are organised into gene-dense regions separated by heterochromatic gene-poor blocks [START_REF] Moore | Cereal genome evolution: Grasses, line up and form a circle[END_REF]. Although it would be expected that retrotransposon insertion would be tolerated in non-coding regions of the genome [START_REF] Bennetzen | Transposable element contributions to plant gene and genome evolution[END_REF][START_REF] Wicker | Analysis of a contiguous 211kb sequence in diploid wheat (Triticum monococcum L.) reveals multiple mechanisms of genome evolution[END_REF], in-situ hybridisations show that retroelements are present (with minor exceptions) throughout the entire genome (Pearce et al., 1996b;[START_REF] Waugh | Genetic distribution of BARE-1-like retrotransposable elements in the barley genome revealed by sequence-specific amplification polymorphisms (S-SAP)[END_REF]. On a finer scale, and in particular where the distribution of individual retroelements is studied in more detail, the distribution of elements is less uniform with extremely high densities of retrotransposons existing in particular genomic locations [START_REF] Vitte | LTR retrotransposons and flowering plant genome size: emergence of the increase/decrease model[END_REF]. A number of studies have shown that insertion is common in coding sequences. In particular, in human and mouse genomes, many transposable elements are found in protein-coding genes [START_REF] Nekrutenko | Transposable elements are found in a large number of human protein-coding genes[END_REF] with up to 25% of promoter regions of human genes containing transposable element sequences [START_REF] Jordan | Origin of a substantial fraction of human regulatory sequences from transposable elements[END_REF]. The discovery of retroelements influencing the expression of genes in wheat [START_REF] Kashkush | Transpositional insertion of retrotransposons alters the expression of adjacent genes in wheat[END_REF] and the preference of Tos17 retrotransposon for gene-rich regions of rice [START_REF] Miyao | Target site specificity of the Tos17 retrotransposon shows a preference for insertion within genes and against insertion in retrotransposon-rich regions of the genome[END_REF] show that the interaction of retroelements with genes is common. We can therefore see that although retroelements are present in almost every part of

A c c e p t e d m a n u s c r i p t

9 the genome, each particular retroelement may have its own preferred location and although may not be abundant in the whole genome may be abundant locally [START_REF] Vitte | LTR retrotransposons and flowering plant genome size: emergence of the increase/decrease model[END_REF].

Considering the apparent non-random distribution of retroelements in genomes, it is interesting to speculate whether high copy number retroelements are directed to, or away from, particular areas of the genome or whether post-transposition mechanisms have operated to remove inserted elements in particular areas. [START_REF] Miyao | Target site specificity of the Tos17 retrotransposon shows a preference for insertion within genes and against insertion in retrotransposon-rich regions of the genome[END_REF] suggest that aggregation of retrotransposon insertions in host genomes is common and although mechanisms for retroelement removal were disputed for some time, these are now firmly established [START_REF] Bennetzen | Transposable elements, gene creation and genome rearrangement in flowering plants[END_REF]. In rice the estimated half-life of a retroelement is 790,000 years [START_REF] Wicker | Genome-wide comparative analysis of copia retrotransposons in Triticeae, rice and Arabidopsis reveals conserved ancient evolutionary lineages and distinct dynamics of individual copia families[END_REF]. Recombination and deletions between the long terminal repeats of the same retroelement has been proposed to occur through unequal homologous recombination [START_REF] Bennetzen | Transposable elements, gene creation and genome rearrangement in flowering plants[END_REF]. The efficiency of removal through this mechanism is expected to increase with the size of the LTR [START_REF] Bennetzen | Transposable elements, gene creation and genome rearrangement in flowering plants[END_REF] and also with an inverse of the length of the intervening sequence, and could well be the mechanism behind the removal of close facing retroelement insertions. In humans, Alu repeats constitute 10% of the genome and even though they have a different mechanism of integration to the copialike retroelements, closely spaced and inverted Alu elements are extremely unstable [START_REF] Lobachev | Inverted Alu repeats unstable in yeast are excluded from the human genome[END_REF][START_REF] Rowold | Alu elements and the human genome[END_REF][START_REF] Stenger | Biased distribution of inverted and direct Alus in the human genome: Implications for insertion, exclusion and genome stability[END_REF]. Similarly the findings of this study suggest that close proximity retroelement insertion is a similarly unstable event.

Although a growing amount of information is available on the distribution of repetitive sequences in eukaryotic genomes, this data is concentrated on a few species of high economic significance, and even with continued expansion of this activity, information on the distribution of repetitive elements in wild or uneconomic species will not become readily available. As IRAP banding numbers are dependant on retroelement copy number and genomic distribution, we have profiles and copy numbers with a random insertion model we observe that the insertion of the two high copy number elements studied here is broadly in line with the model but that close retroelement facing insertions are extremely uncommon. Although close opposed insertions may represent a small fraction of total insertion events for randomly inserting high copy number elements, this effect will be much greater for any retroelements which specifically target particular genomic locations and may be a significant mechanism for reducing their overall copy number.

Materials and methods

Plant materials:

Agave tequilana (var. Weber azul) DNA was prepared from a leaf sample by DNeasy plant mini kit (Qiagen).

Inter-Retrotransposon amplified Polymorphism (IRAP) method.

Primers were designed to complement the 5'LTR of Teq1 and Teq17 retrotransposons [START_REF] Bousios | Isolation and characterization of Ty1-copia retrotransposon sequences in the blue agave (Agave tequilana Weber var. azul) and their development as SSAP markers for phylogenetic analysis[END_REF]. In principle, retrotransposons can integrate into the genome in either orientation, giving head to head, tail to tail, or head to tail orientations. In this work single outward facing primers were used, thus only adjacent retrotransposons with 5'LTRs facing were detected (figure 1).

IRAP was performed in a 20μl reaction mixture containing 20ng DNA, 1 x PCR buffer, 5pmol primer, 200nM dNTP, 1u Taq polymerase (NEB). The PCR programme consisted of 1 cycle at 94˚C for 5 minutes followed by 30 cycles of 1 minute at 94˚C, 1 minute at 64˚C and 2 minutes at Size determination of IRAP products was achieved by excising fragments from the dried gel followed by re-amplification, subcloning into the TOPO TA vector (Invitrogen) and sequencing.

Slot Blotting

Agave tequilana genomic DNA was quantified spectrophotometrically, serially diluted, denatured in 0.2M NaOH/2M NaCl for 5 min, and transferred to a Hybond N + membrane (Amersham) using a vacuum slot blotter. Control double stranded plasmid DNAs containing the individual agave RNaseH fragments specific to each probe [START_REF] Bousios | Isolation and characterization of Ty1-copia retrotransposon sequences in the blue agave (Agave tequilana Weber var. azul) and their development as SSAP markers for phylogenetic analysis[END_REF] were quantified by ethidium bromide agarose gel electrophoresis spectrophotometrically, serially diluted and transferred as above. Cloned RNaseH PCR fragments were direct labelled with alkaline phosphatase, using AlkPhos direct labelling (Amersham). These were hybridized to the filters and washed with varying levels of stringency with the highest being equivalent to 0.1 x SSC at 60˚C. Hybridisation was detected by chemiluminescence using CDP-Star detection reagent (Amersham) using Kodak MXB film. (For each panel) Row A contains cloned target DNA with 2.207 x 10 9 molecules (left slot), 2.207

A c c e p t e d m a n u s c r i p t 12

x 10 10 molecules (middle slot) and 2.207 x 10 11 molecules (right slot). Row B contains A.

tequilana genomic DNA with 2.207 x 10 4 genomes (middle slot) and 2.207 x 10 5 genomes (right slot).

  the distribution of distances between adjacent retrotransposons, and the number of IRAP bands produced, given the number of retrotransposons and the size of the genome.

  Appendix B we present a full description of a mathematical model for the insertion of retrotransposons into the genome. The essence of the model is that retrotransposons have a uniform probability of insertion into all locations of the genome. The mathematical model plays the role of a null hypothesis in this work.

  n retrotransposons within a genome of length L, the distribution (probability density) of distances between adjacent retrotransposons is found to be very accurately given by the exponential distribution ( )

  predicted a size distribution of IRAP bands that is effectively uniform over the size range detected by IRAP (20-1500bp) -see frequency distribution U in Figure2.

  technique, in combination with an essentially null theoretical model, can be used to investigate the insertional dynamics of a particular retroelement. By comparing IRAP

  P dATP was included in the reaction to enable products to be visualized by autoradiography after resolution by denaturing acrylamide gel electrophoresis.

Figure 1 :

 1 Figure 1: Detection of close facing retroelement insertions by IRAP.

Figure 2 :

 2 Figure 2: Size distribution of IRAP products.

Figure 3 :

 3 Figure 3: Copy number of Teq1 and Teq17 Ty1-copia retrotransposons in Agave tequilana.

  Fig 1

Table 1

 1

	also contains the expected

combined with the estimate of n and the values of D and L.

A c c e p t e d m a n u s c r i p t 6 Table 1

Table 1 :

 1 This table contains the A. tequilana genome size L (in base pairs), the maximum size D (in base pairs) of the products of PCR, the copy number n, of Teq1 and Teq17 retroelements (see Appendix A for the method of copy number estimation), the expected number of IRAP bands

A c c e p t e d m a n u s c r i p t

In this Appendix, we provide details of the estimation of copy number of elements within the genome.

Estimation of total copy number was carried out by quantitative Southern blotting (see Figure 3).

Specific retroelement probes, corresponding to 164bp sections of the RNaseH gene, were generated from clones containing larger fragments of the Teq1 and Teq17 retroelements [START_REF] Bousios | Isolation and characterization of Ty1-copia retrotransposon sequences in the blue agave (Agave tequilana Weber var. azul) and their development as SSAP markers for phylogenetic analysis[END_REF]. A series of concentrations of control target DNA (row A) and A. tequilana genomic DNA (row B) were analysed. As the genome size of A. tequilana is known, namely 2C = 8.4 x 10 9 bp [START_REF] Bennett | Angiosperm DNA C-values database (release 5.0[END_REF] the exact number of genomes loaded can be calculated (1G = 1μg = 2.207 x 10 5 genomes). The amount of target was chosen to correspond to different numbers of 1G genomes, ranging from 10 4 to 10 6 . Direct comparison of the intensity of genome hybridization and that of the control sequences provide a good estimate of the number of elements in the genome. For both Teq1 and Teq17, the intensity of the genome hybridization (1G) approximately coincides with the 10 5 control lane, thereby indicating this number of element copies, per genome, in each case. This result is accurate within an overall factor of 2, hence the retrotransposon copy number has a value in the range to . 5. The fraction of the genome occupied by tips of transposable element is sufficiently small that there is negligible chance of tips of different retrotransposons inserting within one another.

Under these assumptions it is possible to predict the distribution of distances between the tips of adjacent retrotransposons and make an explicit prediction about the mean number of IRAP bands observed (see Appendix C for details).

A c c e p t e d m a n u s c r i p t 19

In particular, when there are n retrotransposons within the genome, the distribution (probability density) of distances between adjacent retrotransposons, neglecting the length of the tips, is found (see Appendix C) to be given by and when n is large (n>>1) this distribution is very accurately given by the exponential distribution

The distribution labelled U, in Figure 2, is based on Eq. ( 1).

The number of observed IRAP bands in a particular genome is a random variable that we denote In this Appendix, we give details of the calculations for the mathematical model adopted in this paper.

We assume there are n transposable elements in the genome and that these may be treated as point entities which have positions arising from a uniform probability of insertion into the genome. Thus their positions lie in the continuous range 0 to L (L = length of the genome). where the factor 1/4 follows from the assumption that all 4 four possible orientations of adjacent elements are equally likely, so, on average, only 1/4 of all adjacent retrotransposons have 5'LTRs facing.

Equation ( 5) can also be written as and using Eq. ( 4) yields which is our second theoretical result.


PhD, FRCP Kerry R Mills
email: prof.krmills@mac.com

Detecting fasciculations in amyotrophic lateral sclerosis: duration of observation required

Keywords: amyotrophic lateral sclerosis, fasciculation, EMG Word count: 1209

A practical issue in the diagnosis of amyotrophic lateral sclerosis (ALS) is how long must the EMG be observed before a muscle can be declared free of fasciculations with some degree of certainty. To answer this question, the intervals between fasciculation potentials (FPs) were recorded from 53 muscles of 19 ALS patients. The distribution of the FP intervals found across the sample showed that to record a single fasciculation with a probability approaching unity, up to 90s observation may be required.

Background

Clinical Neurophysiologists bear a heavy responsibility in the diagnosis of ALS, there being no reliable biological marker of the disease. Acute and chronic partial denervation seen on EMG provide diagnostic information of lower motor neuron abnormality and its distribution. Although they are found in many other conditions, fasciculations are the hallmark of the condition and it is difficult to be sure of the diagnosis in their absence. Both denervation and FPs may be seen in muscles without clinical abnormality emphasizing the importance of EMG [START_REF] De Carvalho | Pathophysiological significance of fasciculations in the early diagnosis of ALS[END_REF] .

FPs may be recorded by needle or surface electrodes [START_REF] Howard | Surface EMG in the recording of fasciculations[END_REF] and represent the spontaneous discharge of a motor unit or a part thereof. Needle recorded FP waveforms can be highly complex probably due to distal multifocal triggering and intermittent intramuscular axonal block 3 4 or may be simple resembling the waveform of motor units recruited by voluntary action. Their discharge rate is usually described as irregular or random [START_REF] Buchthal | Spontaneous electrical activity: an overview[END_REF] . A recent multi-electrode surface EMG study [START_REF] Kleine | Firing pattern of fasciculations in ALS: evidence for axonal and neuronal origin[END_REF] has shown that intervals may follow either a Poisson distribution as might be expected from a random rare event or a more symmetrical distribution of shorter intervals (3-15 ms), generally referred to as double FPs. The former was postulated to be related to spinal motoneuron hyperexcitability whereas the latter was thought to reflect the phase of axonal superexcitabilty 7 8 following the first discharge. Thus study of FP discharge rate may give information allowing the neural structure generating the spontaneous potential to be identified.

The electromyographer when examining a patient with potential ALS is often faced with the question: how long should I observe this needle recording before I conclude that fasciculations are absent. The question may be less critical if denervation is clearly present or the muscle is weak, but in early cases where the presentation is in a single limb or purely bulbar, then the finding of FPs in clinically uninvolved muscles assumes much greater importance. This short report attempts to answer this question by providing in a prospective study data on the probability of detecting FPs in relation to observation duration. Many electromyographers would be uncomfortable assigning significance to a single FP and would wait until more were seen to be confident of the finding.

Hence probability values have been calculated for observing 1-5 fasciculations.

Methods

FPs were recorded with a concentric needle electrode from 53 muscles (biceps: 8, tibialis anterior: 27, first dorsal interosseous: 17, trapezius: 1) in 19 patients with definite ALS according to the El Escorial criteria 9 10 modified by the Awaji consensus [START_REF] De Carvalho | Electrodiagnostic criteria for diagnosis of ALS[END_REF] . Patients were selected after EMG had been performed on the above criteria from a sequential group of 61 referred for EMG in whom ALS was a potential diagnosis. Muscles for investigation were selected on clinical grounds with the requirement that muscles from all 4 limbs and a muscle innervated by a cranial nerve were studied in all patients. An arbitrary criterion was needed as to how long to monitor a muscle for FPs before it was included in the study. If a muscle showed no FPs on EMG after 1 minute (the current clinical practice of the author), then that muscle was excluded from analysis. It is acknowledged that this may have excluded some muscles with long inter-fasciculation intervals but had the virtue of practicality. The duration of disease in the 19 patients was between 3 and 36 months (median: 18 months, interquartile range: 10.25 months). Recording at a single needle site continued until at least 50 FPs had been recorded. The EMG signal, bandpass filtered between 5 Hz and 10 kHz, was digitized at a sampling rate of 20 kHz. FPs triggered data collection using Signal 3 software (Cambridge Electronic Design, Cambridge UK). The time of occurrence of each FP from the start of the recording was logged to the nearest ms and the inter-fasciculation intervals computed. For each recording, if t 1 … t n represent the times of occurrence of FPs, the intervals (I 1 … I 5 ) were calculated from t n+1 -t n , t n+2 -t n , t n+3 -t n , t n+4 -t n and t n+5 -t n for all values of n. The maximum values of each of these series therefore represent the longest intervals in which 1, 2, 3, 4 or 5 FPs occurred. The cumulative frequency distribution of I 1 … I 5 from the whole data set was then constructed (see Fig. 1), allowing the probability of 1 … 5 fasciculations occurring with respect to the duration of the recording to be estimated. In other words, after the needle has been inserted, what range of times may elapse before 1, 2, 3, 4 or 5 fasciculations have occurred.

Results

Of the 53 muscles, 16 had normal strength and showed no evidence of acute or chronic partial denervation on qualitative EMG, 8 were severely denervated with fibrillation and marked chronic neurogenic changes. The remainder ( 29

Discussion

Secure diagnosis of ALS relies heavily on EMG. It is required to demonstrate lower motoneuron abnormalities in a distribution that cannot be explained by nerve, plexus or root lesions. In early cases where clinical abnormalities are limited to one or two regions, then EMG assumes a greater importance being able to demonstrate EMG changes in clinically unaffected muscles. In limb onset ALS, it is also particularly important to demonstrate EMG abnormality in muscles innervated by cranial nerves [START_REF] Sonoo | Utility of trapezius EMG for diagnosis of amyotrophic lateral sclerosis[END_REF] . It remains debatable whether FPs can occur in muscles showing no other EMG abnormality. Some contend that by using quantitative techniques such as macroEMG [START_REF] Guiloff | Voluntary activation and fiber density of fasciculations in motor neuron disease[END_REF] , FPs occur exclusively in partially denervated muscle. Others 1 using qualitative assessment, as has been used in the present study, find FPs in muscles with no other EMG abnormality.

Nevertheless, the confident detection of FPs in ALS is of great practical importance.

The distribution of intervals between FPs in ALS is known to be highly skewed with many short intervals and a very few long intervals [START_REF] Kleine | Firing pattern of fasciculations in ALS: evidence for axonal and neuronal origin[END_REF] . Discharge rate may be correlated to the complexity of the FP waveform [START_REF] Janko | Fasciculations in motor neuron disease: discharge rate reflects extent and recency of collateral sprouting[END_REF] , suggesting that the degree and recency of collateral reinnervation may be related to discharge rate. At a single needle site several different FPs may be seen, recognizable by the similarity of waveform components. It would be of interest to study the firing characteristics of individually identified FPs, but the purpose of the present study was a practical one: to determine the longest interval between FPs of any waveform since these will have the same significance for diagnosis.

The length of time an electromyographer should observe the EMG for fasciculation clearly depends on the clinical weighting to be put on the finding and on the experience of the electromyographer. Practicality is also important; watching EMG for FPs in a perfectly relaxed muscle for more than a few minutes would be difficult for both patient and physician. In early ALS cases where clinical abnormality is limited to one limb, the finding of FPs in the contralateral limb would have greater significance and would therefore warrant a longer search.

Any guidance therefore must take into account the clinical situation. However, from this data it may be suggested that to be certain of recording, say, 2 fasciculations then observation should be for 2 minutes before the muscle can be confidently be declared fasciculation free (Fig. 1). Similarly, if the muscle shows denervation and the electromyographer is happy to identify a single fasciculation, then the longest time needed to observe the muscle would be 90 s. Of course in most instances, FPs will be found much more quickly, the above times being the longest needed before deciding that fasciculations are almost certainly absent.

It should also be emphasized that FPs should not be taken in isolation as evidence for ALS [START_REF] Rosenfeld | Fasciculations without fibrillations: the dilemma of early diagnosis[END_REF] ; acute and/or chronic denervation in addition in other muscles will be required.

Fig. 1 .

 1 Fig.1. From this can be seen, for example, that to record 5 fasciculations with a


MD Ilse M P Arts

MD Sebastiaan Overeem

MD S Pillen

MD H Jurgen Schelhaas

MD Machiel J Zwarts
email: m.zwarts@neuro.umcn.nl

Prof J Machiel

MD Zwarts

Muscle ultrasonography to predict survival in amyotrophic lateral sclerosis

Keywords: ultrasonography, survival rate, Amyotrophic Lateral Sclerosis, prognosis

We assessed the prognostic value of muscle ultrasonography in 31 patients with amyotrophic lateral sclerosis (ALS) and compared it with accepted prognostic variables like functional capacity (measured with the ALSFRS-R rating scale) and muscle strength. Each patient was examined once. The following ultrasonography parameters were determined: muscle thickness, echo intensity (EI) and the presence of fasciculations. Correlations between baseline measurements, preslope values and survival were calculated. EI, disease duration, muscle strength preslope and ALSFRS-R preslope correlated with survival. Using a stepwise multivariate analysis, the combination of EI preslope and ALSFRS-R preslope was shown to have the best predictive value for survival.

INTRODUCTION

Survival is an important and frequently used outcome measure for clinical trials in amyotrophic lateral sclerosis (ALS). [START_REF] Gordon | Efficacy of minocycline in patients with amyotrophic lateral sclerosis: a phase III randomised trial[END_REF][START_REF] Lacomblez | Dose-ranging study of riluzole in amyotrophic lateral sclerosis. Amyotrophic Lateral Sclerosis/Riluzole Study Group II[END_REF] The median survival time from diagnosis is estimated to be around 15 months. [START_REF] Zoccolella | Analysis of survival and prognostic factors in amyotrophic lateral sclerosis: a population based study[END_REF] However, more than 25% of patients live longer than 4 years. [START_REF] Zoccolella | Analysis of survival and prognostic factors in amyotrophic lateral sclerosis: a population based study[END_REF] Due to this large interindividual variation, the use of survival as a primary outcome measure leads to studies with large sample sizes and long duration. [START_REF] Zoccolella | Analysis of survival and prognostic factors in amyotrophic lateral sclerosis: a population based study[END_REF][START_REF] De Carvalho | Can selection of rapidly progressing patients shorten clinical trials in amyotrophic lateral sclerosis?[END_REF][START_REF] Kimura | Progression rate of ALSFRS-R at time of diagnosis predicts survival time in ALS[END_REF][START_REF] Kollewe | ALSFRS-R score and its ratio: A useful predictor for ALS-progression[END_REF] One possible solution to compensate for population heterogeneity is the use of prognostic variables for stratification.

Examples of such variables are age, time from symptom onset to diagnosis, and functional capacity as measured with the revised ALS Functional Rating Scale (ALSFRS-R). [START_REF] Zoccolella | Analysis of survival and prognostic factors in amyotrophic lateral sclerosis: a population based study[END_REF][START_REF] Kimura | Progression rate of ALSFRS-R at time of diagnosis predicts survival time in ALS[END_REF][START_REF] Kollewe | ALSFRS-R score and its ratio: A useful predictor for ALS-progression[END_REF][START_REF] Gordon | Outcome measures for early phase clinical trials[END_REF] Ultrasonography, a painless and relatively simple technique, can be used to visualize structural muscle changes caused by a neuromuscular disorder. [START_REF] Pillen | Muscle ultrasound in neuromuscular disorders[END_REF] Affected muscles show a diminished muscle thickness and appear whiter (i.e. have an increased echo intensity; EI). [START_REF] Pillen | Muscle ultrasound in neuromuscular disorders[END_REF] Aging also results in diminished muscle thickness and increased EI. [START_REF] Pillen | Muscle ultrasound in neuromuscular disorders[END_REF] By quantitative analysis, with use of age and gender based normal values, it is possible to discriminate between neuromuscular and non-neuromuscular diseases with high sensitivity and specificity. [START_REF] Pillen | Muscle ultrasound in neuromuscular disorders[END_REF] Furthermore, due to its capability to register movements, ultrasonography is a highly sensitive tool to reveal fasciculations. Recently we showed that even in an early stage of ALS, atrophy, increased EI and fasciculations could be detected. [START_REF] Arts | Quantitative muscle ultrasonography in amyotrophic lateral sclerosis[END_REF] Here, we studied whether muscle ultrasonography can be used as an additional prognostic tool to predict survival in patients with ALS.

METHODS

Subjects

We included 31 patients with probable, possible or definite ALS according to the El Escorial Criteria.[10] Patients were examined once. Recordings included measurements of muscle strength, ALSFRS-R and muscle ultrasonography. Outcome was retrieved from the medical records. The survival time, defined as the time from symptom onset, was estimated using medical records. The local ethical committee approved the study and all patients gave written informed consent.

Ultrasonography

Using a standard protocol, [START_REF] Arts | Quantitative muscle ultrasonography in amyotrophic lateral sclerosis[END_REF] transverse ultrasound scans were made of the following muscles or muscle groups on both sides: sternocleidomastoid, biceps brachii including the underlying brachialis muscle, forearm flexor group, quadriceps femoris and tibialis anterior.

Measurements were performed using either a Philips IU22 with a 5-17 MHz linear broadband transducer or a HP Sonos 2000 phased-array real-time scanner with a 7.5-MHz transducer.

The system settings and conversion procedure to combine results from the two scanners have been described elsewhere. [START_REF] Pillen | Quantitative gray-scale analysis in skeletal muscle ultrasound: A comparison study of two ultrasound devices[END_REF] Thickness and EI of each muscle were determined using standardized procedures. [START_REF] Arts | Quantitative muscle ultrasonography in amyotrophic lateral sclerosis[END_REF] Furthermore, each muscle(group) was screened for fasciculations during 10 seconds. [START_REF] Arts | Quantitative muscle ultrasonography in amyotrophic lateral sclerosis[END_REF] Muscle strength and ALSFRS-R Muscle strength of 10 different muscle groups of the neck, upper and lower limbs was manually tested on both sides, and scored with the modified scale of the Medical Research Council. [START_REF]Aids to the examination of the Peripheral Nervous System[END_REF] The total sum score was calculated, with a maximum score of 100. The ALSFRS-R was used to quantify activities of daily living. [START_REF] Cedarbaum | The ALSFRS-R: a revised ALS functional rating scale that incorporates assessments of respiratory function. BDNF ALS Study Group (Phase III)[END_REF] The ALSFRS-R consists of 12 items with scores ranging from 0 to 4, with a maximum achievable score of 48 (0 = total disability, 48 = normal).

Analysis

Data were analyzed using SPSS V16.0. Data are shown as mean ± SD (range). Muscle thickness and EI were corrected for age, gender and handedness and transformed into zscores, which were then summed to yield a composite z-score (with a normal value of 0).

The prognostic value of each variable was expressed in terms of a hazard ratio using univariate analysis. For the factors reaching significance in this test, adjusted hazard ratios were estimated by multivariate analyses (forward conditional method) using a Cox proportional hazards regression model. Tested variables included age, disease duration, number of muscles with fasciculations, and baseline values of EI, muscle thickness, ALSFRS-R and strength. In addition, preslope values of EI, muscle thickness, ALSFRS-R and strength were calculated and tested as prognostic variable. Preslope values -the estimated rate of change of outcome measures prior to first examination-have shown to be good predictors of survival in various studies. [START_REF] Kimura | Progression rate of ALSFRS-R at time of diagnosis predicts survival time in ALS[END_REF][START_REF] Kollewe | ALSFRS-R score and its ratio: A useful predictor for ALS-progression[END_REF][START_REF] Czaplinski | Predictability of disease progression in amyotrophic lateral sclerosis[END_REF] Preslope values were computed as follows: (measured value -normal value) / disease duration. The relation between the different variables was assessed using Spearman's correlation coefficients.

RESULTS

A total of 31 patients were included (19 males), with an age of 60.4 years ± 12.4 (37-80) and weight of 77.1 kg ± 10.7 (53-93). Six patients were classified as possible, 22 as probable, 2 as definite and 1 patient as familial ALS. They all used riluzole. The diagnosis was established 13.8 months ± 8.6 (3-39) after symptom onset. Ultrasonography examination was performed 21.2 months ± 10.4 (5.2-51.4) after presenting symptoms. At that time the mean ALSFRS-R score was 35.4 ± 6.8 (18-46) and the mean muscle strength 84.7 ± 11.6 (57-99). At the time of analysis -approximately 6 months after inclusion of the last patient-16 out of 31 patients (51.6%) were deceased.

Patients showed a reduced muscle thickness, total z-score: -11.6 ± 13.1 (-38.8 -12.2) and increased EI: total z-score 22.0 ± 13.9 (5-68.5). The average number of muscles showing fasciculations was 6.6 (range 1-10).

Table 1 displays the results of the univariate and multivariate analysis to predict survival. As reported before, [START_REF] Zoccolella | Analysis of survival and prognostic factors in amyotrophic lateral sclerosis: a population based study[END_REF][START_REF] Kimura | Progression rate of ALSFRS-R at time of diagnosis predicts survival time in ALS[END_REF][START_REF] Kollewe | ALSFRS-R score and its ratio: A useful predictor for ALS-progression[END_REF] shorter disease duration, and lower ALSFRS-R preslope were individually associated with a shorter survival. Furthermore, lower muscle strength preslope was correlated with a worse prognosis. Of the ultrasound measures, higher EI preslope values significantly predicted shorter survival as well. to the fact that EI is the result of intramuscular fatty or fibrous tissue, an irreversible consequence of denervation. [START_REF] Pillen | Skeletal muscle ultrasound: correlation between fibrous tissue and echo intensity[END_REF][START_REF] Froes | Muscle morphometry in motor neuron disease[END_REF]. In contrast, muscle atrophy is not only caused by denervation but is also the result of other factors like disuse. Training and reinnervation can temporarily prevent muscle atrophy, but will not have considerable effects on the amount of intramuscular fatty or fibrous tissue.

The significant results in this relatively small study further show the strength of EI preslope as a predictor for survival. However, precise cut-off values need to be established in a large follow-up trial. In plot D, 33th and 67 th percentiles were used as cut off. Abbreviations: EI-p -EI preslope, ALSFRS-p -ALSFRS-R preslope.

FIGURE LEGEND

  and ALSFRS-R preslope as prognostic variables. Although being an individual predictor of survival, muscle strength preslope was not of additional value to the multivariate predictive model. This may be due to the strong correlation between muscle strength and ALSFRS-R preslope (correlation coefficient 0.74, p <0.001). The correlation between EI and ALSFRS-R preslope was relatively weak (correlation coefficient -0.41, p = 0.05), indicating that these are independent predictors. The correlation coefficients of the other preslope parameters were as follows: ALSFRS-R -muscle thickness: 0.18 (p = 0.344), muscle thickness -EI: -0.45 (p = 0.011) and EI -strength preslope: -0.51 (p = 0.004).The prognostic value of both EI and ALSFRS-R preslope is further illustrated by Kaplan-Meier plots, showing clear differences in survival when stratified using EI or ALSFRS-R preslope or the combination of both(Figure).

Kaplan-

  Meier plot describing overall cumulative survival (A) and cumulative survival after stratification by ALSFRS-R preslope (B), EI preslope (C) or the combination of EI and ALSFRS-R preslope (D). The small vertical lines represent censored cases. In plot B and C median values of ALSFRS-R and EI preslope were used as cut-off values respectively.

Table 1 . Cox proportional hazard analysis to predict survival Univariate Multivariate

 1

		Hazard ratio (95% CI) p	Hazard ratio (95% CI) p
	Age	1.02 (0.97 -1.06)	0.471
	Gender	1.27 (0.47 -3.38)	0.637
	Weight	0.98 (0.93 -1.02)	0.318
	Disease onset -diagnosis	0.93 (0.86 -1.00)	0.052
	Disease onset -ultrasound 0.94 (0.89 -0.99)	0.020
	Number of muscles with	1.09 (0.89 -1.33)	0.388
	fasciculations
	Total muscle thickness	1.00 (0.96 -1.06)	0.830
	Total EI	1.01 (0.98 -1.05)	0.574
	Total strength	0.98 (0.94 -1.01)	0.171
	Total ALSFRS-R	0.96 (0.90 -1.02)	0.162
	EI preslope	4.72 (2.02 -11.03)	0.001	4.88 (1.06 -22.44)	0.042
	Muscle thickness preslope 0.89 (0.38 -2.07)	0.774
	Strength preslope	0.08 (0.03 -0.25)	< 0.001
	ALSFRS-R preslope	0.021 (0.003 -0.15)	< 0.001 0.019 (0.002 -0.23)	0.002
	Stepwise multivariate analysis revealed a predictive model that only included the EI preslope

Acknowledgements:

We thank Mirjam Heykers, Henny Janssen and Wilma Raijmann, for performing the ultrasonography measurements.

DISCUSSION

Our study shows that an increased muscle EI preslope predicts survival in patients with ALS.

Moreover, increased EI had additional prognostic value when combined with the ALSFRS-R preslope, an already accepted predictor for survival. [START_REF] Gordon | Efficacy of minocycline in patients with amyotrophic lateral sclerosis: a phase III randomised trial[END_REF][START_REF] Kimura | Progression rate of ALSFRS-R at time of diagnosis predicts survival time in ALS[END_REF][START_REF] Kollewe | ALSFRS-R score and its ratio: A useful predictor for ALS-progression[END_REF] The independent value of both measures is supported by their rather low correlation, suggesting they represent different aspects of the disease process in ALS. The number of fasciculations did not predict survival.

Although muscle strength preslope was strongly predictive for survival in the univariate analysis is, it did not emerge as an additional predictive parameter in the stepwise multivariate analysis. Most likely, this was due to the strong correlation between muscle strength and the ALSFRS-R preslope, so adding muscle strength to the model did not increase its predictive value.

Preslope values were better predictors of survival than baseline values or disease duration.

Preslope values combine these two aspects, probably making these a better indicator of the rate of "preclinical" disease progression. A practical advantage of preslope values is the fact that they are based on just one evaluation point that can be routinely acquired during the diagnostic process. In addition, preslope values are convenient in clinical trials, since these do not require a lead-in period. [START_REF] De Carvalho | Can selection of rapidly progressing patients shorten clinical trials in amyotrophic lateral sclerosis?[END_REF] However, it would be useful to directly compare preslope ultrasonography measures with a longitudinal repeated measures design as a final 'validation'.

Competing interests, Funding:

None of the authors has financial or other conflicts of interest

Licence for Publication:

The Corresponding Author has the right to grant on behalf of all authors and does grant on behalf of all authors, an exclusive licence on a worldwide basis to the BMJ Publishing Group Ltd to permit this article (if accepted) to be published in JNNP and any other BMJPGL products and sublicences such use and exploit all subsidiary rights, as set out in our licence.