Marc Arnaudon

Frank Nielsen
email: frank.nielsen@acm.org

MEDIANS AND MEANS IN FINSLER GEOMETRY

Medians and means in Finsler geometry

Introduction

The geometric barycenter of a set of points is the point which minimizes the sum of the squared distances to these points. It is the most traditional estimator is statistics that is however sensitive to outliers [START_REF] Hampel | Robust Statistics The Approach Based on Influence Function[END_REF]. Thus it is natural to replace the average distance squaring (power 2) by taking the power of p for some p ∈ [START_REF] Afsari | Riemannian L p center of mass : existence, uniqueness, and convexity[END_REF][START_REF] Amary | Methods of Information Geometry Translations of Mathematical Monographs[END_REF]. This leads to the definition of p-means. When p = 1, the minimizer is the median of the set of points, very often used in robust statistics [START_REF] Hampel | Robust Statistics The Approach Based on Influence Function[END_REF]. In many applications, pmeans with some p ∈ (1, 2) give the best compromise. For existence and uniqueness in Riemannian manifolds under convexity conditions on the support of the measure, see Afsari [START_REF] Afsari | Riemannian L p center of mass : existence, uniqueness, and convexity[END_REF].

The Fermat-Weber problem concerns finding the median e 1 of a set of points in an Euclidean space. Numerous authors worked out algorithms for computing e 1 . The first algorithm was proposed by Weiszfeld in [START_REF] Weiszfeld | Sur le point pour lequel la somme des distances de n points donnés est minimum[END_REF] (see also [START_REF] Vardi | The multivariate L 1 -median and associated data dept[END_REF]). It has been extended to sufficiently small domains in Riemannian manifolds with nonnegative curvature by Fletcher and al. in [START_REF] Fletcher | The geometric median on Riemannian manifolds with application to robust atlas estimation[END_REF]. A complete generalization to manifolds with positive or negative curvature (under some convexity conditions in positive curvature), has been recently given by Yang in [START_REF] Yang | Riemannian median and its estimation[END_REF].

The Riemannian barycenter or Karcher mean of a set of points in a manifold or more generally of a probability measure has been extensively studied, see e.g. [START_REF] Karcher | Riemannian center of mass and mollifier smoothing[END_REF], [START_REF] Kendall | Probability, convexity and harmonic maps with small image I: uniqueness and fine existence[END_REF], [START_REF] Kendall | Convexity and the hemisphere[END_REF], [START_REF] Emery | Sur le barycentre d'une probabilité dans une variété[END_REF], [START_REF] Picard | Barycentres et martingales dans les variétés[END_REF], [START_REF] Arnaudon | Barycenters of measures transported by stochastic flows[END_REF], [START_REF] Corcuera | Riemannian barycentres and geodesic convexity Math[END_REF], where questions of existence, uniqueness, stability, relation with martingales in manifolds, behavior when measures are pushed by stochastic flows have been considered. The Riemannian barycenter corresponds to p = 2 in the above description. Computation of Riemannian barycenters by gradient descent has been performed by Le in [START_REF] Le | Estimation of Riemannian barycentres[END_REF].

The aim of this paper is to extend to the context of Finsler manifolds the results on existence and uniqueness of p-means of probability measures, as well as algorithms for computing them. Some convexity is needed, and as we shall see the fact that comparison results for triangles as Alexandroff and Toponogov theorems do not exist impose more restrictions on the support of the probability measure. As a consequence, the sharp results on existence and uniqueness established by Afsari [START_REF] Afsari | Riemannian L p center of mass : existence, uniqueness, and convexity[END_REF] and the algorithm for computing means of Yang in [START_REF] Yang | Riemannian median and its estimation[END_REF] do not extend to Finsler manifolds.

The motivation for this work primarily comes from signal filtering and denoising in the context of Diffusion Tensor Imaging (DTI), High Angular Resolution Imaging (HARDI, see [START_REF] Shen | Riemann-Finsler geometry, with applications to information geometry[END_REF], [START_REF] Pchaud | Brain connectivity using geodesics in hardi[END_REF], [START_REF] Astola | Finsler geometry on higher oredre tensor fields and applications to high angular resolution diffusion imaging[END_REF]), Orientation Distribution Function (ODF), active contours [START_REF] Melonakos | Finsler active contours[END_REF]. Applications with experimental results of an implementation will be reported in forthcoming papers.

Information geometry at its heart considers the differential geometry nature of probability distributions induced by a divergence function. In probability theory, invariance by monotonic re-parameterization and sufficient statistics yields the class of f -divergences [START_REF] Amary | Methods of Information Geometry Translations of Mathematical Monographs[END_REF]

I f (p, q) = p(x)f ( q(x) p(x)
)dx that includes the Kullback-Leibler (KL) information-theoretic divergence KL(p, q) = p(x) log p(x) q(x) dx as its prominent member (for f (t) = -log t). It is well-known that the KL divergence (better known as the relative entropy) yields a dually flat structure [START_REF] Amary | Methods of Information Geometry Translations of Mathematical Monographs[END_REF] generalizing the (self-dual) Euclidean space.

Because divergences are usually asymmetric and violate the triangle inequality they have not been extensively considered from an algorithmic point of view. Indeed, the triangle inequality property is often used in computational geometry to design efficient algorithms by allowing various "pruning" techniques [START_REF] Elkan | Using the Triangle Inequality to Accelerate k-Means[END_REF][START_REF] Fuhry | Efficient skyline computation in metric space[END_REF]. Computational geometry has thus mostly considered metric spaces for keeping the triangle inequality properties.

One can metrize divergences. The KL divergence can be symmetrized either into the Jeffreys divergence J(p, q) = KL(p, q) + KL(q, p) or the Jensen-Shannon (JS) divergence:

JS(p, q) = KL(p, p + q 2 ) + KL(q, p + q 2 ) = (p(x) log 2p(x) p(x) + q(x) + q(x) log 2q(x) p(x) + q(x) )dx.
The latter is preferred in practice because it is bounded, and its square root yields a metric that can be embedded into a Hilbert space [START_REF] Fuglede | Jensen-Shannon divergence and Hilbert space embedding[END_REF]. Finsler distances, arising from the underlying the Finsler metrics, are attractive as they preserve the triangle inequality [START_REF] Shen | Riemann-Finsler geometry, with applications to information geometry[END_REF] for efficient algorithmics but potentially model asymmetric distances.

In information geometry, the regular divergence D associated to a Finslerian metric distance d can be defined as D(p, q) = d 2 (p, q). Observe that the Finslerianbased divergence looses then the triangle inequality property [START_REF] Shen | Riemann-Finsler geometry, with applications to information geometry[END_REF]. (eg., the squared Euclidean distance does not satisfy the triangle inequality).

Preliminaries

Let M be a smooth manifold. On M we consider a Finsler structure F :

T M → R + . For any x ∈ M , V, X, Y, Z ∈ T x M such that V = 0, let (2.1) g V (X, Y ) := 1 2 ∂ 2 ∂s∂t (s,t)=(0,0) F 2 (V + sX + tY ).
(we shall also use the notation <X, Y > V = g V (X, Y )) and

(2.2) <X, Y, Z> V := 1 4

∂ 3 ∂r∂s∂t (r,s,t)=(0,0,0) F 2 (V + rX + sY + tZ).
We have

(2.3) <X, Y, Z> V = 1 2 ∂ ∂r r=0 g V +rX (Y, Z)
and in particular since

F 2 is 2-homogeneous and V → g V (X, Y ) is 0-homogeneous, (2.4) <V, Y, Z> V = 0.
Let V be a non-vanishing vector field on M . The Chern connection ∇ V is torsionfree and almost metric, and can be characterized by

(2.5) X<Y, Z> V = <∇ V X Y, Z> V + <Y, ∇ V X Z> V + 2<∇ V X V, Y, Z> V .
More precisely, parameterizing locally T M by coordinates (x 1 , . . . , x m , y 1 = dx 1 , . . . , y m = dx m ), defining the geodesic coefficients as

(2.6) G i (y) = 1 4 g ik (y) 2 ∂g jk ∂x l - ∂g jl ∂x k y j y l , y ∈ T M \{0}, letting (2.7) N i j = ∂G i ∂y j , δ δx i = ∂ ∂x i -N k i (y)
∂ ∂y k ∈ T y (T M \{0}) , then the Christoffel symbols of the Chern connection are given by

(2.8) Γ k ij = 1 2 g kl δg lj δx i + δg il δx j - δg ij δx l
(see [START_REF] Chern | Riemann-Finsler geometry[END_REF]). Note that defining (2.9) δy i = dy i + N i j (y)dx j we have for a smooth function f : T M \{0} → R (2.10) df = δf δx i dx i + ∂f ∂y i δy i . The Chern curvature tensor is defined by the equation

(2.11) R V (X, Y )Z := ∇ V X ∇ V Y Z -∇ V Y ∇ V X Z -∇ V [X,Y ] Z, and the flag curvature is (2.12) K (V, W ) := <R V (V, W )W, V > <V, V > V <W, W > V -<V, W > 2 V , for two non collinear V, W ∈ T x M .
We say that M has nonpositive flag curvature if for all V, W , K (V, W ) ≤ 0.

The tangent curvature of two vectors V, W ∈ T x M is defined as

(2.13) T V (W ) = <∇ W W W -∇ V W W , V > V
where W is a vector field satisfying Wx = W . For a nonnegative constant δ ≥ 0 we say that T ≥ -δ or T ≤ δ if respectively

(2.14) T V (W ) ≥ -δF (V )F (W ) 2 or T V (W ) ≤ δF (V )F (W ) 2 .
For x ∈ M we define

(2.15) C (x) = sup v,w∈TxM\{0} <v, v> v <v, v> w , D(x) = sup v,w∈TxM\{0}
<v, v> w <v, v> v .

Remark 2.1. For applications in active contours, a "Wulff shape" is given which does not depend on x and defines the Finsler structure. From this shape C and D can easily be calculated. See e.g [START_REF] Melonakos | Finsler active contours[END_REF] and [START_REF] Zach | Globally Optimal Finsler Active Contours[END_REF].

A geodesic in M is a curve t → c(t) satisfying for all t, ∇ ċ(t) ċ(t) ċ = 0. It is well known that a geodesic has constant speed, and that it locally minimizes the distance ( [START_REF] Chern | Riemann-Finsler geometry[END_REF]). If so, letting ρ(x, y) the forward distance from x to y, then (2.16)

ρ 2 (x, y) = < ċ(0), ċ(0)> ċ(0)
where t → c(t) is the minimal geodesic satisfying c(0) = x and c(1) = y. By definition, the backward distance from x to y is ρ(y, x).

For x ∈ M and v ∈ T x M , we let whenever it exists exp x (v) := c(1) where t → c(t) is the geodesic satisfying ċ(0) = v.
If M is complete, analytic, simply connected and has nonpositive flag curvature (we say that M is an analytic Cartan-Hadamard manifold), then exp x : T x M → M is an homeomorphism ([6] theorem 4.7). Under these assumption, letting for x, y ∈ M , -→ xy = exp -1 x (y), we have (2.17)

ρ 2 (x, y) = < -→ xy, -→ xy>-→ xy .
For x 0 ∈ M and R > 0,, let us denote by B(x 0 , R) (resp. B(x 0 , R)) the (forward) open (resp. closed) ball with center x 0 and radius R:

(2.18) B(x 0 , R) = {y ∈ M, ρ(x 0 , y) < R} (resp. B(x 0 , R) = {y ∈ M, ρ(x 0 , y) ≤ R}) Now let (t, s) → c(t, s) a family of minimizing geodesics t → c(t, s), t ∈ [0, 1], parametrized by s ∈ I, I an interval in R. Define (2.19) E(s) = 1 2 ρ 2 (c(0, s), c(1, s)).
The computation of E ′ (s) and E ′′ (s) is well-known, see e.g. [START_REF] Bao | An introduction to Riemann-Finsler geometry[END_REF]. We have

(2.20) E ′ (s) = <∂ s c(1, s), ∂ t c(1, s)> ∂tc(1,s) -<∂ s c(0, s), ∂ t c(0, s)> ∂tc(0,s) .
As for the second derivative, letting c = c(•, 0), and for X, Y vector fields along c

(2.21) I(X, Y ) = 1 0 <∇ T T X, ∇ T T Y > T -<R T (X, T )T, Y > T dt
the index of X and Y , we have

E ′′ (0) = <∇ ∂tc(1,0) ∂sc(1,0) ∂ s c(1, •), ∂ t c(1, 0)> ∂tc(1,0) -<∇ ∂tc(0,0) ∂sc(0,0) ∂ s c(0, •), ∂ t c(0, 0)> ∂tc(0,0) + I(∂ s c(•, 0), ∂ s c(•, 0)). (2.22)
Assuming s → c(0, s) and s → c(1, s) are geodesics, we obtain

E ′′ (0) = T ∂tc(0,0) (∂ s c(0, 0)) -T ∂tc(1,0) (∂ s c(1, 0)) + I(∂ s c(•, 0), ∂ s c(•, 0)). (2.23)
We are interested in the situation where c(1, s) ≡ z a constant. In this case we have

E ′′ (0) = T ∂tc(0,0) (∂ s c(0, 0)) + I(∂ s c(•, 0), ∂ s c(•, 0)). (2.24) For p ≥ 1, define (2.25) D p (s) = ρ p (c(0, s), z) Proposition 2.2. Assume K ≤ k, T ≥ -δ, C ≤ C for some k, δ ≥ 0, C ≥ 1. Let p > 1.
Then writing r = ρ(x, z),

(2.26) D ′′ p (0) ≥ pr p-2 min p -1, √ kr cos( √ kr) sin( √ kr) C -2 -δr . If z and x = c(0, 0) satisfy ρ(x, z) < R(p, k, δ, C) with (2.27) R(p, k, δ, C) = min p -1 C 2 δ , 1 √ k arctan √ k C 2 δ
and the injectivity radius at x is strictly larger than R(p, k, δ, C), then

D ′′ p (0) > 0. Remark 2.3. Note if p ≥ 2 then R(p, k, δ, C) = R(2, k, δ, C) = 1 √ k arctan √ k C 2 δ . Proof. Define T (t) = ∂ t c(t, 0), J(t) = ∂ s c(t, 0), J T (t) = 1 F (T (t)) 2 <J(t), T (t)> T (t) T (t), J N (t) = J(t) -J T (t).
Using successively [START_REF] Bao | An introduction to Riemann-Finsler geometry[END_REF] Lemma 9.5.1 which compares the index I(J, J) with the one of its "transplant" into a manifold with constant curvature k 2 and the index lemma [START_REF] Bao | An introduction to Riemann-Finsler geometry[END_REF] Lemma 7.3.2 which compares the index of the transplant to the one of the Jacobi field with same boundary values, we get, letting r = ρ(x, z) = D 1 (0), (2.28)

I(J, J) ≥ √ kr cos( √ kr) sin( √ kr) <J N (0), J N (0)> T (0) + <J T (0), J T (0)> T (0) .
Using the expression (2.24) for E ′′ (0) we obtain

(2.29) E ′′ (0) ≥ -δr + √ kr cos( √ kr) sin( √ kr) <J N (0), J N (0)> T (0) + <J T (0), J T (0)> T (0) .
We have from (2.20)

(2.30) E ′ (0) 2 = r 2 <J T (0), J T (0)> T (0) . Now from D 1 (s) = 2E(s) we get (2.31) D ′ 1 (s) = E ′ (s) D 1 (s) , D ′′ 1 (s) = E ′′ (s) D 1 (s) - E ′ (s) 2 D 3 1 (s)
, and this yields

D ′′ p (0) = pD 1 (0) p-2 (p -1)D ′ 1 (0) 2 + D 1 (0)D ′′ 1 (0) = pr p-2 (p -2)<J T (0), J T (0)> T (0) + E ′′ (0) ≥ pr p-2 (p -1)<J T (0), J T (0)> T (0) -δr + √ kr cos( √ kr) sin( √ kr) <J N (0), J N (0)> T (0) ≥ pr p-2 min p -1, √ kr cos( √ kr) sin( √ kr) <J(0), J(0)> T (0) -δr ≥ pr p-2 min p -1, √ kr cos( √ kr) sin( √ kr) C -2 -δr .
¿From this bound the rest of the proof follows easily.

Similarly, we can obtain an upper bound for D ′′ p (0): Proposition 2.4. Assume the sectional curvatures K have a lower bound -β 2 for some β > 0, and T ≤ δ ′ for some δ ′ > 0, D ≤ D for some D ≥ 1. Again let r = ρ(x, z), assume that the injectivity radius at x is larger than r. Then

(2.32) D ′′ p (0) ≤ pr p-2 D 2 max (p -1, βr coth(βr)) + δ ′ r . Proof.
We have by (2.24) and (2.31) together with the fact that (2.33)

I(J, J) = <J T (0), J T (0)> + I(J N , J N ), D ′′ p (0) = pr p-2 p -1)<J T (0), J T (0)> T (0) + I(J N , J N ) + T T (J)
. Let t → X(t) the parallel vector field along t → c(t, 0) with initial condition J N (0), and for t ∈ [0, 1], let

G(t) = cosh(rβt) -coth (rβ) sinh(rβt).
This is the solution of G ′′ = rβG with conditions G(0) = 1 and G(1) = 0. The vector field t → Y (t) along t → c(t, 0) defined by (2.34)

Y (t) = G(t)X(t)
has same boundary values as t → J N (t), so by the index lemma [START_REF] Bao | An introduction to Riemann-Finsler geometry[END_REF]

I(J N , J N ) ≤ I(Y, Y ).
On the other hand

I(Y, Y ) = 1 0 G ′ (t) 2 <J N (0), J N (0)> T (0) -G(t) 2 <R T (X(t), T (t))T (t), X(t)> T (t) dt ≤ <J N (0), J N (0)> T (0) 1 0 G ′ (t) 2 + r 2 β 2 G(t) 2 dt = <J N (0), J N (0)> T (0) [G ′ (t)G(t)] 1 0 + 1 0 G(t) -G ′′ (t) 2 + r 2 β 2 G(t) dt = <J N (0), J N (0)> T (0) rβ coth (rβ) . So D ′′ p (0) ≤ pr p-2 (p -1)<J T (0), J T (0)> T (0) + rβ coth(rβ)<J N (0), J N (0)> T (0) + δ ′ r ≤ pr p-2 max ((p -1), rβ coth(rβ)) <J(0), J(0)> T (0) + δ ′ r ≤ pr p-2 D 2 max ((p -1), rβ coth(rβ)) + δ ′ r (2.36) since F (J(0)) = 1.
For x ∈ M , let ℓ x : T x M → T * x M be the Legendre transformation, defined as (2.37)

ℓ x (V ) = g V (V, •) if V = 0, ℓ x (0) = 0.
It is well-known that ℓ x is a bijection. The global Legendre transformation on T M is defined as

(2.38) L (V ) = ℓ π(V ) (V )
where π : T M → M is the canonical projection. If we define the dual Minkowski norm

F * on T * x M as (2.39) F * (ξ) = max{ξ(y), y ∈ T x M, F (y) = 1}, then (2.40) F = F * • L and for non zero V ∈ T M and α ∈ T * M , (2.41) L (V ), V = F (V ) 2 , α, L -1 (α) = F * (α) 2
(see e.g. [START_REF] Alvarez Paiva | Some problems on Finsler geometry[END_REF]).

For f a C 1 function on M we may define the gradient of f

(2.42) grad f = L -1 (df ).

Forward p-means

Let µ be a compactly supported probability measure in M . For p > 1 and x ∈ M we define

(3.1) E µ,p (x) = M ρ p (x, z) µ(dz).
The (forward) p-mean of µ is the point e p of M where E µ,p reaches its minimum whenever it exists and is unique.

In this paper we will consider forward p-means and we will call them p-means. Similarly we could define the backward p-mean ←e p as the point which minimizes

x → ← - E µ,p (x) := M ρ p (z, x) µ(dz).
Depending on the context, forward or backward mean is more appropriate. One should note that defining the reverse (or adjoint) Finsler structure ← -

F (v) = F (-v), v ∈ T M ,
it is easy to check that the associated distance ←ρ satisfies ←ρ (z, x) = ρ(x, z), and forward p-mean for ← -F is backward p-mean for F . So without loss of generality we can consider only the forward p-means.

One should also note that in High Angular Resolution Imaging the Finsler structure is symmetric, so both notions coincide. It is not the case for the application concerning active contours where it is natural to consider non symmetric F .

Even if it is in a non-Finslerian context, one can give the example of rightsided and left-sided Kullback-Leibler divergences for families of Gaussian probability densities (see [START_REF] Nielsen | Sided and symmetrized Bregman centroids[END_REF]). The left sided centroid focuses on the highest mode (it is zero-forcing), and the right-sided centroid tries to cover the support of both normals (it is zero-avoiding as depicted in Fig. 2 of [START_REF] Nielsen | Sided and symmetrized Bregman centroids[END_REF]). Proposition 3.1. Assume there exists C > 0 such that C (x) ≤ C for all x ∈ M , where C (x) is defined in (2.15). Assume furthermore that supp(µ) ⊂ B(x 0 , R) for some x 0 ∈ M and R > 0. Then x → E µ,p (x) has at least one global minimum in B(x 0 , C(1 + C)R).

Proof. We begin with establishing that for all y 1 , y 2 ∈ M ,

(3.2) 1 C ρ(y 2 , y 1 ) ≤ ρ(y 1 , y 2 ) ≤ Cρ(y 2 , y 1 ).
It is sufficient to establish the second inequality and then to exchange y 1 and y 2 . If t → ϕ(t) is a path from y 1 = ϕ(0) and y 2 = ϕ(1) then its length L(ϕ) satisfies

L(ϕ) = 1 0 < φ(t), φ(t)> φ(t) dt = 1 0 < -φ(t), -φ(t)> φ(t) < -φ(t), -φ(t)> -φ(t) < -φ(t), -φ(t)> -φ(t) dt ≤ 1 0 C (ϕ(t)) < -φ(t), -φ(t)> -φ(t) dt ≤ C 1 0 < -φ(t), -φ(t)> -φ(t) dt = CL( φ)
where φ is the path from y 2 to y 1 defined by φ(t) = ϕ(1t). Minimizing over all paths φ from y 2 to y 1 we get

(3.3) ρ(y 1 , y 2 ) ≤ Cρ(y 2 , y 1 ). Now if supp(µ) ⊂ B(x 0 , R) then E µ,p (x 0 ) ≤ R p . On the other hand, if x ∈ B(x 0 , C(1 + C)R) then for all y ∈ B(x 0 , R) ρ(x, y) ≥ ρ(x, x 0 ) -ρ(y, x 0 ) ≥ 1 C ρ(x 0 , x) -Cρ(x 0 , y) ≥ (1 + C)R -CR = R
and this clearly implies that E µ,p (x) ≥ R p . From this we get the conclusion.

Concerning the uniqueness of the global minimum of E µ,p , we also have the following easy result. Proposition 3.2. Assume that µ is supported by a compact ball B(x 0 , R), and that for all z ∈ B(x 0 , R), the function x → ρ p (x, z) is strictly convex in B(x 0 , C(1 + C)R). Then µ has a unique p-mean in B(x 0 , C(1 + C)R).

Proof. If x → ρ p (x, z) is strictly convex for all z in the support of µ then E µ,p is strictly convex, and this implies that it has a unique minimum, which is attained at a unique point e p .

Corollary 3.3. Assume K ≤ k, T ≥ -δ, C ≤ C for some k, δ ≥ 0, C ≥ 1. Let p > 1. Again let R(p, k, δ, C) = min p -1 C 2 δ , 1 √ k arctan √ k C 2 δ
If µ is supported by a geodesic ball B(x 0 , R) with

(3.4) R ≤ 1 C(C + 1) 2 R(p, k, δ, C)
and the injectivity radius at any x ∈ B(x 0 , C(1 + C)R) is strictly larger than R(p, k, δ, C) then µ has a unique p-mean e p satisfying

(3.5) e p ∈ B x 0 , 1 C + 1 R(p, k, δ, C) . Proof. If x, z ∈ B(x 0 , C(1 + C)R) then ρ(x, z) ≤ ρ(x, x 0 ) + ρ(x 0 , z) ≤ (1 + C) 2 CR ≤ R(p, k, δ, C).
Using proposition 2.2, we obtain that E µ,p is strictly convex on B(x 0 , C(1 + C)R).

So by proposition 3.2 µ has a unique p-mean in B(x 0 , C(1 + C)R).

Remark 3.4. Letting x 0 ∈ M , D be a relatively compact neighborhood of x 0 , then K and C are bounded above on D by, say k D and C D , and T is bounded below on D by -δ D . Using these bounds instead of k, C and δ, we can find R sufficiently small so that the conditions of corollary 3.3 are fulfilled. So we can say any measure µ with sufficiently small support has a unique p-mean.

Remark 3.5. If M is a Cartan-Hadamard manifold, we recover the fact that we can take R(p, k, δ, C) as large as we want. More generally, in the Riemannian case, Afsari [START_REF] Afsari | Riemannian L p center of mass : existence, uniqueness, and convexity[END_REF] proved existence and uniqueness of p-means, p ≥ 1 on geodesic balls with radius r <

1 2 min inj(M ), π 2 √ k if p ∈ [1, 2), and r < 1 2 min inj(M ), π √ k if p ≥ 2.
Even taking δ = 0 and C = 1 in Corollary 3.3 the support of µ has half the size of the one in [START_REF] Afsari | Riemannian L p center of mass : existence, uniqueness, and convexity[END_REF] for p ∈ (1, 2) due to the fact that we have an additional condition (3.4) coming from the non optimality of Proposition 3.2 in the Riemannian context. As for p ≥ 2 another factor two is gained in [START_REF] Afsari | Riemannian L p center of mass : existence, uniqueness, and convexity[END_REF] with repeated use of Toponogov and Alexandroff theorems which are not available in our context.

Remark 3.6. The condition on injectivity radius is the same as in the Riemannian case. The cut locus of any point of x ∈ B(x 0 , C(1+C)R) has to be at distance larger than R(p, k, δ, C). As for Riemannian manifold there is no general condition which insures this property for cut points, but for conjugate points the same condition holds, due to Rauch comparison theorem, see Theorem 9.6.1 in [START_REF] Bao | An introduction to Riemann-Finsler geometry[END_REF]. In the particular case when M is a Cartan-Hadamard Finsler manifold, i.e. it has nonpositive flag curvature and it is simply connected, then the injectivity radius in everywhere infinite (see Theorem 9.4.1 in [START_REF] Bao | An introduction to Riemann-Finsler geometry[END_REF]).

Proposition 3.7. Let a → x(a) solve the equation

(3.6) x(0) = x 0 and for a ≥ 0 x ′ (a) = grad x(a) (-E µ,p ).
Under the conditions of Corollary 3.3, the path a → x(a) converges as a → ∞ to the p-mean of µ.

Proof. If f (a) = (-E µ )(x(a)) we have as soon as grad x(a) (-E µ,p ) = 0,

f ′ (a) = d x(a) (-E µ,p ), x ′ (a) = d x(a) (-E µ,p ), grad x(a) (-E µ,p ) = d x(a) (-E µ,p ), L -1 (d x(a) (-E µ,p )) = F * d x(a) (-E µ,p )
2 by (2.40) and (2.41).

On the other hand, we have f (0) ≥ -R p and f is nondecreasing. This implies that for all a ≥ 0, x(a) ∈ B(x 0 , C(1 + C)R), since for all x ∈ B(x 0 , C(1 + C)R), E µ,p (x) ≥ R p . As a consequence x(a) has limit points as a goes to infinity, and since f (a) converges, any limit point is a critical point of x → E µ,p (x). But by Proposition 3.2 E µ,p has a unique critical point in B(x 0 , C(1 + C)R) which is the mean e p of µ. So we can conclude that x(a) converges to e p .

Forward median

Let µ be a compactly supported probability measure in M . For x ∈ M we define (4.1)

F µ (x) = M ρ(x, z) µ(dz).
The median of µ is the point in M where F µ reaches its minimum whenever it exists and is unique. Again we have the following result.

Proposition 4.1. Assume there exists C > 0 such that C (x) ≤ C for all x ∈ M . Assume furthermore that supp(µ) ⊂ B(x 0 , R) for some x 0 ∈ M and R > 0. Then x → F µ (x) has at least one global minimum in B(x 0 , C(1 + C)R).

Proposition 4.2. Assume that µ is supported by a compact ball B(x 0 , R), that the support of µ is not contained in a single geodesic and that for all z ∈ B(x 0 , R), the forward distance to z is convex, and strictly convex in any geodesic of B(x 0 , C(1 + C)R) which does not contain z. Then µ has a unique median m ∈ B(x 0 , C(1+C)R).

Proof. Clearly under these assumptions F µ is strictly convex, so it has a unique local minimum, this minimum is global and is attained at a unique point m ∈ B(x 0 , C(1 + C)R).

Remark 4.3. Contrarily to the case of p-means for p > 1, we cannot say at this stage that any probability measure µ with sufficiently small support has a unique median, since we don't know whether F µ is strictly convex or not. In the next proposition we give a sufficient condition for strict convexity of F µ .

Proposition 4.4. Assume K ≤ k and T ≥ -δ for some k, δ > 0. Assume that the injectivity radius at any point of B(x 0 , C(1

+ C)R) is larger than (C 2 + C + 1)R. Define η = min M √ k cot √ kρ(π(v), z) <v N , v N >---→ π(v)z µ(dz), v ∈ T M satisfying π(v) ∈ B(x 0 , C(1 + C)R), F (v) = 1 (4.2)
where v N is the normal part of v with respect to the vector ---→ π(v)z and the scalar product

<•, •>---→ π(v)z . If η -δ > 0 then F µ is strictly convex on B(x 0 , C(1 + C)R).
More precisely, for all x ∈ B(x 0 , C(1 + C)R) and for all unit speed geodesic γ starting at x,

(4.3) (F µ • γ) ′′ (0) ≥ η -δ.
Proof. With the notations of section 2, from (2.31) we have (4.4)

D ′′ 1 (0) = r -1 E ′′ (0) -<J T (0), J T (0)> T (0)
. Let γ(s) = c(0, s), the unit speed geodesic with initial condition v = J(0), and f (s) = F µ (γ(s)). Equation (4.4) together with (2.29) gives (4.5)

f ′′ (0) ≥ -δ + M √ k cot √ kρ(π(v), z) <v N , v N >---→ π(v)z µ(dz) ≥ η -δ.
From this we get the condition for the strict convexity of F µ .

For x ∈ M define the measure µ x = µµ({x})δ x . Then the map y → F µx (y) is differentiable at y = x. Since and for v ∈ T x M , F µ is differentiable in the direction v with derivative (4.7)

dF µ , v = dF µx , v + µ({x})F (-v),
we see that x is a local minimum of F µ if and only if for all nonzero v ∈ T x M

(4.8) µ({x})F (-v) ≥ dF µx , -v which is equivalent to (4.9) µ({x}) ≥ (F * (dF µx )) 2 F (L -1 (dF µx )) (take -v = L -1 (dF µx ) F (L -1 (dF µx ))
). But since Alternatively, (4.12)

F * = F • L -1 ,
H(x) = L -1 M\{x} L - 1 ρ(x, z) -→ xz µ(dz) .
Let a → x(a) be the path in M defined by x(0) = x 0 and

(4.13) ẋ(a) = -H(x(a)) if for all a ′ ≤ a, µ({x(a ′ )}) < F * dF µ x(a ′ ) ; ẋ(a) = 0 if for some a ′ ≤ a, µ({x(a ′ )}) ≥ F * dF µ x(a ′ ) . Define (4.14) f (a) = F µ (x(a)).
We have for the right derivative of f when x(a) is not a minimal point of F µ :

f ′ + (a) = d x(a) (F µ x(a) ), ẋ(a) + µ({x(a)})F (-ẋ(a)) = -d x(a) (F µ x(a) ), L -1 (d x(a) (F µ x(a) )) + µ({x(a)})F L -1 (d x(a) (F µ x(a) )) = -F * d x(a) (F µ x(a) ) 2 + µ({x(a)})F L -1 (d x(a) (F µ x(a) )) .
We get

(4.15) f ′ + (a) = -F * d x(a) (F µ x(a) ) F * d x(a) (F µ x(a)
)µ({x(a)}) which is negative as soon as x(a) is not a minimal point of F µ . From this we get the following Proposition 4.6. Assume that µ is supported by a compact ball B(x 0 , R), that the support of µ is not contained in a single geodesic and that for all z ∈ B(x 0 , R), the forward distance to z is convex, and strictly convex in any geodesic of B(x 0 , C(1 + C)R) which does not contain z. Then the path a → x(a) converges to the median m of µ.

Proof. Similar to the proof of proposition 3.7

5. An algorithm for computing p-means

Lemma 5.1. Assume K ≥ -β 2 , T ≤ δ ′ , D ≤ D with β > 0, δ ′ ≥ 0 D ≥ 1. For p > 1, r > 0, define (5.1) H(r) = H p,β,D,δ ′ (r) := pr p-2 D 2 max ((p -1), rβ coth(rβ)) + δ ′ r .
If µ is a probability measure on M with bounded support and x ∈ M , define

(5.2) H µ (x) = H µ,p,β,D,δ ′ (x) := M H p,β,D,δ ′ (ρ(x, y)) dµ.
If t → γ(t) is a unit speed geodesic then for all t

(5.3) (E µ,p • γ) ′′ (t) ≤ H µ (γ(t)).
Proof. For x, y ∈ M , r = ρ(x, y), s → γ(s) = c(0, s) a unit speed geodesic started at x = c(0, 0), t → c(t, s) the geodesic satisfying c(1, s) = y, we have The main result is the following (see [START_REF] Le | Estimation of Riemannian barycentres[END_REF] for a similar result in a Riemannian manifold).

D ′′ p (0) ≤ pr p-2 D 2 max ((p -
Proposition 5.3. Assume -β 2 ≤ K ≤ k, -δ ≤ T ≤ δ ′ , C ≤ C and D ≤ D for some β, k, δ, δ ′ > 0 and C, D ≥ 1. Let p > 1.
Assume the support of µ is contained in B(x 0 , R) and E µ,p is strictly convex on B(x 0 , C(C + 1)R). Assume furthermore that the function H µ = H µ,p,β,D,δ ′ is bounded on B(x 0 , C(C + 1)R) by a constant C H > 0, and that the injectivity radius at any point of B(x 0 , C(C + 1)R) is larger than C 2 + C + 1. Define the gradient algorithm as follows:

Step 1 Start from a point x 1 ∈ B(x 0 , C(C + 1)R) such that E µ,p (x 1 ) ≤ R p (take for instance x 1 = x 0 ) and let k = 1.

Step 2 Let

(5.4) v k = grad(-E µ,p (x k ))) F (grad(-E µ,p (x k )))

, t k = F (grad(-E µ,p (x k ))) C H .

and let γ k be the geodesic satisfying γ k (0) = x k , γk (0) = v k . Define

(5.5)

x k+1 = γ k (t k )
then do again step 2 with k = k + 1.

Then the sequence (x k ) k≥1 converges to e p .

Proof. We first prove that the sequence (E µ,p (x k )) k∈N is nonincreasing. For this we write (5.6)

This proves that the sequence is nonincreasing. As a consequence, for all k ≥ 1, x k ∈ B(x 0 , C(C + 1)R), since E µ,p (x k ) ≤ R p and for all x ∈ B(x 0 , C(C + 1)R), E µ,p (x) > R p . Next we prove that E µ,p (x k ) converges to E µ,p (e p ). We know that E µ,p (x k ) converges to a ≥ E µ,p (e p ). Extracting a subsequence x k ℓ converging to some x ∞ ∈ B(x 0 , C(C + 1)R), this implies that t k ℓ converges to 0. But this is possible only if x ∞ = e p , which implies that a = E µ,p (e p ). As a consequence, any converging subsequence has e p a limit, and this implies that x k converges to e p . Remark 5.4. For this result we need the Hessian of E µ,p to be bounded, and the subgradient algorithm in Riemannian manifolds as developed in [START_REF] Yang | Riemannian median and its estimation[END_REF] does not work. The reason is that for this algorithm, we would need to take v k = grad---→ x k ep (-E µ,p (x k ))) F grad---→

x k ep (-E µ,p (x k ))

where grad---→ x k ep denotes the gradient with respect to the metric <•, •>---→ x k ep . So we would need to know e p !

(4. 6 )

 6 F µ (y) = F µx (y) + µ({x})ρ(y, x)

Remark 5 . 2 .

 52 [START_REF] Afsari | Riemannian L p center of mass : existence, uniqueness, and convexity[END_REF], rβ coth(rβ)) + δ ′ r . Integrating with respect to y this equation gives the result. If p ≥ 2 or µ has a smooth density then the function H µ is bounded on all compact sets.

E 2 k 2 ≤H 2 F 2 = 2 F 2 .

 222222 µ,p (γ k (t k )) ≤ E µ,p (γ k (0)) + dE µ,p , v k t k + C H t E µ,p (γ k (0)) -F (grad(-E µ,p (x k ))) 1 C H F (grad(-E µ,p (x k ))) + C (grad(-E µ,p (x k ))) C H E µ,p (γ k (0)) -C H (grad(-E µ,p (x k ))) C H

Corollary 5 . 5 . 1 C 2

 5512 Let p = 2. If R ≤ or M has nonpositive flag curvature, then the algorithm of Proposition 5.3 can be applied with the appropriate constants Proof. With this assumption, by Proposition 2.2 the function E µ,2 is strictly convex on B(x 0 , C(C + 1)R).

  we get Proposition 4.5. A point x in M is a local minimum of F µ if and only if (dF µx ) .

	(4.10) µ({x}) ≥ F Note that for the Riemannian case this result is due to Le Yang [34].
	Define the vector
	(4.11)	H(x) = grad y (F µx (y)) | y=x .

*


Christophe Krzeminski
email: christophe.krzeminski@isen.fr

Emmanuel Dubois

Xiaohui Tang

Nicolas Reckinger

André Crahay

Vincent Bayot

Optimisation and Simulation of an Alternative nano-flash Memory: the SASEM device

de niveau recherche, publiés ou non, émanant des établissements d'enseignement et de recherche français ou étrangers, des laboratoires publics ou privés.

Since the conventional floating gate device is believed to be hardly scalable below the 65 nm technology node, alternative storage structures for nonvolatile memories are strongly needed. The feasibility of a silicon-on-insulator nano flash memory device based on the differential oxidation rate resulting from gradients in the arsenic doping concentration has previously been reported [1,[START_REF] Tang | Fabrication, characterization and simulation of SOI single-electron devices[END_REF]. The key processes involved in the fabrication are arsenic implantation, lithography and wet oxidation. The resulting device is a triangular MOSFET with a nanocrystal floating gate embedded in the gate oxide. Our objective is to improve the reliability of the process and to ensure the presence of the memory dot for various conditions.

We investigate the wet oxidation step in details. A clear understanding of the mechanisms of dot formation and of the influence of self-limited oxidation effects on the final device should be undertaken. The various geometrical parameters (e.g dot surface, dot-channel distance ...) of the nano-device have been simulated as a function of the oxidation temperature and the duration of this oxidation step.

FABRICATION OF THE SASEM DEVICE.

We first brefly summarize the fabrication and the resulting S ¯elf A ¯ligned S ¯ingle E ¯lectron M ¯emory (SASEM). A SOI wafer with a 200 nm thick silicon overlay is used as a starting substrate. A high dose of Arsenic is implanted (1 ×10 15 cm -2 ) to create a very localized and a very doped area. Since the oxidation rate is doping dependent [START_REF] Ho | [END_REF]4], it will be possible to realize an anisotropic oxidation. Next a thin oxide film and a nitride layer is deposited on the substrate. E-beam lithography and reactive ion etching is used to pattern the device. Figure 1 presents the layout of the device. The central part of the device consists of a 150 nm by 150 nm square. Two 100 nm wide constrictions are defined in order to realise the connection between the source and drain regions and the central part. The various dimensions are defined such as a dot can be created in the central part of the device and is separated from the source/drain regions. Finally, a wet oxidation step is performed in order to create the dot. Figure 2 shows a SEM photography of the central part where a nano-floating gate is observed on top of the channel.

PROCESS SIMULATION ENVIRONMENT.

The main steps of the SASEM process have been described and the most critical step is the wet oxidation since it governs the creation of the silicon dot on top of the channel.

Oxidation is simulated by a standard Deal and Grove generalised viscoelastic approach [5,6] using a commercial simulator [START_REF]ISE DIOS, Process Simulator V6.0[END_REF]. We use the Ho and Plummer model [START_REF] Ho | [END_REF]4] in order to describe the enhancement of the oxide growth rate through the arsenic dopants (anisotropic oxidation). We used a standard set of parameters for the oxidation model which is known to give reasonable results on various configurations [START_REF] Senez | [END_REF]. Parameters are not so far from the standard parameters of the simulation tool [START_REF] Senez | [END_REF]. The objective is not a get a perfect agreement by adjusting the parameters on a specific experimental configuration of the SASEM device but to give raisonnable and interesting trends. Moreover, we must keep in mind that this process is very complex to simulate since the refinement of the grid must be preserved during the oxidation and the diffusion and segregation of arsenic have to be correctly described. We start with an anisotropic grid with a very high level of refinement in the region of high doping where the dot is expected to be created. This condition is very important. If the mesh describing the dot geometry is very relaxed, the study of the evolution of the dot properties could be very problematic. To keep the quality of the initial grid during oxidation, very small timesteps are performed ( inf than 1ms) in order to generate reduced and smooth adaptation of the grid. The main problem of the process simulation in general is the difficulty to access to specific geometrical information during the oxidation process. Furthermore, it is not fully adapted to the specificity of the SASEM memory cell.

To overcome this problem, a monitoring tool has been developed. The flowchart of the different steps is presented in figure 3. The monitoring tool generates a complex script with a large amount of commands. The basic idea is to perform a spatial sampling of the shape of the MESA-structure during the oxidation step. After the process simulation, the monitoring tool is able to treat all the flow of data generated to access to the evolution of the various parameters (dot length, dot width, dot-channel distance ...). We will see that the coupling between the process simulator and the monitoring tool is very efficient and allows a broad range of simulations.

SIMULATIONS.

Example of shape monitoring. The profile for zero minute of wet oxidation is not strictly straight since a short dry oxidation step (2 minutes) is performed before the wet oxidation. As shown by the two figures, the oxidation is very fast in the region where the arsenic concentration is high (the maximum for the arsenic profile is located at 240 nm). The dot is created after 52 minutes in the configuration. As shown by the figure 4, the dot is consumed mostly at the bottom where the arsenic concentration is large. The top of the dot is clearly well protected by the nitride mask. In comparison, the bottom of the channel is clearly much oxidized.

Evolution of dimensional parameters.

Various parameters have been studied. We report here only the evolution of the dot properties (figure 5). The thickness at the dot creation is nearly the same ∼ 50nm. This result indicates that the position of the maximum of the arsenic concentration governs the initial size of the silicon dot. There is a sharp decrease at the time of dot creation and next the dot size decreases almost linearly with a slope of about 0.5 nm/min. The evolution of the maximum dot width is also reported on the right figure.

Influence of the oxidation temperature on the final shape.

The evolution of the dot shape as a function of the oxidation temperature has been simulated. Figure 6 reports the evolution of both the channel and the dot shape. The time of oxidation for each temperature has been set just at dot creation. We clearly observe that the size of the dot is strongly reduced with an increasing temperature. In figure 6, the maximal half dot width is about 30 nm at 750 • C and reduce to less than 10 nm for 950 • C. For higher temperature, no dot is created. Furthermore, the oxidation temperature has clearly an impact on the shape that tends to be flat for a high temperature of oxidation. This result indicates that the oxidation temperature has an important influence on the final structure and that a wet low oxidation temperature is recommended.

CONCLUSIONS

Only a short part of the simulations performed has been described here. Thanks to an extensive set of simulations, we have now access to the evolution of many parameters and evolutions which otherwise would have required many experiments. More understanding of the dot formation is the key to improve the process. To conclude, we stress that the optimisation of the process is of crucial interest in order to allow the downscaling of this alternative nano flash memory. The simulation and the definition of a strategy for scaling the device is the next step to adress.

Analysis of results

Process simulator Script Specific file of results

Monitoring of SASEM creation during wet oxidation

Figure 4

 4 Figure 4 presents an example of the monitoring of the silicon shape during a wet oxidation at 800 • C. The silicon shape is represented for [0, 20, 40, 60, 80] minutes of wet oxidation.

FIG. 1 :

 1 Figures

FIG. 3 :

 3 FIG.3: Flowchart of the interaction between the monitoring tool and the process simulator.

FIG. 4 :

 4 FIG. 4: (Left) Initial configuration of the nanowire. (Right) Zoom on the evolution of the silicon shape of the mesa structure during a wet oxidation step. The shape variations can be easily monitored thank to the monitoring tool.

FIG. 6 :

 6 FIG.6: Evolution of the MESA structure (only the silicon part) and for various oxidation temperatures. The time of oxidation is set at the dot creation time.

ACKNOWLEDGMENTS.

This research is supported by the European Union through the IST-2001-32674 SASEM project.


Yaseen Thaer
email: y.thaer@iamb.it

Fatima Ayoub

Massimo Reverberi

Alessandra Ricelli

Anna Maria D'onghia

Thaer Yaseen

Early detection of Aspergillus carbonarius and A. niger on table grapes: a tool for quality improvement

Keywords: PCR, HPLC Additives/Contaminants: Mycotoxins -fungi, Ochratoxin Aspergillus carbonarius, Aspergillus niger, Ochratoxin A, Polyketide synthase, Real time PCR, Table grapes

Aspergillus carbonarius and A. niger aggregate are the main fungal contaminants of table grapes. Besides their ability to cause the Black rot, they can produce Ochratoxin A (OTA), a mycotoxin that has attracted increasing attention worldwide. The objective of this work was to set up a simple and rapid molecular method for the early detection of both fungi in table grapes before fungal development becomes evident. PCR-based assays were developed by designing species-specific primers based on the polyketide synthases (PKS S ) sequences of A. carbonarius and A. niger that have recently been demonstrated to be involved in OTA biosynthesis. Three table grape varieties (Red globe, Crimson seedless and Italia) were inoculated with A. carbonarius and A. niger aggregate strains producing OTA. The extracted DNA from control (non-inoculated) and inoculated grapes was amplified by PCR using ACPKS2F-ACPKS2R for A. carbonarius and ANPKS5-ANPKS6 for A. niger aggregate.

Both primers allowed a clear detection even in symptomless samples. Since the high specificity and sensitivity of PCR-based methods, they are considered to be a good alternative to traditional diagnostic means for the early detection of fungi in complex matrix. The results obtained could be useful for the definition of a "quality label" for tested grapes to improve the safety measures taken to guarantee the production of fresh table grapes.

Introduction

Many fungi may occur on table grapes in the vineyard, but the main concern from the point of view of the mycotoxin contamination is represented by the group of black Aspergillus species, in particular Aspergillus carbonarius and Aspergillus niger aggregate [START_REF] Ailsa | Fungi and mycotoxins in vineyards and grape products[END_REF].

They are Deuteromycotina fungi (Imperfect fungi, asexual reproduction) that belong to the order of Hyphomycetes, family of Moniliaceae, section Nigri [START_REF] Atoui | Approche de la mycotoxinogenese chez Aspergillus ochraceus et Aspergillus carbonarius: Etudes moleculaire et physiologique[END_REF]. In addition to their ability to cause Black rot diseases in the field, they can produce ochratoxin A (OTA), a mycotoxin that is attracting increasing attention worldwide (Varga et al, 2006).

Ochratoxin A, the most toxic of the ochratoxins, is produced by several fungal species belonging to Aspergillus and Penicillium genera (Varga et al. 2006). This mycotoxin has a significant economic impact and poses a serious problem for human health. It is named after Aspergillus ochraceus, the first rot agent from which it was isolated (Van der Merwe et al, 1965), although later on other genera were reported to be capable of producing this toxin. The compound 7-carboxyl-5-chloro-3,4-dihydro-8-hydroxyl-3R-methylisocoumarin-7-Lphenylalanine (C 20 H 18 ClNO 6 ), usually known as OTA, can be found in several food and drink commodities, including cereals, coffee and wine, at different levels [START_REF] Battaglia | Occurrence and significance of ochratoxin A in food[END_REF].

Several studies have been performed to investigate the contamination with black aspergilli on wine grapes [START_REF] Battilani | Ochratoxin A production by Aspergillus carbonarius on some grape varieties grown in Italy[END_REF][START_REF] Leong | Survival and growth of Aspergillus carbonarius on wine grapes before harvest[END_REF]), but only few works were carried out on table grapes [START_REF] Guzev | Occurrence of ochratoxin A producing fungi in wine and table grapes in Israel[END_REF]. In this work A. carbonarius and A. niger aggregate strains able to produce OTA were frequently isolated from table grapes.

Most research efforts were concentrated on the means for prevention of OTA biosynthesis as the best defence for protecting the consumer. However, prevention is not always possible, especially from OTA biosynthesis occurring under field conditions. In fact, OTA remains in processed food and the recovery of infected products is only possible through their decontamination (Bellì, 2006). Conversely, in the case of table grapes, the prevention remains the only possible way to control OTA formation; this point can be achieved only thanks to the early detection of the OTA-producer fungus. Traditional diagnostic methods used in food mycology are based on macroscopic and microscopic features and culture in appropriate media. On the one hand, the development of fruiting structures requires 2 to 10 days of culture on different media [START_REF] Raper | The genus Aspergillus[END_REF], increasing considerably the time of analysis. On the other hand, these methods have shown a low degree of sensitivity and are difficult to standardize (Zhao et al, 2001).

Misidentification can often occur because some fungi may be poorly characterized or because a considerable expertise is required [START_REF] Patiño | PCR detection assays for the ochratoxin-producing Aspergillus carbonarius and Aspergillus ochraceus species[END_REF].

The increasing use of molecular methods in fungal diagnosis has provided tools for answering taxonomic questions that morphological procedure have left unsolved [START_REF] Taylor | The evolutionary biology and population genetics underlying fungal strain typing[END_REF][START_REF] Parenicova | Integration of modern taxonomic methods of Penicillium and Aspergillus classification[END_REF]. PCR-based methods that target on DNA are considered a good alternative for rapid and early diagnosis because of their high specificity and sensitivity [START_REF] Bluhm | Multiplex polymerase chain reaction assay for the differential detection of trichothecene-and fumonisin-producing speciesof Fusarium in cornmeal[END_REF].

Reduction of OTA contamination may be accomplished by the early detection of OTAproducing fungal species in order to intervene at the right time to reduce the source of OTA.

Several approaches have been followed for the molecular detection of the OTA-producing fungi. [START_REF] Atoui | Targeting a polyketide synthase gene for Aspergillus carbonarius quantification and ochratoxin A assessment in grapes using real-time PCR[END_REF]

Materials and methods

Fungal strains and culture conditions

Each grape sample (~500 g) was washed with 500 ml of sterile distilled water on a shaker for 30 min. Aliquots (100µl) of washing water and it's serial dilution were spread on Czapek-Dox-Agar medium (CDA) with 0.02% of Streptomycin sulphate.. The plates were incubated at 25°C

and daily checked until fungal colonies formation. Each colony was transferred on a Petri dish with the same medium and incubated for 7 days at 25°C. Aspergillus strains were identified according to their morphological characteristics (macroscopically and microscopically), following the rules and descriptions provided by [START_REF] Raper | The genus Aspergillus[END_REF] and Pitt and Hocking 1997 and in the other hand following molecular tools. The strains were also tested to determine their ability to produce OTA. The black aspergilli strains have been preserved in slants with CDA medium at 4°C.

All strains used in this study were OTA producers. The different strains were isolated from grapes in Apulia, Italy. Fungal cultures were maintained on CDA medium and stored on the same medium at 4°C.

Grapes infection

Berries from three glass beads (425-600 µm diameter) and two sterile 5 mm stainless steel balls. This mixture was vortexed at 3,000 rpm for 10 min, and centrifuged for 15 min at 13,000 rpm. The supernatant was mixed with an equal volume of Chloroform/Isoamyl alcohol (24:1 v/v), vortexed at 3,000 rpm for 2 min, centrifuged for 5 min at 13,000 rpm and precipitated with two volumes of 100% cold isopropanol (-20°C). The precipitated DNA was washed with 70% (-20°C) ethanol, dissolved in 50 µl of TE buffer (10 mM Tris-HCL, 0.1 Mm EDTA, pH 7.6) and quantified using a spectrophotometer (Yaseen et al. 2004).

DNA extraction from grape samples

An aliquot of 250 mg of lyophilized homogenized grapes were suspended in 1 ml of lysis buffer (200 mM of Tris HCl pH 8, 25 mM of NaCl, 25 mM of EDTA and 0.5% of SDS) and vortexed at 3,000 rpm for 10 min. After 1 hour at 65°C, and 10 minutes in ice, the sample was centrifuged at 13,000 rpm at 4°C for 15 minutes. The aqueous phase was mixed with 200 µl of sodium acetate 4M, put in ice for 30 minutes and centrifuged at 13,000 rpm at 4°C for 15 minutes. The recovered supernatant was homogenized with an equal volume of Phenol/Chloroform/Isoamyl alcohol (25:24:1 v/v/v), vortexed for 5 minutes and centrifuged at 13,000 rpm at 4°C for 10 min. An equal volume of chloroform 100% was added to the aqueous phase, then gently vortexed for 5 min and centrifuged at 12,000 rpm at 4°C for 15 min.

The DNA contained in the aqueous phase was precipitated by adding half a volume of Isopropanol and stored for 20 min at -20°C. The pellet obtained after 15 minutes of centrifuge at 13,000 rpm at 4 °C was washed with ethanol 70% and centrifuged at 13,000 rpm for 15 min at 4°C. The pellet obtained was dried, dissolved in 28 µl of sterile distilled H 2 O and 2 µl of RNAse (20 mg/ml) and stored at 37°C for 30 min.

Primer design and PCR reaction

Aspergillus carbonarius-specific primer PKS genes encode the polyketide synthase (PKS) enzyme, which is involved in the biosynthesis of secondary metabolites. In the gene bank, there are only five partial sequences of A. carbonarius PKS gene [START_REF] Atoui | Targeting a polyketide synthase gene for Aspergillus carbonarius quantification and ochratoxin A assessment in grapes using real-time PCR[END_REF]). On the basis of these sequences, a degenerated primer (PKS1) for Aspergillus spp. and a specific primer for Aspergillus carbonarius (ACPKS2) were designed (Table 1).

Aspergillus niger-specific primer

The degenerated primer for Aspergillus spp. PKS1 was used to amplify a partial sequence of the PKS gene of A. niger (Yaseen et al, 2009). On the basis of this sequence, a couple of primers have been designed to amplify a DNA fragment of 157 bp: ANPKS 5 (5' ACGGTAAACGTCCTGGATGA 3') and ANPKS 6 (3' CGTGCTGTTGAAGCCACTT 5').

PCR amplification

Specific PCR assays have been developed in this study for the detection of both A. [START_REF] Atoui | Targeting a polyketide synthase gene for Aspergillus carbonarius quantification and ochratoxin A assessment in grapes using real-time PCR[END_REF]. At the end of each run, and to ensure that only the specific PCR product was present, a Melting curve was performed which was programmed as 81 cycles of 30 seconds at 55 °C.

Results and discussion

Specificity of the PCR reaction

A single fragment of 330 bp of the target DNA was obtained only when genomic DNA from A. carbonarius was used with ACPKS2 (F/R) and of 157 bp from A. niger with ANPKS (5/6) primers (Figure 1 A/B). The absence of any aspecific band using DNA extracted from other Aspergillus spp. and from filaments fungi demonstrate the specificity of primers, and the possibility of application of quantities real time PCR.

Application of the qPCR method for A. carbonarius and A. niger in artificially

contaminated samples by interpolating its Ct value against the standard curve. The results were expressed in ng DNA/µl. The amount of DNA of both fungi did not statistically change until 26 days of post-harvest storage at 4°C, then the quantity of DNA of both fungi significantly increased after 40 days of storage in case of Italia and Red Globe but not in Crimson (Figure 3). These results could be partly explained considering the post-harvest aging processes. In fact, during storage time, the grape berries become more susceptible to the attack by these fungi. In addition, these results are correlated with those of visual infection,

showing a high percentage of infection: up to 70% in berries stored for 40 days then incubated at 30°C for 12 days. Hence, the development of this fungus was much rapid reaching high percentage of infection after 26 days of post-harvest cold storage. A gradual increase in the DNA content of A. carbonarius and A. niger in each variety during the incubation periods at 30°C (0, 3, 6, 9 and 12 days) was monitored after 40 days of post-harvest cold storage (Figure 4 A/B).The target fungal DNA, of both A. carbonarius and A. niger was detected even at time 0 (without incubation), the lowest detection limit achieved was 2.37 10 -2 ng/µl showing the high sensitivity of the detection method (Figure 4 A/B).

Optimization of a Real-time PCR assay is very important, mainly when SYBR Green dye is used. The analysis of the melting curve and the parameters calculated from the standard curve showed that the designed method was highly optimized. The melting curve for A. carbonarius showed that only one fragment was amplified and primers did not form other aspecific bands (data not shown). It was evident that both A. carbonarius and A. niger were detected in all samples, nevertheless a superficial sterilization was performed. These results correlate with the visual observations (Figure 5); fungal growth was observed not only nearby the zone of inoculation, but also around the stalk of each berry, where, usually, the fungus starts to develop in non inoculated samples. Moreover, according to [START_REF] Battilani | Occurrence of ochratoxin A-producing fungi in grapes grown in Italy[END_REF], A. carbonarius, as well as A. niger can enter the berries even without skin damage.

Several amplification methods have been developed with conventional PCR for the analysis of phytopathogenic fungi (competitive PCR) [START_REF] Hadidi | Molecular Methods in Plant Pathology[END_REF][START_REF] Mahuku | Quantifying Verticillium dahliae in soils collected from potato fields using a competitive PCR assay[END_REF].

However, these methods are laborious and insufficiently accurate because the amplification efficiency decreases in later PCR cycles [START_REF] Ginzinger | Gene quantification using real-time quantitative PCR: an emerging technology hits the main stream[END_REF]. An accurate, reliable, and highthroughput quantification of target DNA requires a real-time PCR approach [START_REF] Lie | Advances in quantitative PCR technology: 5' nuclease assays[END_REF][START_REF] Schmittgen | Real-time quantitative PCR[END_REF]. Real-time PCR quantifies initial concentrations of DNA by the threshold cycle (Ct), the amplification cycle at which the concentration of DNA produced in the PCR reaction gives a visualized fluorescence climbs above a baseline "threshold". Therefore, a lower Ct corresponds to higher amount of DNA in the beginning of the reaction.

Concentration values obtained by Real-time PCR in samples with a starting quantity of 5 ng of total DNA extracted from samples were related with the initial amounts of A. carbonarius and A. niger DNA. Moreover, the application of the method to quantify these fungi in artificially contaminated matrices was confirmed and detection was executed at longer incubation times at 30 °C. Primer Forward 5'-3' Reverse 3'-5' Fragment

  Aspergillus carbonarius and A. niger aggregate strains were cultured in 100 ml Erlenmeyer flasks containing 20 ml Potato Dextrose-Broth with 4-5 mycelial disks cut from the border of 4 to 5-days old colonies on PDA plates, and incubated at 26°C on shaker at 120 rpm. Mycelia from 4 days-old cultures were harvested by filtration through Whatman paper n° 1 and kept at -20°C for DNA extraction. Five to ten mg of the conserved mycelia were collected with a spatula, suspended in 400 µl of breaking buffer and extracted with 400 µl of Phenol/Chlorophorm/Isoamyl alcohol (25:24:1 v/v/v) in the presence of 25 mg of acid washed

  was obtained by plotting the Ct value versus the logarithm of the concentration of each DNA dilution. The linear correlation coefficient of the standard curve was R 2 =0.963 and R 2 =0.993 for A. carbonarius (Figure 2/A) and A. niger (Figure 2/B) respectively, demonstrating the accuracy of PCR-based quantification. Consequently, the use of the generated standard curves allows quantifying DNA of A. carbonarius and A. niger from contaminated table grape vine.

  Quantitative detection by Real-time PCRA Real-time PCR method was developed using SYBR Green dye to detect A. carbonarius and A. niger OTA-producing strains on table grapes. DNA extractions were performed from 191 samples (from the three tested varieties), and subject to real-time PCR amplification.

	Quantitative real-time PCR reaction was incubated in a programmable thermal cycler (iQ5-
	Bio-rad) starting from 5 min denaturation at 95°C followed by 35 cycles at 95°C for 10 s,
	annealing at 62°C for 10 s for A. carbonarius and 52°C for A. niger, extension at 72°C for 15 F s. Reactions were performed in 25 µl mixture containing: 12.5 µl of 1X iQSYBR Green o r Supermix (Bio-Rad), 10µM forward and reverse primers, 5 ng of DNA template and
	P
	e
	e r
	R
	e
	v i e
	w
	the following PCR conditions: 95°C for 300 s; 35 cycles at 95°C for 10 s, 62°C for 10 s for A. carbonarius and 52°C for A. niger, 72°C for 15 s; plus a final extension cycle at 72°C for 300 O n l
	s. To visualize the amplification of the expected size fragment, 8 µl of each PCR sample was y
	analyzed by agarose gel electrophoresis. Amplicons were separated by electrophoresis in
	1.5% agarose gels in TAE buffer (1x) at 80V for 10 600 s followed by 100 V for 45 min, and
	stained for 10 min with ethidium bromide (2 µg ml -1 ).
	7

carbonarius and A. niger, the main source of OTA contamination on table grapes. To demonstrate the specificity of these primers, PCR amplification of genomic DNA from various fungal strains has been performed; PCR was executed in 25 µl mixture containing 100ng of genomic DNA, 5 µl of 5x Colorless Go Taq Flexi buffer (Promega, Madison WI-U.S.A), 100 µM of dNTPs buffer (Promega, Madison WI-U.S.A), 2.5 mM MgCl 2 (Promega, Madison WI-U.S.A), 10 µM each of forward and reverse primers and 2.5 unit of GoTaq® DNA polymerase (5 U/µl Promega, Madison WI-U.S.A). The reactions were performed under deionised sterilized water up to 25 µl. Standard curve was generated from the amplification of ten-fold dilutions of A. carbonarius and A. niger DNA ranging from 10 ng to 0.01ng. These DNA amounts were subject to real-time PCR conjugated with SYBR Green dye using the specific primers ACPKS2 F/R for A. carbonarius and ANPKS 5/6 for A. niger. The parameter used for data interpretation was the threshold cycle (Ct value), which represents the cycle number at which the fluorescence emission exceeds the fixed threshold. Quantification of the unknown DNA samples was carried out by interpolating the samples Ct values with the standard curves

Table 1 .

 1 The specificity and high degree of sensitivity of real time PCR detection method developed for A. carbonarius and A. niger provide a good tool for the early detection of this OTA producing fungi in table grapes and to prevent OTA entering the food chain. Detection of these fungi, in the case of table grapes, becomes particularly critical around harvest time when the OTA production is considered high.In this work we described a real-time PCR method that lead to the identification and quantification of Aspergillus carbonarius and Aspergillus niger aggregate in three varieties of table grapes without visible symptoms of rotting. This method was rapid, specific with high sensitivity, providing a useful tool for early detection of both fungus and prediction of OTA contamination in table grapes susceptible to be colonized by these species. Moreover, the results obtained could be useful to define a "quality label" for the assayed grapes in order to improve the safety measures taken to guarantee the production of fresh table grapes. Agarose gel electrophoresis of amplification products from genomic DNA of various fungal isolates obtained with primer pairs ACPKS2 (F/R) specific for A. carbonarius, amplifying a fragment of 330 bp (A); primer pairs ANPKS 5/6 specific for A. niger, amplifying a fragment of 157bp. Lanes from 1 to 17 are Alternaria sp., Botrytis cinerea, Cladosporium sp., Cryphonectria parasitica, Cylindrocarpon sp., Gliocladium sp., Cladosporium sp., Colletotrichum sp., Penicillium sp., Phomopsis viticola, Sclerotinia sclerotiorum, Dothiorella ribis, Septoria tritici, Trichoderma sp., Phoma sp., Macrophomina sp. and A. niger respectively. Lane 18/A: A. carbonarius and lane 19/A: Non-template control, lane 18/B: Non-Template control and lane 19/B: A. niger. Correlation between Log starting quantity of DNA ng/µl and cycle threshold (Ct) value as assessed by SYBR green Real-time PCR of A. carbonarius (A) and of A. niger (C). Standard curve generated from the amplification of ten-fold dilutions (from 10 ng/µl to 0.01 ng/µl) from target genomic DNA of A. carbonarius (B) and of A. niger (D), the coefficient of linear regression equation of Ct with Log DNA concentrations was R 2 = 0.963 and R 2 = 0.993 respectively. Detection of DNA quantity of Aspergillus carbonarius and A. niger by SYBR green Real-time PCR extracted from inoculated and non-inoculated table grapes of Red globe, Italia and Crimson seedless, starting from the day of harvest and during 40 days of post-harvest cold storage. Columns report average of DNA content from four different DNA extractions and amplifications. Bars represent standard deviation. and 40 (D) days of postharvest cold storage during different incubation periods (0, 3 , 6, 9 and 12 days) at 30 °C. Un-inoculated grapes from each cultivar were used as controls (t). Sequence of degenerated primers for Aspergillus spp. (PKS1) and specific primers for A. carbonarius (ACPKS2).

	F o r
	P
	e
	e r
	R
	e
	v i e
	w
	O n l
	y
	10

Figure 5. Percentage of infected grapes Red globe (R), Italia (I) and Crimson seedless (C) cultivars with A. carbonarius (c) and A. niger (n), after 0 (A), 13 (B), 26 (C)

http://mc.manuscriptcentral.com/tfac Email: fac@tandf.co.ukFood Additives and Contaminants