Vincent Briane

Charles Kervrann

Myriam Vimond

A STATISTICAL ANALYSIS OF PARTICLE TRAJECTORIES IN LIVING CELLS

Keywords: MSC 2010 subject classifications: Primary 62G10, 60J60; secondary 92B05 Three-Decision Test, Multiple Hypothesis Testing, Diffusion Processes, Intracellular Particle Trajectories 1

 to fit with the three-decision test setting, in order to apply the test procedure to a collection of independent trajectories. The performance of our procedure is much better than the MSD method as confirmed by Monte Carlo experiments. The method is demonstrated on real data sets corresponding to protein dynamics observed in fluorescence microscopy.

Introduction.

A cell is composed of lots of structures in interaction with each other. They continuously exchange biological material, such as proteins, directly via the cytosol or via networks of polymerised filaments namely the microtubules, actin filaments and intermediate filaments. The dynamics of these proteins determine the organization and function of the cell (Bressloff, 2014, chapter 9). The underlying traffic is known to be oriented and it is established that local dynamics of proteins obey to biophysical laws, including subdiffusion (diffusion in a closed domain or in an open but crowded area), free diffusion (or Brownian motion) and superdiffusion (active transport along the microtubules). Then, inference on the modes of mobility of molecules is central in cell biology since it reflects the interaction of the structures of the cell. For instance the postsynaptic AMPA-type glutamate receptor (AMPAR) is a protein involved in the fast excitatory synaptic transmission : it plays a crucial part in many aspects of brain functions including learning, memory and cognition. The dynamics of AMPARs determine synaptic transmission : aberrant AMPAR trafficking is implicated in neurodegenerative process, see [START_REF] Henley | Routes, destinations and delays: recent advances in AMPA receptor trafficking[END_REF]. [START_REF] Hoze | Heterogeneity of AMPA receptor trafficking and molecular interactions revealed by superresolution analysis of live cell imaging[END_REF] model their motion with diffusions confined in a potential well. As an other example, [START_REF] Lagache | Quantitative analysis of virus and plasmid trafficking in cells[END_REF] modelled the dynamics of a virus invading a cell to infer its mean arrival time to the cell nucleus where it replicates. In their model, the dynamic of the virus alternates between superdiffusion and Brownian motion.

The biophysical literature uses the word diffusion to describe the dynamics of intracellular particles. This word encompasses a great variety of models. First, we have to distinguish diffusion at the macroscopic and microscopic level. The two are closely related in the sense that diffusion at the macroscopic level gives the dynamic of population (or concentration) of particles undergoing microscopic diffusion. In the sequel, we are only interested in microscopic diffusion i-e diffusion of individual trajectories. Secondly, among the variety of mathematical models associated to microscopic diffusion alone, we focus on models that reflects one of the three dynamics : free diffusion, superdiffusion and subdiffusion. Usually, in the biophysical literature, the definition of these dynamics is related to the criterion of the mean square displacement (MSD) [START_REF] Qian | Single particle tracking. Analysis of diffusion and flow in two-dimensional systems[END_REF]. According to the shape of the MSD curve, the literature distinguishes four main types of diffusion [START_REF] Saxton | Single-particle tracking: applications to membrane dynamics[END_REF][START_REF] Monnier | Bayesian approach to MSD-based analysis of particle motion in live cells[END_REF]: superdiffusion, free diffusion, confined diffusion and anomalous diffusion. A particle evolving freely inside the cytosol or along the plasma membrane is modelled by free diffusion. Its motion is due to the constant collisions with smaller particles animated by thermal fluctuations. Then, the particle does not travel along any particular direction and can take a very long time to go to a precise area in the cell. Active intracellular transport can overcome this difficulty so that motion is faster and direct specific. The particles (called in this context cargo) are carried by molecular motors along microtubular filament tracks. Superdiffusions model the motion of molecular motors and their cargo. Subdiffusion, which includes confined diffusion and anomalous diffusion, are the translations of several biological scenarios. Confined or restricted diffusion [START_REF] Metzler | The random walk's guide to anomalous diffusion: a fractional dynamics approach[END_REF][START_REF] Hoze | Heterogeneity of AMPA receptor trafficking and molecular interactions revealed by superresolution analysis of live cell imaging[END_REF] is characteristic of trapped particles : the particle encounters a binding site, then it pauses for a while before dissociating and moving away. Anomalous diffusion includes parti- , the blue trajectory is Brownian, the red one is from a fractional Brownian motion (parameter h > 1/2), in cyan it is from a Ornstein-Uhlenbeck process and in green from a fractional Brownian motion (parameter h < 1/2). The parameters of the processes are given in Table (4).

. cles which encounters dynamic or fixed obstacles [START_REF] Saxton | Anomalous diffusion due to obstacles: a Monte Carlo study[END_REF][START_REF] Berry | Anomalous diffusion due to hindering by mobile obstacles undergoing Brownian motion or Ornstein-Ulhenbeck processes[END_REF], or particles slowed by the contrary current due to the viscoelastic properties of the cytoplasm. In this paper, we will not distinguish confined and anomalous diffusion and consider that both are subdiffusion. In order to classify the observed trajectories into the three aforementioned types of diffusion i-e subdiffusion, superdiffusion or Brownian motion (see Figure 1), a popular statistic is the pathwise Mean Square Displacement (MSD). The simplest rule to classify a trajectory with the MSD is based on the least-squares estimate of the slope β of the log-log plot of the MSD versus time [START_REF] Feder | Constrained diffusion or immobile fraction on cell surfaces: a new interpretation[END_REF]. [START_REF] Didier | The asymptotic distribution of the pathwise mean-square displacement in single-particle tracking experiments[END_REF] study the limiting distribution of the pathwise MSD according to the true value of β. Nevertheless, MSD has several limitations :

• The MSD statistic is a summary statistic, and does not suffice to characterize the dynamics of the trajectory. [START_REF] Gal | Particle tracking in living cells: a review of the mean square displacement method and beyond[END_REF] present several other statistics which can be associated to MSD for trajectory analysis. [START_REF] Lund | SpatTrack: An imaging toolbox for analysis of vesicle motility and distribution in living cells[END_REF] propose a decision tree for selection motion model combining MSD, Bayesian information criterion and the radius of gyration. [START_REF] Lysy | Model comparison and assessment for single particle tracking in biological fluids[END_REF] present a likelihood-based inference as an alternative to MSD for the comparison between two models of subdiffusions : fractional Brownian motion and a generalized Langevin equation. They consider a Bayesian model to estimate the parameter of the diffusion and they use the Bayes factor to compare the models.

• The variance increases with the time lag (see Figure 2) : only the first few points of the MSD may be used to estimate the slope. Moreover the MSD variance is also severely affected at short time lags by dynamic localization error and motion blur. [START_REF] Michalet | Mean square displacement analysis of single-particle trajectories with localization error: Brownian motion in an isotropic medium[END_REF] details an iterative method, known as the Optimal Least Square Fit (OLSF) for determining the optimal number of points to obtain the best fit to MSD in the presence of localization uncertainty.

• In order to take account of the variance of the MSD estimate, several authors use a set of independent trajectories rather than single trajectories. These trajectories may have different lengths but are assumed to have the same kind of motion. For instance [START_REF] Pisarev | Numerical analysis of particle trajectories in living cells under uncertainty conditions[END_REF] consider weighted-least-square estimate for β by estimating the variance of pathwise MSD. Their motion model selection is then based on the modified Akaike's information criterion. [START_REF] Monnier | Bayesian approach to MSD-based analysis of particle motion in live cells[END_REF] propose a Bayesian approach to compute relative probabilities of an arbitrary set of motion models (free, confined, anomalous or directed diffusion).

In a general situation, this averaging process can lead to oversimplification and misleading conclusions about the biological process [START_REF] Gal | Particle tracking in living cells: a review of the mean square displacement method and beyond[END_REF].

In this paper, we propose a measure that circumvents some limitations of the MSD and which is efficient for classifying single trajectories. In our experiments, the raw data consist in a time series of images obtained by time-lapse fluorescent microscopy. They depict small spots or particles trafficking at the scale of a cell over a fixed time period. Single-and multiple-particle tracking procedures allow us to extract individual molecule and protein trajectories, see Figure 1 (b), which are the data on which we infer. These tracking procedures involve mathematical models and scientific computing to track the single-particle trajectories, see [START_REF] Chenouard | Objective comparison of particle tracking methods[END_REF]; [START_REF] Maroulas | Tracking rapid intracellular movements: A Bayesian random set approach[END_REF]; [START_REF] Brunel | A tracking approach to parameter estimation in linear ordinary differential equations[END_REF] and references therein. Our procedure is a threedecision test [START_REF] Shaffer | Control of directional errors with stagewise multiple test procedures[END_REF]. We previously proposed a test with promising results to address this classification problem [START_REF] Briane | An adaptive statistical test to detect non Brownian diffusion from particle trajectories[END_REF]. The null-hypothesis is that the observed trajectory is generated from a Brownian motion and the two distinct alternatives are the subdiffusion and superdiffusion. The test statistic T n is the standardized largest distance covered by the particle from its starting point. We interpret this measure as follows : i/ if the value T n is low, it means that the process stayed close to its initial position and the particle may be trapped in a small area or hindered by obstacles (subdiffusion); ii/ if the value T n is high, the particle went far to its initial position and the particle may be driven by a motor in certain direction (superdiffusion). We propose a theoretical study of the distribution of the test statistic under the null hypothesis, and we assess the power of the test under parametric models illustrating superdiffusion and subdiffusion. At the end, we derive a multiple test procedure in order to apply simultaneously the test on a collection of independent trajectories which are tracked inside the same living cell. This procedure is an adaptation of the procedure of [START_REF] Benjamini | On the adaptive control of the false discovery rate in multiple testing with independent statistics[END_REF].

The present paper is organized as follows. In Section 2, we describe the inference model and provide some examples of subdiffusion and superdiffusion. Our testing procedure is defined in Section 3. In Section 4, we derive a multiple testing procedure for a collection of trajectories. We carry out a simulation study and illustrate the method on real data in Section 5. We focus on the analysis of the Rab11a GTPase protein. This protein is involved in the trafficking of molecules from the endosomes located inside the cell to the cell plasma membrane. The data are TIRF microscopy images depicting the last steps of exocytosis events observed in the region very close the plasma membrane [START_REF] Schafer | Rab11-FIP2 interaction with MYO5B regulates movement of Rab11a-Containing recycling vesicles[END_REF]. The proofs are postponed to the appendix.

2. Diffusion models for particle trajectories. We observe the successive positions of a single particle in a two-dimensional space at times t 0 , t 1 , . . . , t n . We suppose that the lag time between two consecutive observations is a constant ∆. The observed trajectory of the particle is,

X n = (X t 0 , X t 1 , . . . , X tn ) , where X t i = (X 1 t i , X 2 t i ) ∈ R 2
is the position of the particle at time t i = t 0 + i∆, i = 0, . . . , n. This discrete trajectory is generated by a stochastic process (X t ) t 0 ≤t≤tn with continuous path, solution of the stochastic differential equation (SDE) :

(2.1)

dX i t = µ i (X i t )dt + σdB h,i t , i = 1, 2.
where B h,i t are unobserved independent 2D-fractional Brownian motions of unknown Hurst parameter h, σ > 0 is the unknown diffusion coefficient and (µ 1 (x 1 ), µ 2 (x 2 )) : R 2 → R 2 is the unknown drift term.

Assumption 1. We assume that µ i fulfils the linear growth hypothesis :

(2.2) ∃K > 0, ∀x ∈ R, |µ i (x)| ≤ K(1 + |x|),
and the Lipschitz condition :

(2.3) ∃M > 0, ∀(x, y) ∈ R 2 , |µ i (x) -µ i (y)| ≤ M |x -y|.
We denote by L the set of functions verifying Assumption 1. Assumption 1 is sufficient to ensure that SDE (2.1) admits a strong solution (see [START_REF] Nualart | Regularization of differential equations by fractional noise[END_REF] for the case 0 < h ≤ 1/2, and Mishura (2008, Chapter 3) for the case 1/2 < h < 1). For a given fractional Brownian motion, we say that (X t ) is a strong solution of the SDE (2.1) if (X t ) verifies (2.1), has continuous paths and that, at time t, X t depends only on X t 0 and on the trajectory of the fractional Brownian motion up to time t. In the following, P h,µ,σ denotes the measure induced by the stochastic process (X t ) solution of (2.1). This measure comprises all the finite dimensional distributions of the process that is the distribution of the vectors (X t 0 , . . . , X tn ), n ∈ N * and t 1 <, . . . , < t n . We also note P = {P h,µ,σ : 0 < h < 1, µ ∈ L, σ > 0} the set of solutions of the SDE (2.1).

Remark 2.1. In the following, we adopt the large-sample scheme to derive asymptotic properties of our procedure, that is the inter-observation time ∆ remains fixed and the number of observations n tends to infinity. In the experimental context of microscopic sequences, ∆ is the resolution of the camera while n is the number of frames during which we track the particle. Obviously the resolution of the camera does not change during the experiment. In an ideal situation, we track the particle during an infinite time of observation therefore the number of frames n tends to infinity. Then the large-sample scheme is the most realistic scheme in our context. Other schemes exist (see (Fuchs, 2013, Section 6.1.3)) as the high-frequency scheme for which ∆ tends to zero while the duration of observation is fixed. (Fuchs, 2013, Section 6.1.3) emphasizes that the large-sample scheme is the most realistic in real applications while the high-frequency scheme is convenient from a theoretical point of view.

Heuristically, a SDE models the motion of a particle in a fluid submitted to a deterministic force due to the fluid and a random force due to random collisions with other particles. That is why we model efficiently the motion of intra-cellular particles with these processes. In Equation (2.1), the velocity of the fluid is given by the drift µ while the term σdB h t expresses the random component of the motion due to random collisions.

2.1. Mean Square Displacement and Modes of Motion. The mean square displacement is a second order moment defined as :

(2.4) MSD(t) = E h,µ,σ X t -X t 0 2 2 ,
where E h,µ,σ is the expectation associated to the measure P h,µ,t . The empirical mean square displacement of the trajectory X n at lag j is given by :

(2.5)

MSD(j∆) = 1 n -j + 1 n-j k=0 X t k+j -X t k 2 .
According to the shape of the MSD curve, it is common to classify the motions into four classes (rather than into the three groups we presented until now) : Brownian motion (B), anomalous diffusion (A), confined diffusion (C) and directed diffusion (D) [START_REF] Monnier | Bayesian approach to MSD-based analysis of particle motion in live cells[END_REF] Saxton andJacobson (1997) Pisarev et al. (2015)),

MSD B (t) = 2σ 2 t, (2.6) MSD A (t) = 2σ 2 t β , (2.7) MSD C (t) = r 2 c a (1 -be -cσ 2 /(2r 2 c ) ), (2.8) MSD D (t) = v 2 t 2 + 2σ 2 t, (2.9)
where β, a, b, c, r c are positive constants, and v ∈ R 2 . In the following subsections, we present diffusion models defined by a SDE (2.1) which have a MSD of the type (A), (B), (C) or (D). We also associate each situation (A), (B), (C) and (D) to one of the three group of diffusions we defined in first place.

Figure 2 illustrates one of the limitation of the MSD : the variance of the estimator (2.5) increases with the lag j. From the MSD model (2.7), [START_REF] Feder | Constrained diffusion or immobile fraction on cell surfaces: a new interpretation[END_REF] states that the trajectory is subdiffusive if β < 0.9, superdiffusive if β > 1.1 and Brownian if 0.9 < β < 1.1. If β < 0.1 it states that the particle does not move. Figure 2 suggests that the classification of [START_REF] Feder | Constrained diffusion or immobile fraction on cell surfaces: a new interpretation[END_REF] overdetects subdiffusion and superdiffusion while it is Brownian motion. 2.2. Free diffusion. Free diffusion or Brownian motion is the most popular process for describing particle motion suspended in a liquid [START_REF] Einstein | Investigations on the Theory of the Brownian Movement[END_REF]. It suits particularly well for describing intracellular particle motion as the interior of the cell is mainly made of a fluid called the cytosol. Brownian motion allows dissolved macromolecules to be passively transported without [START_REF] Feder | Constrained diffusion or immobile fraction on cell surfaces: a new interpretation[END_REF], t → t β , β = 0.9 and 1.1. The dashed lines are the pointwise high probability interval of 95% associated to the empirical MSD curve for a standard Brownian motion trajectory of length n = 30. The bounds of the interval are the 2.5% and 97.5% empirical quantile of (2.5) and are computed by Monte Carlo simulation from 10 001 Brownian trajectories of size n = 30.

any input of energy. In the SDE (2.1), it matches with the situation where the drift µ i = 0 and h = 1/2. In biophysics, it is the process of reference. As a consequence they define the concepts of subdiffusion and superdiffusion from the MSD of Brownian motion (2.6).

2.3. Subddiffusion. Anomalous and confined diffusion respectively defined by their MSD (2.7) and (2.8) are both subdiffusions since their MSD is sublinear. There are several ways to model confined diffusion.

On the first hand, we can consider that the domain has non permeable boundaries. Then the motion can be modelled by the SDE (2.1) adding boundary conditions. Equation (2.8) is based on the first two terms of the exact series solution of the MSD of a Brownian particle trapped in a square corral or circular corral defined by the actin cytoskeleton (see [START_REF] Kusumi | Confined lateral diffusion of membrane receptors as studied by single particle tracking (nanovid microscopy). Effects of calciuminduced differentiation in cultured epithelial cells[END_REF][START_REF] Saxton | Lateral diffusion in an archipelago. Single-particle diffusion[END_REF]. The parameter r c quantifies the area of confinement of the trap while its geometry is described by parameters a, b and c. The diffusion process in this restricted area is called confined diffusion.

On the other hand, we can state that the particle is attracted by an external force modelled by a potential well. We can then use the SDE (2.1) with a specific form for the drift : µ i (x) = -∇U i (x)/γ i where -∇U i is the external force of the fluid and γ i is the frictional coefficient. For instance, we may consider the Ornstein-Uhlenbeck process :

(2.10)

dX i t = -λ i (X i t -θ i )dt + σdB 1/2,i t , i = 1, 2
where λ i > 0. Here the particle is assumed to be trapped in a single domain, the potential U is uni-modal and is approximated by a polynomial of order 2 :

U i (x) = (1/2)k i (x i -θ i ) 2 .
The parameter k i = λ i γ i measures the strength of attraction of the potential (related to the potential depth) while θ = (θ 1 , θ 2 ) is the equilibrium position of the particle. The Ornstein-Uhlenbeck process is a confined diffusion according to the MSD criterion since its MSD is sublinear,

(2.11) M SD(t) = 2σ 2 (1 -e -λt ) λ ≤ 2σ 2 t,
here it is written in the case λ i = λ for simplicity.

Anomalous diffusion can occur for two main reasons. First the particle can bind to an immobile trap that can generate long jump times [START_REF] Saxton | Anomalous diffusion due to binding: a Monte Carlo study[END_REF]. In this situation, its motion can be modelled by a continuous time random walk [START_REF] Metzler | The random walk's guide to anomalous diffusion: a fractional dynamics approach[END_REF]. We will not consider this model here as it is not solution of the SDE (2.1). Secondly, the particle can be hindered by mobile or immobile obstacles as the interior environment of cells are crowded with solutes and macromolecules [START_REF] Bressloff | Stochastic models of intracellular transport[END_REF]. Then, a popular model is the fractional Brownian motion [START_REF] Jeon | In vivo anomalous diffusion and weak ergodicity breaking of lipid granules[END_REF]. It corresponds to the case 0 < h < 1/2 and µ i = 0 in (2.1), (2.12)

dX i t = σdB h,i t , i = 1, 2.
Its MSD is given by (2.7) with β = 2h < 1.

2.4. Superdiffusion. At the macroscopic level, the main type of active intracellular transport involves molecular motors which carry particles (cargo) along microtubular filament tracks. The molecular motors and their cargo undergo superdiffusion on a network of microtubules in order to reach a specific area quickly. The molecular motor moves step by step along the microtubules thanks to a mechanicochemical energy transduction process. A single step of the molecular motor is modelled by the so-called Brownian ratchet [START_REF] Reimann | Brownian motors: noisy transport far from equilibrium[END_REF]. When we observe the motion of the molecular motor along a filament on longer time-scales (several steps), its dynamic can be approximated by a Brownian motion with constant drift (also called directed Brownian) (see [START_REF] Peskin | Coordinated hydrolysis explains the mechanical behavior of kinesin[END_REF][START_REF] Elston | A macroscopic description of biomolecular transport[END_REF].

The Brownian motion with drift is solution of the SDE :

(2.13)

dX i t = v i dt + σdB 1/2,i t , i = 1, 2, where v = (v 1 , v 2 ) ∈ R 2
is the constant drift parameter modelling the velocity of the molecular motor. The MSD of the directed Brownian motion is given by (2.9). It is a superdiffusion according to the criterion of MSD. Superdiffusion can also be modelled by the fractional Brownian motion with Hurst parameter 1/2 < h < 1 since its MSD, as we have already said, is given by (2.7) with β = 2h > 1. However, we note that in the biophysic literature the use of the fractional Brownian motion is mainly related to subdiffusion.

3. A statistical test procedure for a single trajectory. We suppose that the trajectory X n = (X t 0 , . . . , X tn ) is generated from some unknown diffusion process (X t ) solution of the SDE (2.1). Our procedure allows to test from which type of diffusion the observed trajectory is generated.

We derive two hypothesis test procedures : one for testing H 0 "(X t ) is a free diffusion" versus H 1 "(X t ) is a subdiffusion", the second for testing H 0 "(X t ) is a free diffusion" versus H 2 "(X t ) is a superdiffusion". Then we aggregate the two procedures to build a three-decision procedure.

3.1. The test statistic. Let us consider the standardized maximal distance T n of the process from its starting point :

(3.1) T n = D n (t n -t 0 )σ 2 n
where D n is the maximal distance of the process from its starting point,

(3.2) D n = max i=1,...,n X t i -X t 0 2
and σn is a consistent estimator of σ. The choice of σ is discussed in Section 3.4. If T n is low, it means the process stays close to its initial position during the period [t 0 , t n ] : it is likely that it is a subdiffusion. On contrary, if T n is large, it means the process goes away from its starting point as a superdiffusion does with high probability. This new measure introduces an order in the diffusion processes solution of the SDE (2.1). Then, it allows to classify them into the different classes of diffusion i-e free diffusion, superdiffusion and subdiffusion. We want to build a test whose null hypothesis is that the trajectory comes from a Brownian motion, the gold standard process in biophysics. As a consequence T n must be a pivotal statistic under the hypothesis H 0 that is the trajectory is Brownian.

Lemma 3.1. Let σn be a consistent estimator of σ such that the distribution of σn /σ does not depend on σ. If (X t ) is a Brownian Motion, the distribution of T n does not depend on σ.

Let q n (α) the quantile of T n of order α ∈ (0, 1) when (X t ) is a Brownian motion. From Lemma 3.1, q n (α) does not depend on σ.

3.2.

Two hypothesis test procedures derived from the test statistic. First we define φ 1,α the hypotheses test associated to H 0 versus H 1 at level α ∈ (0, 1). The procedure φ 1,α is defined through its critical region,

(3.3) R 1,α = {T n < q n (α)} , as the following, φ 1,α (X n ) = 1 if X n ∈ R 1,α 0 otherwise.
Then T n has probability α to lie in the critical region (3.3). According to Lemma 3.1, the level of the test φ 1,α is α,

(3.4) sup σ>0 P 1/2,0,σ (T n < q n (α)) = α.
In a similar way, we can perform the test φ 2,α by replacing subdiffusion by superdiffusion in the alternative hypothesis. The associated critical region is :

(3.5) R 2,α = {T n > q n (1 -α)} .

3.3.

A three-decision test procedure. From the two tests φ 1,α/2 and φ 2,α/2 , we define a new procedure φ as follows,

   we decide H 1 if X n ∈ R 1,α/2 , we decide H 2 if X n ∈ R 2,α/2 , we do not reject H 0 otherwise. (3.6)
This procedure is well defined since the intersection of the critical region R 1,α and R 2,α is empty. This procedure is a three-decision test procedure and admits three kinds of errors, see Table 1.

The first kind of errors is to reject the null hypothesis H 0 while H 0 is actually true. The probability that this error occurs is the level of the test which is defined as, We only control the occurrence of this first kind of error. Then we draw attention that acceptance of H 0 "(X t ) is a free diffusion" does not necessarily demonstrate that H 0 is true. It only means that data does not show any evidence against the null hypothesis. At the end, we reject this assumption in direction to one of the alternatives at level α/2. The second type of errors occurs when we do not reject the null hypothesis while one of the alternatives is true.

(3.7) sup σ>0 E 1/2,0,σ (φ 1,α + φ 2,α ) = α.
The last type of errors is to reject the null hypothesis in favour to a wrong alternative. In the literature of three-decision test such an error is called a Type III error, see for example [START_REF] Rasch | Hypothesis testing and the error of the third kind[END_REF] and references therein.

3.4. Choosing the estimator of σ. Ideally, we would like to find an estimator of σ which is consistent according to the large-sample scheme under the hypotheses H 0 , H 1 and H 2 , and satisfies the assumption that the distribution of σn /σ is free of σ under H 0 . However, the large-sample scheme is not favourable to get an estimator with such properties. For instance, [START_REF] Florens-Zmirou | Approximate discrete-time schemes for statistics of diffusion processes[END_REF] shows that the naive maximum likelihood estimator for the drift parameter has an asymptotic bias of the order of lag time ∆. Then, the high-frequency scheme and the rapidly increasing design turns out to be more convenient to provide consistent estimators. In fact, in the limit, these schemes correspond to the situation in which we have a continuous observation of the process on the time interval of observation. [START_REF] Jiang | A nonparametric approach to the estimation of diffusion processes, with an application to a short-term interest rate model[END_REF] propose non parametric estimators of both the drift and the diffusion coefficient. The consistency of these estimators is proven under the high-frequency scheme only. Therefore, in this section, we discuss about the estimation of the diffusion coefficient under the large-sample asymptotic.

The first proposition to estimate σ may be :

(3.8) σ2 1,n = 1 2n∆ n j=1 X t j -X t j-1 2 2
Even if the estimator (3.8) is strongly consistent under the high-frequency scheme for every process (X t ) solution of (2.1) (Basawa and Prakasa Rao, 1980, Lemma 4.2, p 212), Proposition 1 tells us that it is not the case under the large-sample scheme.

Proposition 1.

• Under H 0 , σ1,n is strongly consistent and the distribution of σ1,n /σ is free of σ.

• If (X t ) is an Ornstein-Uhlenbeck process (2.10), σ2 1,n /σ 2 converges in probability to (1 -e -λ∆ )/(λ∆). • If (X t ) is a Brownian motion with drift (2.13), σ2 1,n /σ 2 converges al- most surely to ∆ v 2 2 /(2σ 2 ) + 1. • If (X t ) is a fractional Brownian motion (2.12), σ2 1,n /σ 2 converges al- most surely to ∆ 2h-1 .
A proof of Proposition 1 is given in Appendix A.1. Proposition 1 states that σ1,n is adequate to our procedure under the null hypothesis. However σ1,n is asymptotically biased under some alternatives. Notice that if (X t ) is an Ornstein-Uhlenbeck process (2.10), then σ2

1,n underestimates σ 2 in average since (1 -e -x )/x < 1 for x > 0. Then T n might be overvalued with this estimator, increasing Type II or type III error rate in our procedure. If (X t ) is a Brownian motion with drift (2.13), σ2

1 overestimates σ 2 in average. Then T n might be overvalued with this estimator, increasing Type II or type III error rate. Similarly, if (X t ) is a fractional Brownian motion (2.12), σ2

1 underestimates σ 2 if h < 1/2, and overestimates σ 2 if h < 1/2.
The second suggestion to estimate σ may be based on the second order differences rather than the first order differences,

(3.9) σ2 2,n = 1 2n∆ n-1 j=1 (X t j+1 -X t j ) -(X t j -X t j-1 ) 2 2 .
As σ2 1,n , σ2 2,n fulfils the assumption of Lemma 3.1 under H 0 . This estimator has the advantage of decreasing the bias under some alternatives. For instance it removes the bias in the case of the Brownian motion with drift.

3.5. Approximation of the distribution of the statistic under the null hypothesis and asymptotic behaviour of our procedure. Theorem 3.1 gives the asymptotic behaviour of our procedure under the null hypothesis.

Theorem 3.1. Let (X t ) be a Brownian Motion on R 2 . Let σn be a consistent estimator of the diffusion parameter σ of (X t ). The test statistic T n converges in distribution to S 0 = sup 0≤s≤1 W s 2 as n → ∞. Here (W t ) is a standard 2D Brownian motion that is the Brownian motion of variance I 2 and initialization W 0 = (0, 0) .

A proof of Theorem 3.1 is given in Appendix A. The limit distribution of the test statistic under H 0 admits an analytical form (see Borodin and Salminen, 1996, Formulae.1.1.4, p. 280):

x ∈ (0, +∞) → ∞ k=1 2e -j 2 0,k /(2x 2 ) j 0,k J 1 (j 0,k ) ,
where x ≥ 0, J ν the Bessel function of order ν and 0 < j ν,1 < j ν,2 < . . . the positive zeros of J ν . Replacing the quantiles q n (α) by the quantiles of S 0 in our test procedure provides us a test of asymptotic level α. Furthermore, Proposition 2 gives the asymptotic behaviour of the test statistic under parametric alternatives when the estimator σ1,n is considered (see Appendix A.2 for a proof). More generally, as long as the estimator σn of the diffusion coefficient is such that σn /σ converges in probability to a positive constant whatever the dynamic of (X t ), then Proposition 2 holds.

Proposition 2. Assume that we consider the estimator (3.8) in our procedure (3.1).

• If (X t ) is an Ornstein-Uhlenbeck process (2.10), then T n converges in probability to 0. • If (X t ) is a fractional Brownian motion (2.12) with 0 < h < 1/2, then T n converges in probability to 0. • If (X t ) is a fractional Brownian motion (2.12) with 1/2 < h < 1, then

T n converges in probability to +∞. • If (X t ) is a Brownian motion with drift (2.13), then T n converges in probability to +∞.

Note that Theorem 3.1 and Proposition 2 allow us to control the error rates of type II and type III under parametric alternatives : the associated error rates converges to 0 with n.

However, as in practice n may be small, the asymptotic approximation of the quantiles of T n may not be accurate. Then the level of the test is no longer α. Since we are able to draw a sample from the distribution of T n under H 0 (see Algorithm 1), we propose a Monte Carlo estimate of the quantile q n (x), 0 < x < 1. This estimate is defined as the [xN ] th order statistic, q 2 shows that there is a significant difference between asymptotic and non asymptotic quantiles. As expected, as n → ∞, q n (α) converges to q(α).

(N ) n (x), of the sample (T (1) n , . . . , T (N ) n ). Table
Input: n, α, L // the length n of the trajectory // the probability α ∈ (0, 1) // the number N of Monte Carlo experiments Result: q (N ) n (α). for i=1 to N do // Simulation of a Brownian trajectory of size n, of variance σ = 1 and with resolution time

∆ = 1. initialization Y (i) 0 = (0, 0) ; for j=1 to n do Draw ∼ N (0 , I2); Y (i) j = Y (i) j-1 + ; end // Computation of the test statistic Compute the ratio T (i) n = D (i) n /σ (i) n from (Y (i) 0 , . . . , Y (i) n ); end Algorithm 1: Simulation of a N -sample (T (1) n , . . . , T (N ) n ) of the distribu- tion of the statistic T n under H 0 . Table 2
Estimation of the quantiles of order α/2 and 1 -α/2 (α = 5%) for different trajectory lengths n, using algorithm 1 with N = 1 000 001. In dealing with a test, we can also be interested in computing the p-value. The p-value of the test H 0 vs H 1 (subdiffusion as the alternative) is defined as :

Estimated quantiles

(3.10)

p 1,n = F n (T n ),
where F n denotes the cumulative distribution function (cdf) of T n under H 0 . The p-value of the test H 0 vs H 2 (superdiffusion as the alternative) is defined as :

(3.11) p 2,n = 1 -F n (T n ).
Testing the hypothesis H 0 vs the hypotheses H 1 or H 2 is more tricky as we use a two-sided test with a non-symmetric distribution. In this case we can define the p-value as :

(3.12) p n = 2 min {p 1,n , p 2,n } .

Doubling the lowest one-tailed p-value can be seen as a correction for carrying out two one-tailed tests.

We estimate F n with the standard empirical distribution function estimated by Monte Carlo simulations using Algorithm 1 :

(3.13) Fn (x) = N -1 N i=1 1(T (i) n ≤ x).
Then we estimate the p-value (3.12) substituting Fn to F n .

4. Multiple test procedure for a collection of trajectories. Trackers compute a collection of particle trajectories from a sequence of images. Then, it is desirable to decide the modes of mobility for a collection of particle trajectories. From now, we consider a collection X m of m trajectories which are simultaneously observed. We denote by X (k) n k the observations associated to the k th particle :

X (k) n k = X (k) t 0 , . . . , X (k) tn k , k = 1, . . . , m X m = X (k) n k , k = 1, . . . , m.
In this section, we denote by P the probability distribution of the m-uplet stochastic processes (X (k) t ), k = 1 . . . m and by E its associated expectation. We assume that the observed trajectories are independent, that means P belongs to the tensorial product of probabilities P, (defined in Section 2) P ∈ P ⊗m . For all trajectories k = 1 . . . m, we derive our trichotomy hypothesis test procedure :

H (k) 0 "(X (k) t ) is a free diffusion" versus H (k) 1 "(X (k) t ) is a subdiffusion" or H (k) 2 "(X (k)
t ) is a superdiffusion". We are faced with the problem of simultaneous tests when the rejections of null hypotheses H (k) 0 are accompanied by claims of the direction of the alternative (H

(k) 1 or H (k)
2 ). In this setup, multiple test procedures are preferable than single test procedures. Indeed, applying the procedure at level α for each trajectory produces in average a number of mα type I errors. A multiple testing procedure aims to control the number of false discoveries. We refer the reader to [START_REF] Shaffer | Multiple hypothesis testing[END_REF]; [START_REF] Roquain | Type I error rate control in multiple testing: a survey with proofs[END_REF]; [START_REF] Grandhi | Multiple Testing Procedures for Complex Structured Hypotheses and Directional Decisions[END_REF] for a review.

A multiple testing procedure of m null hypotheses against two alternative hypotheses is a

rule R 1 (X m ) × R 2 (X m ), where R 1 (X m ) and R 2 (X m ) are disjoint subsets of {H (1) 0 , . . . H (m) 0 }. For i = 1, 2, R i (X m ) is the set of the rejected hypotheses H (k) 0
to the benefit of the alternative H (k) i .

Table 3

Outcomes in testing m null hypotheses against two-alternatives. For i = 1, 2, Ri is the cardinal of Ri(Xm). The variables (Si)i=1...4, (Ti)i=1,2, U, (Vi)i=1,2 are not observed and depend on Xm and P.

True situation Decision

Accept H 0 Accept H 1 Accept H 2 Total H 0 U V 1 V 2 m 0 (P) H 1 T 1 S 1 S 3 m 1 (P) H 2 T 2 S 4 S 2 m 2 (P) Total m -R 1 -R 2 R 1 R 2 m
We may commit three kinds of errors in such a multiple testing procedure. Let us introduce the following notations before listing these errors. For a given P ∈ P ⊗m , we denote by I(P) the subset of indexes {1, . . . m} for which the hypothesis (H (k) 0 ) is actually true and by m 0 (P) the unknown cardinal of the set I(P). We denote by R = R 1 + R 2 the observed number of null hypotheses which are rejected by the multiple testing procedure. Table 3 summaries the number of errors which may occur following a multiple testing procedure.

• We make a type I error on H (k) 0 when we reject H (k) 0 while it is a true null hypothesis. In this case, k belongs to the set

I(P) ∩ (R 1 (X m ) ∪ R 2 (X m )). The number of errors of first kind is V = V 1 + V 2 .
• Type II error occurs when we do not reject a null hypothesis H 0,k while H 0,k is false (k / ∈ I(P)). The number of errors of second kind is

T = T 1 + T 2 .
• The type III errors are directional errors : the index k / ∈ I(P) is correctly rejected (k ∈ R 1 (X m ) ∪ R 2 (X m )), but for the wrong alternative. We mix up the alternatives deciding one while it is the other. The number of errors of third kind is S = S 3 + S 4 .

To measure the type I error rate, it is common to consider the k-familywise error rate (k-FWER) or the false discovery rate (FDR), see [START_REF] Roquain | Type I error rate control in multiple testing: a survey with proofs[END_REF] and references therein. In our settings, controlling the type I error rate is a first step, but it would be necessary to control type III errors as well. In the literature, the sum of the number of errors of first and third kind is controlled using the mixed-directional-family-wise error rate (md-FWER) or the mixed-directional-false discovery rate (mdFDR), see [START_REF] Grandhi | Multiple Testing Procedures for Complex Structured Hypotheses and Directional Decisions[END_REF]. To our knowledge, the mdFWER and mdFDR are only controlled for the problem of testing null hypotheses against two-sided alternatives for finite-dimensional parameters, see for example [START_REF] Guo | On stepwise control of directional errors under independence and some dependence[END_REF] and references therein.

Biologists are interested in the proportions of each dynamic (subdiffusion, superdiffusion and Brownian motion) and their geographic location in the cell. In this context, controlling the FWER, that is the probability to make a single false discovery, is not relevant. That is why we focus on a procedure which enables to control the FDR. [START_REF] Guo | On stepwise control of directional errors under independence and some dependence[END_REF], Section 5) also present several multiple test procedures associated to three-decision problems which aim to control the FDR. Their approach is different since the problem is rewritten as a problem which carries out 3m null hypotheses. Their proposed procedures control strongly the FDR only on 2m null hypotheses among the 3m under the dependence or independence of the test statistics. In this section, we propose to adapt the multiple testing procedures of [START_REF] Benjamini | Controlling the false discovery rate: a practical and powerful approach to multiple testing[END_REF] and [START_REF] Benjamini | On the adaptive control of the false discovery rate in multiple testing with independent statistics[END_REF] controlling the FDR that is the average proportion of false discoveries among the discoveries. We stress that our model is non-parametric. Then we will consider the control of the mdFDR or mdFWER for a next issue.

Let p (k) , p the associated null hypotheses.

The adaptation of the Benjamini-Hochberg (BH) procedure is described in Procedure 1.

Procedure 1 (Adaptation of the Benjamini-Hochberg (BH) procedure).

1. Use the Benjamini-Hochberg procedure on the p-values (p (k) ) k=1...m :

Let k be the largest k for which p

(k:m) ≤ k m α. R α (X m ) is the set of all hypotheses H (k:m) for k = 1, . . . , k . 2. Let R 1,α (X m ) be the subset R α (X m ) such that p (k) 1 < p (k) 2 . 3. Let R 2,α (X m ) be the subset R α (X m ) such that p (k) 1 > p (k) 2 .
The set R α (X m ) is the set of all rejected null hypotheses for our trichotomy test. According to [START_REF] Finner | On the false discovery rate and expected type I errors[END_REF], we have,

∀P ∈ P ⊗m , FDR(R α (X m ), P) = E V max(R, 1) = m 0 (P) m α.
Then the FDR of Procedure 1 is controlled by α. Moreover the p-values p = 1/2) never occurs since such null hypothesis will not be rejected at the step 1 of the Procedure 1.

Actually, we may also use the adaptive BH procedure of [START_REF] Benjamini | On the adaptive control of the false discovery rate in multiple testing with independent statistics[END_REF] as the first step of Procedure 1. Then the Procedure 1 will be referred to as the adaptive (respectively standard) Procedure 1 when we use the adaptive (respectively standard) BH procedure as the first step. The adaptive BH procedure is more powerful than the standard BH procedure. It uses an estimation of the number of true null hypotheses m 0 (P) to increase the power of the BH procedure. [START_REF] Benjamini | On the adaptive control of the false discovery rate in multiple testing with independent statistics[END_REF] simply define the adaptive BH procedure by replacing m by an estimator m0 of m 0 in the BH procedure. The associated FDR is (m 0 / m0 )α and is less than α if m0 ≤ m 0 almost surely. The estimator m0 is based on the plot of the ordered p-values p (i:m) versus i. The ordered p-values corresponding to the true null hypothesis are on the right side of the plot. There exists an index i such that the points (i, p (i:m) )), (i + 1, p (i+1:m) ), . . . , (m, p (m:m) ) can be fit by a straight line with slope κ = 1/(m 0 + 1). Then they use a method to determine i and estimate m 0 from the estimate of κ. The procedure to estimate m 0 presented in [START_REF] Benjamini | On the adaptive control of the false discovery rate in multiple testing with independent statistics[END_REF] is made for m0 to be upward biased. This bias favours the control of the FDR at level α. Due to the fact that m0 does not fulfil the condition m0 ≤ m 0 almost surely, we can not say that the adaptive BH procedure controls the FDR at level α theoretically. However simulations from [START_REF] Benjamini | On the adaptive control of the false discovery rate in multiple testing with independent statistics[END_REF] suggest that the adaptive BH procedure controls the FDR at level α.

5. Simulation study and real data applications. We assess the power of our single test procedure (on a single trajectory) and our multiple test procedure (on a collection of trajectories) by Monte Carlo simulations. We consider parametric alternatives : the Ornstein-Uhlenbeck (2.10) and the fractional Brownian motion 0 < h < 1/2 with Hurst index for subdiffusion processes (H 1 ); the Brownian motion with drift (2.13) and the fractional Brownian motion with Hurst index 1/2 < h < 1 for superdiffusion processes (H 2 ). At the end, we apply our procedure on real data comparing our results with those obtained thanks to a method based on the mean square displacement (MSD). 5.1. Power of the test procedure for a single trajectory. In Section 3, we study the asymptotic distribution of the test statistic under the null hypothesis and parametric alternative hypotheses. More precisely Proposition 2 states that the power of the test under parametric alternatives converges to 1 with n. Figure 3 shows the Monte Carlo estimates of the power un-der the parametric alternatives aforementioned in Proposition 2. For a fixed step of time ∆ and a fixed diffusion coefficient σ, we vary the values of the other parameters and the length n of the trajectories. For each parametric alternatives of Proposition 2, we can use exact simulation schemes.

If (X t ) is an Ornstein-Uhlenbeck process (2.10) which is entered in its stationary regime, then the distribution of the test statistic does not depend on θ (see Appendix A.3). Figure 3(b) shows the plot of the power regarding the values of λ which models the strength of the restoring force toward the equilibrium position θ. Stronger is the force, more powerful is the test.

Furthermore if (X t ) is a Brownian motion with drift with parameters (v, σ) such that v √ ∆ > σ, then the particle goes toward the direction of v while the Brownian random part of the SDE (2.13) does not affect much its trajectory (see Appendix A.3). The bigger is the norm of the drift parameter v, more powerful is the test, see Figure 3(a).

Finally if (X t ) is a fractional Brownian motion, then the distribution of T n depends only on the Hurst index h (see Appendix A.3). Then the test procedure is equivalent to test the null hypothesis "h = 1/2" versus "h = 1/2", see Figure 3(c). 5.2. The Average Power and the mdFDR of the multiple test procedure for a collection of trajectories. The simulation settings are described as follows. According to experience, we choose the number of trajectories to be m = 100 or m = 200. All trajectories are assumed to have the same size n = 30, since this size is reasonable regarding real data. The diffusion coefficient σ and the lag-time ∆ are set to 1. The collection of trajectories X m is composed of :

• m 0 < m Brownian trajectories (H 0 ); • (m -m 0 )/2 subdiffusive trajectories (H 1 ), half from an Ornstein-Uhlenbeck process with parameter λ > 0, half from a fractional Brownian motion with Hurst index 0 < h < 1/2; • (m -m 0 )/2 superdiffusive trajectories (H 2 ), half from a Brownian motion with drift v ∈ R 2 , half from a fractional Brownian motion with Hurst index 1/2 < h < 1.

The parameters to simulate these trajectories are given in Table 4. We take the parameters corresponding to a power of the single test procedure of 80%. Such parameters are used to produce Figure 1 (a). This choice seems coherent in regards to trajectories from real data, see Figure 1 (b). For a given m, the proportion of true null hypotheses H 0 varies : m 0 /m ∈ {0, 0.2, 0.4, 0.6, 0.8}. The mdFDR is a rate which controls the error of type I and type III. It is defined as E((V + S)/ max(R, 1)) (see Table 3). Table 5 shows that the Procedure 1 also controls the mdFDR. The mdFDR and FDR appear to be very close meaning that the number of type III errors is extremely low. Furthermore, the adaptive Procedure 1 (where m 0 is estimated) is less conservative than the standard Procedure 1. As expected, the FDR and mdFDR increase as the proportion of true null hypotheses increases.

Table 4

Parameters used for simulating the alternative hypotheses in the simulation for multiple testing. These choices of parameters correspond to a power of 80% on the power curve of the single test procedure. For simplicity we took σ = 1 for all processes (including Brownian motion). We choose ∆ = 1. To assess the performance of our multiple test procedure, we use the average power [START_REF] Grandhi | Multiple Testing Procedures for Complex Structured Hypotheses and Directional Decisions[END_REF] :

Hypothesis

(5.1)

E S i m i , i = 1, 2
where m i is the number of true alternatives H i and S i (i = 1, 2) is defined in Table 3. In our simulation scheme, we set m i = (m -m 0 )/2. The average power is the expected proportion of hypotheses accepted as H i among all true alternatives H i . Average powers of the different simulations corresponding to different values of m 0 /m and m are shown on Figure 4. First, we can see that the powers of H 1 and H 2 are not very sensitive to the number of hypotheses m for both the standard Procedure 1 and the adaptive Procedure 1. Secondly, the adaptive Procedure 1 is more powerful than the standard Procedure 1 (red and blue dashed lines respectively above red and blue solid lines in Figure 4). The benefit of the adaptive Procedure 1 over the standard Procedure 1 decreases as the proportion of true null hypotheses m 0 /m increases (solid and dashed line of same color getting closer as m 0 /m increases in Figure 4). This is due to the fact that, as m 0 /m tends to 1, m 0 and then m0 tend to m. As a result, the adaptive and standard Procedure 1 become similar.

Remark 5.1. We observe that, given a certain procedure (standard or adaptive Procedure 1), the average power of H 1 is lower than the average power of H 2 (for our particular choice of parametric models for the alternatives H 1 and H 2 ), see Figure 4. It is not due to the choice of parameters as both alternatives H 1 and H 2 are simulated to share the same power (80%) with the single test procedure. Actually, it comes from the fact that the pvalues under H 2 are stochastically smaller than the p-values under H 1 (see Figure 7 Appendix B). Then, the true superdiffusive trajectories are more easily detected as non Brownian in the first step of the (adaptive) Procedure 1 than the true subdiffusive trajectories. We note that, if we use other parametric models for modelling subdiffusion (H 1 ) and superdiffusion (H 2 ), we can have the opposite situation. The main idea is that, as long as the distribution of p-values are not similar under the alternatives H 1 and H 2 , one of the two alternative H 1 or H 2 will always be favoured by the (adaptive) Procedure 1.

Finally, we compare the adaptive Procedure 1 to the MSD classification of [START_REF] Feder | Constrained diffusion or immobile fraction on cell surfaces: a new interpretation[END_REF], based on a fit of the MSD curve to t → t β , see Section 2.1. We assess the two methods on a single collection of trajectories X m with m = 200 and m 0 /m = 0.4, composed of a mixture of Brownian motion, subdiffusion and superdiffusion as described at the beginning of this section. We get the confusion matrices Table 7 and 6 for respectively the adaptive Procedure 1 and the MSD method. As suggested by the limiting curves used by [START_REF] Feder | Constrained diffusion or immobile fraction on cell surfaces: a new interpretation[END_REF] (see Figure 2), the MSD method mixes up the Brownian trajectories with both subdiffusion and superdiffusion (see line 1 of Table 6). Another big issue is that 40% of the particles undergoing subdiffusion are considered as immobile by the MSD method. On the other hand, the adaptive procedure 1 detects well subdiffusion and superdiffusion in the setting of this simulation (line 2 and 3 of Table 7). More importantly, it controls the number of false discoveries through the FDR (line 1 of Table 7).

5.3. Real data : the Rab11a protein sequence. Fluorescence imaging and microscopy has a prominent role in life science and medical research. It consists of detecting specific cellular and intracellular objects of interest at the diffraction limit (200 nm). These objects are first tagged with genetically engineered proteins that emit fluorescence. Then, they can be observed using wide field or confocal microscopy. Several image analysis methods have been developed to quantify intracellular trafficking, including object detection and tracking of fluorescent tags in cells [START_REF] Chenouard | Objective comparison of particle tracking methods[END_REF]; [START_REF] Kervrann | A guided tour of selected image processing and analysis methods for fluorescence and electron microscopy[END_REF]).

Here, we are particularly interested in studying the exocytosis process, that is the mechanism of active transport of proteins out of the cell. Small structures, called the vesicles, travel from organelles to the cell membrane, propelled by motor activity. The vesicle fuses with the plasma membrane and delivers the transported protein in the extra-cellular medium. Given computed trajectories, we investigate here the quantification of vesicles dynamics and trafficking. As explained earlier in the paper, the trajectories can be generally classified into three categories : Brownian motion, subdiffusion and superdiffusion.

As a model of exocytosis/recycling, we focus on the Rab11a protein. This protein is a member of the dynamic architecture of the complex molecular assembly which regulates recycling organelles trafficking. It plays an essential role in the regulation of late steps of vesicle recycling to the plasma membrane, namely the tethering-docking process [START_REF] Schafer | Rab11-FIP2 interaction with MYO5B regulates movement of Rab11a-Containing recycling vesicles[END_REF]). During exocytosis, Rab11a is attached to the vesicle membrane. Then, tracking Rab11a amounts to tracking the vesicle during the exocytosis phase. After the fusion of the vesicle to the cell membrane, Rab11a is recycled in the cytosol. During the recycling step, the tracking of Rab11a is not accurate as the proteins are detached from the vesicle and scatter around the cytosol. It is currently under investigation. For that reason, we focus on the exocytosis process until the fusion time with the cell membrane.

An illustration of the Rab11a sequence is shown in Figure 5 where the dark spots correspond to Rab11a vesicles in a "crossbow" micro-patterned shape cell. A typical image extracted from an image sequence is shown Figure 5. The image sequence is composed of 600 images of size 256 × 240 (1 pixel=160nm) acquired at 10 frames/s (∆ = 0.1s) 1 . We tracked 1 561 trajectories with the multiple hypothesis tracking method with default parameters [START_REF] Chenouard | Multiple hypothesis tracking for cluttered biological image sequences[END_REF]), available on the Icy software (http: www.icy.org). However, we discarded too small and to long trajectories corresponding to tracking errors in most cases. We used a second molecular marker (Transferin Receptor (TfR)) to select trajectories related to the transport of vesicles until the fusion time. The transmembrane TfR protein is fluorescently labeled with a pH-sensitive probe, the pHluorin. Before the fusion time, the pH inside the vesicle is acidic, leading to a very low pHluorin photon emission. When the vesicle fuses to the plasma membrane, the pHluorin gets exposed to the neutral extracellular medium and the fluorescence 1 Live-cell imaging is performed using simultaneous dual color Total Internal Reflection Fluorescence (TIRF) microscopy. All imaging acquisition was carried out in full conditioned medium at 37 o C and 5% CO2 unless otherwise indicated. Simultaneous dual color TIRF microscopy sequences were acquired on a Nikon TE2000 inverted microscope equipped with a x100 TIRF objective (NA=1.49), an Azymuthal TIRF module (Ilas2, Roper Scientifc), an image splitter (DV, Roper Scientific) installed in front of an EMCCD camera (Evolve, Photometrics) that can be bypassed or not, depending on the experimental conditions and a temperature controller. GFP and m-Cherry were excited with a 488 nm and a 561 nm laser, respectively (100mW). The system was driven by the Metamorph software (Molecular Devices).

suddenly increases in the TIRF image plane. This feature allows us to detect the fusion time and the end of the exocytosis process [START_REF] Basset | Detection and estimation of membrane diffusion during exocytosis in TIRFM image sequences[END_REF]). Now the steps to select the trajectories of Rab11a undergoing exocytosis are the following ones :

1. We simultaneously observe two sets of trajectories : TfR and Rab11a trajectories. 2. We match each trajectory of Rab11a with the corresponding trajectory of TfR. 3. We cut the trajectory of Rab11a at the time when the matched trajectory of TfR starts (fusion time).

There is an additional step of selection of trajectories based on mathematical considerations. As we model particles motions with the diffusion processes 2.1, the particles are expected to move over time. Then, we have to get rid of the particles that do not move enough and consequently, can not be modelled by diffusion processes. In practice, we analyse only the trajectories with at least 20 distinct positions and the vesicles that stop at the same position less than K = n/10 times (with n the length of the trajectory). In the case of the aforementioned image sequence, we end up with 166 trajectories whose median length is n = 83.

In Figure 5, our results show that the four procedures -adaptive Procedure 1, standard Procedure 1, single test and MSD method -do not produce similar classification results visually. From the simulations, we found that the MSD method tends to wrongly over-detect subdiffusion and superdiffusion (see Tables 6 and7). This is probably true also in the case of real Rab11a sequence. In Table 8, we give the proportion of each type of diffusion for the different methods aforementioned. The adaptive Procedure 1 tends to decrease the number of Brownian trajectories compared to the standard Procedure 1. It is not surprising as the adaptive Procedure 1 is defined to be more powerful than the standard Procedure 1 : it rejects more easily the null hypothesis. This gain in power benefits to the alternative H 1 (subdiffusion). In fact we detect 23% of subdiffusion for the adaptive Procedure 1 against 16% for the standard Procedure 1 while both detect 4% of superdiffusion (see Table 8). The single test procedure detects even less Brownian motion but we know that it can not control the FDR. In Figure 5, the subdiffusion trajectories labelled with the test approach are more located in the center of the cell in a region corresponding to the Endosomal Recycling Compartment which is known to organize Rab11a carrier vesicles [START_REF] Schafer | Rab11-FIP2 interaction with MYO5B regulates movement of Rab11a-Containing recycling vesicles[END_REF]). It is also true for the subdiffusion trajectories labelled with the MSD analysis but we have just said that there is probably an over-detection of the subdiffusion with this method. We note that we carry the classification of trajectories with our different test procedures and the MSD method on multiple sequences of Rab11a protein, see Figure 6 Appendix B.

6. Discussion. In this paper, we proposed a method for classifying the particle trajectories observed in living cells into three types of diffusion : Brownian motion, subdiffusion and superdiffusion. We used a test approach with the Brownian motion as the null hypothesis. More specifically, we de- veloped a non-parametric three-decision test whose alternatives are subdiffusion and superdiffusion. On the one hand we built a single test procedure for testing a single trajectory, on the other hand we proposed a multiple test procedure for testing a collection of trajectories. These procedures control respectively the type I error and the false discovery rate at level α. It is worth noting that the length of the trajectory n is taken into account in our classification rule. Our approach can be considered as an alternative to the MSD method. It gives more reliable results as confirmed by our Monte Carlo simulations and evaluations on real sequences of images depicting protein dynamics acquired with TIRF or SPT-PALM microscopy.

Then according to the strong law of large numbers, σn converges almost surely to ∆ v 2 /2 + σ 2 .

Ornstein-Uhlenbeck process. Let (X t ) be an Ornstein-Uhlenbeck process (2.10). The SDE (2.10) admits a unique solution (Bressloff, 2014, Section 2

.2.3) (A.3) X t -X s = (X s -θ)(e -λ(t-s) -1) + σ t s e -λ(t-u) dB 1/2 u .
Then (X t ) is a stationary Gaussian process where transition density p(s, x, t, y) is the density of

N x + (x -θ)(e -λ(t-s) -1), σ 2 (1 -e -2λ(t-s) )/(2λ)I 2 .
Then we get that,

E( X t+∆ -X t 2 | X t = x) = x -y 2 p(t, x, t + ∆, y)dy, = x -θ 2 (e -λ∆ -1) 2 + σ 2 (1 -e -2λ∆ )/λ.
Moreover the density µ of the stationary distribution of (X t ) is the Gaussian variable N θ, (σ 2 I d )/(2λ) . Then we obtain that,

E( X t+∆ -X t 2 ) = E( X t+∆ -X t 2 | X t = x)µ(x)dx, = σ 2 (e -λ∆ -1) 2 /λ + σ 2 (1 -e -2λ∆ )/λ, = 2σ 2 (1 -e -λ∆ )/λ.
Now, according to Bibby and Sørensen (1995, Lemma 3.1), if (X t ) is a stationary diffusion, σ2 n converges in probability to E( X t+∆ -X t 2 )/(2∆). We deduce the result.

Fractional Brownian Motion. Let (X t ) be a fractional Brownian motion (2.12). Due to the self-similarity property and the stationary increments of the fractional Brownian motion, the following process,

W (n) t = X t 0 +n∆t -X t 0 (n∆) h σ , t ∈ [0, 1],
is a standard fractional Brownian motion. The statistic associated to the quadratic variation of the process (W (n) t ) may be defined as,

V n = 1 n n i=1 W (n) i/n -W (n) (i-1)/n 2 E W (n) i/n -W (n) (i-1)/n 2 -1, = σ2 n σ 2 ∆ 2h-1 -1.
According to Coeurjolly (2001, Proposition 1), V n converges almost surely to 0. Then we deduce that σ2 n /σ 2 tends to ∆ 2h-1 almost surely.

A.2. Proof of Proposition 2 : the asymptotic behaviour of the test statistic under parametric alternatives. Since the diffusion parameter σ is unknown, the test statistic (3.1) is normalized by an estimator of σ. Proposition 1 states that σn /σ converges in probability to a constant. Therefore, it is sufficient to study the asymptotic behaviour of the test statistic as if σ was known. Then, in this subsection, we consider the test statistic T n as :

(A.4) T n = max i=1,...,n X t i -X t 0 2 σ √ t n -t 0 .
Brownian motion with drift (H 2 ). The process (X t ) is a Brownian motion with drift (2.13) and may be rewritten as,

X tn -X t 0 = v(t n -t 0 ) + σ(B tn -B t 0 ).
Using that (B t ) is a Brownian motion, the distribution of B tn -B t 0 is N (0 2 , (t n -t 0 )I 2 ). Then we have :

(A.5) E X tn -X t 0 σ(t n -t 0 ) - v σ 2 = 2 t n -t 0 .
As t n -t 0 = n∆, we deduce that V n = (X tn -X t 0 )/(σ(t n -t 0 )) converges in probability to v/σ. As the euclidean norm is a continuous function, the variable V n converges in probability to v /σ > 0. Then √ n∆V n converges in probability to +∞. Since T n is lower bounded by √ n∆V n = (X tn -X t 0 ) /(σ √ t n -t 0 ), the proof is complete.

The Ornstein-Uhlenbeck process (H 1 ). The process (X t ) is an Ornstein-Uhlenbeck process (2.10). We assume that the process is in its stationary regime, that means X t 0 is drawn from the stationary distribution that is X t 0 ∼ N (θ, σ 2 /(2λ)I 2 ). The SDE (2.10) admits an unique solution (Bressloff, 2014, Section 2.2.3) (A.6) X t -θ = (X t 0 -θ)e -λ(t-t 0 ) + σ t t 0 e -λ(t-u) dB 1/2 u .

Then we may bound the test statistic T n by,

X t 0 -θ σ √ n∆ + 2 i=1 max k=1...n |X i t k -θ i | σ √ n∆ .
Since X t 0 is drawn from the stationary distribution, the term X t 0 -θ / √ n∆ converges in probability to zero. Now we show that the second term in the previous equation tends to zero in probability as well. We introduce the variables (ξ 1 k , ξ 2 k ) defined as,

ξ i k = (X i t k -θ i ) √ 2λ/σ, k = 1 . . . n, i = 1, 2.
Then for i = 1, 2, the sequence (ξ i k ) k is a standardized stationary normal sequence with covariance function,

r k = E ξ i ξ i +k = e -k∆ , k ≥ -.
Let i be in {1, 2}. Then (a n (max k=1...n (ξ i k ) -b n )) n converges in distribution according to (Leadbetter, Lindgren and Rootzén, 1983, Theorem 4.3.3), where a n = 2 log(n) and b n = a n -(2a n ) -1 (log log(n) + log(4π)). We deduce that max k=1...n (ξ i k )/ √ n∆ converges in probability to 0. Moreover, since (ξ i k ) k is a centred Gaussian process, then max k=1...n (-ξ i k )/ √ n∆ converges in probability to 0 by symmetry. Then we conclude that max k=1...n |X i t kθ i |/ √ n∆ converges in probability to 0.

The fractional Brownian Motion (H 1 ). The process (X t ) is a fractional Brownian motion with h ∈ (0, 1/2). From the property of self-similarity and stationarity of increments of the fractional Brownian motion, the following process, (A.7)

Z (n) t = X tn∆+t 0 -X t 0 σ(n∆) h , t ∈ [0, 1],
is a fractional Brownian motion. We rewrite the test statistic as,

T n = 1 (n∆) 1/2-h max k=1...n Z (n) k/n
Then T n is bounded by, 1 (n∆) ). The process Z (n) has a version with continuous path as a result of being γ-Holder continuous for any γ < h. Let i ∈ {1, 2} be fixed. Then the random variable max k=1...n Z i,(n) k/n is bounded by,

M (n) i = sup t∈[0,1] Z i,(n) t
, which possesses an absolutely continuous density on R * + according to [START_REF] Zaïdi | Smoothness of the law of the supremum of the fractional Brownian motion[END_REF]. That means the sequence max k=1...n Z (n) k/n n is tight. Since h < 1/2, we deduce that T n converges in probability to 0.

The fractional Brownian Motion (H 2 ). The process (X t ) is a fractional Brownian motion with h ∈ (1/2, 1). From the property of self-similarity we get that:

(A.8) Y n = X tn -X t 0 2 2
σ 2 (t -t 0 ) 2h ∼ χ 2 (2).

We observe that T 2 n ≥ Y n (n∆) 2h-1 . Let x be a positive constant. We have : P (T n < x) ≤ P Y n (n∆) 2h-1 < x 2 ≤ P Y n < x 2 /(n∆) 2h-1 . (A.9) Since h > 1/2, x 2 /(n∆) 2h-1 converges to 0 as n → ∞. Then the right hand side of (A.9) converges to 0. That means P (T n < x) converges to 0 as n → ∞ : T n converges to +∞ in probability.

A.3. Dependency of the power on the parameters of the parametric alternatives.

Lemma A.1. Let (X t ) be a Brownian motion with drift (2.13). Let σn be the estimator of the diffusion coefficient defined in Equation (3.8). The distribution of T n (3.1) depends only on the parameter v √ ∆/σ and the trajectory size n.

Proof of Lemma A.1. We may rewrite the strong solution of the SDE (2.13) as,

X t k = X t k-1 + v∆ + σ √ ∆ k , k = 1 . . . n,
where √ ∆ k = B t k -B t k-1 , and (B t ) is a standard Brownian motion. Then ( k ) is a sequence of independent Gaussian variables N (0, 1). Furthermore, we have immediately :

X t k -X t 0 = vk∆ + σ √ ∆ k i=1 i , k = 1 . . . n.
Finally the test statistic T n may be rewritten as,

T n = max k=1,...,n k v √ ∆ σ + k i=1 i 1 2 n i=1 v √ ∆ σ + i 2 .
As the distribution of ( k ) is free of the parameters the distribution of T n depends only on v √ ∆/σ.

Lemma A.2. Let (X t ) be a fractional Brownian motion (2.12). Let σn be the estimator of the diffusion coefficient defined in Equation (3.8). The distribution of T n (3.1) depends only on the parameter h and the trajectory size n.

Proof of Lemma A.2. The fractional Brownian motion may be described by its incremental process [START_REF] Taqqu | Fractional Brownian motion and long-range dependence[END_REF] :

(A.10) k = (X t k -X t k-1 )/(σ∆ h ), k ≥ 1,
where ( k ) is a fractional Gaussian noise which is a stationary standardized Gaussian process with autocovariance function E( k k+i ) = (1/2)(|i + 1| 2h -2|i| 2h + |i -1| 2h ). Finally the test statistic T n may be rewritten as,

T n = max k=1,...,n k i=1 i 1 2 n i=1 i 2
.

Then the distribution of T n depends only on the trajectory size n and on h through the distribution of ( k ).

APPENDIX B: SUPPLEMENTARY FIGURES

Br Sub Sup Br Sub Sup Br Sub Sup Br Sub Sup corresponding to each true alternative hypothesises H1 and H2. The green (respectively orange) line is the threshold h = p (k * ) obtained by the first step of Procedure 1 (respectively Procedure 1) .The null hypothesis is rejected if the p-value is lower than h. The black line is the level α = 5%.

Fig 1 .

 1 Fig 1. Representative trajectories from (a) simulated data, (b) a Rab11a protein sequence in a single cell.(Courtesy of UMR 144 CNRS Institut Curie -PICT IBiSA). For the simulated data in Figure (a), the blue trajectory is Brownian, the red one is from a fractional Brownian motion (parameter h > 1/2), in cyan it is from a Ornstein-Uhlenbeck process and in green from a fractional Brownian motion (parameter h < 1/2). The parameters of the processes are given in Table (4).

  Fig 2.A classification rule for motion modes from MSD. The dashdotted lines are the bound defined by[START_REF] Feder | Constrained diffusion or immobile fraction on cell surfaces: a new interpretation[END_REF], t → t β , β = 0.9 and 1.1. The dashed lines are the pointwise high probability interval of 95% associated to the empirical MSD curve for a standard Brownian motion trajectory of length n = 30. The bounds of the interval are the 2.5% and 97.5% empirical quantile of (2.5) and are computed by Monte Carlo simulation from 10 001 Brownian trajectories of size n = 30.

  the p-value (3.12), (3.10) and (3.11) associated to the k th trajectory, k = 1 . . . m. Let p (1:m) ≤ p (2:m) ≤ . . . ≤ p (m:m) be the ordered p-values, and H

Fig 3 .

 3 Fig 3. Monte Carlo estimate of the power of the test at level α = 0.05 according to the trajectory length n and the parameter associated to the following parametric alternatives : (a) Brownian motion with drift (parameter v = (v1, v2) such that v1 = v2); (b) the Ornstein-Uhlenbeck process (parameter λ) and (c) fractional Brownian motion (parameter h). We use 10 001 Monte Carlo replications for computing each point of the power curves.

Fig 4 .

 4 Fig 4. Monte Carlo estimate of the average power against the proportion of true null hypothesis m0/m in the collection of hypotheses. On the left we test m = 100 hypotheses, on the right m = 200.

Fig 5 .

 5 Fig 5. Map of the classification of the trajectories of the Rab11a sequence with (a) standard multiple test procedure 1, (b) its adaptative version, (c) MSD, (d) single test procedure. The following colour code is used : blue (Brownian motion), red (superdiffusion) and green (subdiffusion).

Fig 6 .Fig 7 .

 67 Fig 6.Boxplots of the proportions of Brownian, subdiffusion and superdiffusion computed from 12 Rab11a sequences. In blue proportions obtained with the single test procedure, in cyan with the Procedure 1, in violet with the adaptive Procedure 1 and in orange with the MSD method. Br stands for Brownian, Sub for subdiffusion and Sup for superdiffusion.

Table 1

 1 The three kinds of error in a three-decision test procedure.

			Decision
	Truth	Do not Reject H 0	Decide H 1	Decide H 2
	H 0 True	No error	Type I	Type I
	H 1 True	Type II	No error	Type III
	H 2 True	Type II	Type III	No error

  Monte Carlo estimate of the FDR and mdFDR for both standard and adaptive Procedure 1 at level α = 0.05. The number of replications is 10 001. The error rate estimations are expressed in percentages.

		Process			Parameter Value
	H1	Ornstein-Uhlenbeck	λ	0.53
	H1	Fractional Brownian	h	0.13
	H2	Brownian motion with drift	v	0.66
	H2	Fractional Brownian	h	0.85
				Table 5
			Standard	Adaptive
	m	m0/m FDR mdFDR FDR mdFDR
	100	0	0	0	0	0.2
		0.2	1	1	3.7	3.7
		0.4	2.1	2.1	4.2	4.2
		0.6	3.2	3.2	4.7	4.7
		0.8	4.1	4.1	4.8	4.8
	200	0	0	0	0	0.4
		0.2	1	1	3.4	3.4
		0.4	2.1	2.1	4	4
		0.6	3.2	3.2	4.6	4.6
		0.8	4	4	4.7	4.7

Table 6

 6 Confusion matrix for the MSD method

	Ground truth/Test label Brownian Subdiffusion Superdiffusion Not moving
	Brownian	19	45	36	0
	Subdiffusion	0	60	0	40
	Superdiffusion	3	0	97	0
	Not moving	0	0	0	0
		Table 7
		Confusion matrix for the adaptive Proc.1
	Ground truth/Test label Brownian Subdiffusion Superdiffusion
	Brownian	96	0	4
	Subdiffusion	23	77	0
	Superdiffusion	10	0	90

Table 8

 8 Percentages of Brownian, superdiffusive and subdiffusive trajectories in the Rab11a sequence according to the different methods of classification.

	Method	Brownian Subdiffusion Superdiffusion
	Standard Proc. 1	80	16	4
	adaptive Proc. 1	73	23	4
	Single test	66	28	6
	MSD	16	63	21

Acknowledgements. We thank Jean Salamero (UMR 144 CNRS-Institut Curie) who provided the microscopy image sequences and for his helpful insights and assistance with experiments.

SUPPLEMENTARY MATERIAL

Supplement A: Source code (https://serpico.rennes.inria.fr/doku.php?id=software:thot:index). matlab code source APPENDIX A: PROOFS Proof of Theorem 3.1. Under the null hypothesis, X t /σ = B t is a standard Brownian Motion. Let us introduce the following random variable, (A.1) Tn = max

Since σn is a consistent estimator of σ and using the Slutsky Lemma, it remains to prove that Tn converges in distribution to S 0 . Using the fact that the increments of the Brownian process are independent and Gaussian, R k is the sum of j independent identically N (0, 1)-distributed random variables. We define the following process,

where x denotes the integer part of x ∈ R. Then we get:

.

Due to Donsker's Theorem (Billingsley, 2013, Theorem 8.2), (W

t ) converges in distribution to the Wiener measure as n → ∞ over the space of continuous function on [0, 1]. Since x → sup t∈[0,1] x(t) is a continuous function on the space of continuous functions from [0, 1] to R, Tn converges in distribution to S 0 .

A.1. Proof of Proposition 1 : the convergence of the estimator (3.8) of the diffusion coefficient. Notice that σn = σ1,n is strongly consistent under the null hypothesis due to the strong law of large numbers and the independence of the increments of the Brownian motion.

We focus now on the three alternatives. According to the alternative, we denote by E the expectation associated to the measure P of the solution of the related SDE ((2.12) or (2.10) or (2.13)).

Brownian with drift. We may rewrite the strong solution of the SDE (2.13) as,

, and (B t ) is a standard Brownian motion. Then the random variables Z k = v∆ + σ √ ∆ k 2 , k = 1 . . . n, are positive independent identically distributed random variables, and admit a moment of order 1, E(Z k ) = ∆ 2 v 2 + 2∆σ 2 .


F Peyraut

J.-L Seichepine

C Coddet

M Hertter

Finite element modeling of abradable materials -Identification of plastic parameters and issues on minimum hardness against coating's thickness

Keywords: Finite element method, optimization, HR15Y indentation test, abradable materials

Abradable materials are used to decrease the gas consumption of aircraft engines by minimizing the gap between the blade tips and the stator. The key idea consists in using the blades themselves to machine the gap on the abradable coating. The best compromise between soft and hard coating properties has to be reached to avoid blades wear and prevent coating erosion by gas flux and particles. The plastic parameters of abradable coating were identified by using an optimization process directly connected to FEA. The first order optimization method (conjugate gradient strategy + golden section algorithm) was applied to achieve the optimal solution. A good agreement was found between experimental and numerical results. The plastic parameters were used to study the hardness variability of abradable materials with the coating thickness. Surprisingly, a minimum hardness value was found while it was expected that hardness should be always decreasing with thickness. It has been demonstrated that this minimum is produced by the boundary conditions influence on hardness measurement.

Introduction

Abradable materials are used to decrease the gas consumption of aircraft engines by minimizing the gap between the blade tips and the stator [START_REF] Schmid | An Overview of Compressor Abradables[END_REF]. These materials are located on the static parts of gas turbines, in front of blades, in order to control the overtip leakage. The key idea consists in machining the gap on an abradable coating with the blades themselves. The best compromise between soft and hard coating properties must be reached to avoid blades wear on the one hand and prevent coating erosion by gas flux and particles on the other.

In the last decades, many attempts have been made to model indentation problems with the finite element method. A detailed bibliographical review of papers published in 1997-2000 and dealing with FEM and BEM simulations of indentation problems can be found in [START_REF] Mackerle | [END_REF]. Most of these works use general purpose FE commercial software like ABAQUS [3], ANSYS [4], MARC [5], CASTEM 2000 [6], NASTRAN [7] or ADINA [START_REF] Olaf | Finite-element analysis of indentation experiments in surface-coated materials[END_REF]. To identify plastic parameters by indentation test, material properties are fitted in such a way that measured and computed hardness match. However, this fitting can be performed by various means. Cao et al. [3] have proposed to solve the inverse problem of identifying the plastic parameters by using an analytical dimensionless function relating the indentation response to the material properties. This approach requires preliminary computational work to fit the coefficients of the dimensionless function. Moreover, the results are only available for a restricted class of materials. To cover a wider range of alloys and pure metals, Zhao et al. [9] have extended the work of Cao and Lu to materials involved in engineering design, such as copper, aluminium, tin and tin alloys. Lee et al. [10] have proposed numerical formulas based on the incremental theory of plasticity for material property evaluation. The material properties used for the FE analysis cover the property range of general metals with a Young's modulus varying from 70 to 400 GPa. It does not include abradable materials for which it was found that calculated Young's moduli are closer to about 10 rather than 70 GPa. Tang and Arnell [11] have proposed a new mathematical model for the effective Young's modulus of the combined substrate and coating and the mean pressure as a measure of indentation hardness. But this model, which provides coating modulus and hardness independent on substrate effects, is only applicable for elastic-perfectly plastic materials. In [12], FEM is employed to obtain indentation parameters like maximum and residual penetration depths. These parameters are used in the indentation force equation which is solved by the linear regression method. However, this method only concerns materials exhibiting exponential hardening and can not be extended to others without loss of accuracy. More recently, Beghini et al. [13] have proposed an inverse procedure, starting from experimentally measured pairs of load and depth of penetration, for deducing the material parameters by means of an optimization algorithm. This procedure requires preliminary FEA to build up a database of numerical curves relating the load to the depth of penetration. Moreover, it can only be applied to plastic behaviour involving three parameters such as the Hollomon like power law or a linear strain hardening.

The limited range of applicability is generally considered as the main drawback of methods using response surface even if most classic alloys and metals are accounted. To provide a general approach applicable to any kind of material, including the abradables, an alternative is to directly connect optimization algorithms to FE analysis. Bouzakis et al. [14] have proposed this kind of approach by using FEM to extract material parameters for a multilinear plastic law. The optimization algorithm was developed by using minimization technique supported by the FEM package ANSYS. The material parameters were updated and FE computation repeated until a criterion between measured and calculated penetration forces will be fulfilled. We have selected a similar approach by using the FE code ANSYS [15]. The difference between measured and computed data was minimized to find the appropriate material properties satisfying the similarity of computational and experimental hardness. The first order optimization method (conjugate gradient strategy + golden section algorithm) was applied to achieve the optimal solution. This method is time consuming but converges accurately towards the global minimum of the objective function. Two indentation depths corresponding to two different coating thicknesses have been considered instead of the complete load-depth indentation curve related to a single coating thickness. Since two different coating thicknesses are considered for the reverse analysis, the proposed approach likely works as the dual method [16] used for sharp indentation to provide a single combination of appropriate material parameters. This ensures that the relation between material parameters and hardness is one-to-one. The proposed method can be directly applied to any kind of bilinear plastic material and does not need any preliminary computational work to fit coefficients and material constants as required by other methods.

Thanks to the identification of the plastic parameters, an extensive study of hardness against thickness was performed by using an axisymmetrical FE model of the HR15Y indentation test. Several hardness tests are available to perform hardness measurement, each one providing different hardness value (Vickers, Rockwell and Rockwell superficial, Brinell . . . ). Among these various hardness tests, the HR15Y indentation test [17] is the most widely used by industry manufacturers to control the quality of thermally sprayed abradable coating. The general principle is the measurement of the permanent increase in depth of penetration with a 12.70 mm diameter steel ball indenter, a preliminary test force of 29 N and a total test force of 147 N. The FE results show a fairly good agreement with measured hardness. However, a minimum hardness value was found while it was expected that hardness should be always decreasing with thickness. By focussing the study on the central part of the coating, which is the most significant region to account for plastic effects, we have demonstrated that this minimum hardness is produced by the boundary conditions influence on hardness measurement.

Problem presentation

Abradable coatings are made of complex 3-D heterogeneous composite materials randomly distributed (Fig. 1). A bilinear plastic law was used to model the mechanical behaviour of these coatings (Fig. 2). To take into account the scale of the heterogeneities (micron) and the scale of the coating (millimetre), the elastic properties were computed thanks to a multi-scaled approach [18]. The Young's modulus and the Poisson ratio calculated with this approach for the abradable coating considered in this work are respectively 15 GPa and 0.2.

The plastic parameters were identified by minimizing the difference between computed and measured hardness's. The general principle is the measurement of the permanent increase in depth of penetration under a multi-step load. Various geometrical indenters can be used as spheroconical tip or ball with various diameters. In [19], the effects of different ball materials on HRB hardness measurement values are discussed. In our work, the standardized HR15Y test (Fig. 3) is selected because it is the most widely used by manufacturers of abradable coatings for aero-engines application. This test uses a 12.70 mm diameter steel ball indenter, a preliminary test force of 29 N and a total test force of 147 N (Fig. 4). The preliminary test force is applied to cause an initial penetration h i of the specimen in accordance with recommendations provided by hardness measurement standards [17]. The total force is next applied to provide the maximum depth of penetration h m . The permanent depth of penetration h p is obtained by unloading the system from 147 N to 29 N. By introducing the difference e between the permanent h p and the initial h i depths of penetration with a μm scale, the material hardness H is defined as H = 100-e.

Numerical procedures

Finite element model

The FE code ANSYS was used to build a bi-dimensional model of the coating and the indenter (Fig. 5). Both of them are meshed in the (x, y) plane by four nodes PLANE42 elements used with axisymmetrical option. The conventional notations used by ANSYS were adopted to represent the radius and the rotation axis (respectively x and y directions). The FE model is made of 2176 nodes and 2163 elements for the lowest thickness used in the present work (0.1 mm) and 14893 nodes and 14723 elements for the highest (20 mm). This model includes local non linearity to model frictionless contact between the indenter and the coating and material non linearity to represent the plastic behaviour of abradable materials. Frictionless contact is assumed because it has been observed that friction produces negligible effects on the curve relating the load to the depth of penetration [10]. TARGE169 and CONTA171 elements were used to model contact between the coating and the indenter and the mesh density is refined in the vicinity of the contact area to account for high stress gradient growing in this region. TARGE169 element represents the "target" surface and CONTA171 is a "surface-to-surface" contact element used to represent contact and sliding between 2-D "target" surfaces and a deformable surface defined by this element. The bottom of the coating was constrained against displacement in the y-direction and the vertical symmetry axis, defined by x = 0, was constrained against displacement in the x-direction to model axisymetrical condition. To model the plastic behaviour of the coating, a bilinear isotropic law was used (Fig. 2). The material parameters are the yield stress σ y , the tangent modulus E T and the Young modulus E. A preliminary FEA is performed on a 2D model of the material microstructure in order to compute E [18]. σ y and E T are next identified by using a reverse procedure described in Section 3.2.

Optimization

The optimization process used to identify the yield stress σ y and the tangent modulus E T values is presented in this subsection. The idea is to determine σ y and E T in order to minimize the difference D between measured and computed hardness Min

σ y , ET [D]. (1)
The yield stress σ y and the tangent modulus E T represent the design variables and the objective function D is defined by the distance between experimental and numerical hardnesses. This distance is expressed in terms of percentage by using the usual Euclidian norm

D = 100 × H c -H m H m . ( 2
)
Measured and computed hardnesses are stored with two vectors H m and H c

H m = 24.1 10 ; H c = H c 1.2 H c 2.3 . ( 3
)
Each component corresponds to a different thickness. Two components are needed since two plastic parameters have to be found. In order to extract the most relevant information from experimental data's, the lowest and highest coating thicknesses were logically considered (1.2 and 2.3 mm). The first order optimization method was applied to achieve the optimal solution. This method combines a conjugate gradient strategy to determine the descent direction and a golden section algorithm to perform a line search along the descent direction. It is noted in Figure 6 that nearly twenty hours are required to achieve the optimal solution on a PC (Pentium 4/2.8 GHz). The optimization procedure is thus computationally demanding. However, it is possible to reduce the computation time by about a half (11 h instead of 19 h) without losing accuracy (D = 0.86 instead of 0.55%). Moreover, the method is efficient and provides accurate results as it will be demonstrated in the next section.

Computational results

Plastic parameters

The objective function surface is depicted in Figures 7 and 8. This surface was built by computing hardness against design variables (σ y , E T ) satisfying 0.5 Mpa σ y 10 Mpa and 1% E T /E 50%. The iterative solutions provided by the optimization procedure are represented by squares. It is observed that the procedure converges perfectly well towards a minimum defined by: σ y = 3 MPa and E T = 1.5% E. It is also noticed that the first iterations move from the starting point along a straight line defined by σ y = 5 MPa. The objective function is thus more sensitive to the tangent modulus than the yield stress. The most significant effect is found for E T changing from 0 to 10% E corresponding to the stiffest slope of the objective function. This observation is in agreement with the results obtained by Gan et al. in [20]. It is also noticed that there is only one global minimum in the considered range of parameters variation. Non-uniqueness, i.e. identical hardness's obtained from different combination of plastic parameters [21][22][23][24], is not observed. Since two thicknesses are used to perform the reverse analysis, the proposed approach likely works as the dual method [16] used for sharp indentation to provide a single combination of appropriate material parameters. It is actually noticed that two different thicknesses provide two different load-indentation depth curves (Fig. 9). This ensures that the relation between material parameters and hardness is one-to-one.

Hardness against thickness

The plastic parameters identified in Section 4.1 and the FE model described in Section 3.1 were used to compute hardness against coating thickness. Comparison between measured and computed hardness is presented in Figure 10. A perfect agreement is noticed if the coating thickness is equal to 1.2 or 2.3 mm. This result was expected because these thicknesses were used to equate to zero the difference between tested and computed hardness (see Sect. 3.2). A very good agreement between numerical and experimental results is also observed for other thicknesses since the difference never exceeds 6.7% (Tab. 1). If the thickness is less than 3 mm, it is noted that hardness decreases sharply with the thickness up to a minimum value (Fig. 11). On the contrary, if the thickness is more than 3 mm, the hardness firstly increases and then reaches a thickness independent value. The variability of hardness against thickness has been already observed for real material and is not induced by a numerical artefact [25].

Relationship between hardness measurements and residual stresses in the coating

It was observed in Figure 11 that the hardness reaches a minimum for a thickness value equal to 3 mm. By investigating the relationship between hardness and stress field, we will demonstrate in this section that the minimum hardness is produced by the boundary conditions influence on hardness measurement. It is firstly noticed that the central part of the coating is the most significant region to account for plastic effects. The Von Mises plastic strain is actually concentrated within the range x ∈ [0-2.54 mm] (Fig. 12). In this region, the plastic strain extends radially from the contact zone to the coating bottom. This general trend of radial expanding has been already observed in [5,19,26]. By focusing on the central part of the coating, it is noted that the coating thickness giving the maximum reaction force is equal to 3 mm (Fig. 13). This corresponds to the thickness value providing the minimum hardness (Fig. 11). Since the reaction force results from the constrained displacement applied to the coating bottom, this demonstrates that hardness variability, and particularly occurrence of minimum hardness, are ascribed to the boundary condition influence on hardness measurement. The distribution of the vertical stress component σ yy on the coating bottom was also studied to understand why the maximum force is related to a 3 mm coating thickness. Three stress curves are plotted against x coordinate in Figure 14. These three stress curves are related to a 3 mm thickness and to two thicknesses close to 3 mm (2.5 and 3.5 mm). These curves intersect each other nearly x equal to 1.6 mm. The lowest thickness provides the lowest stress value up to the intersection while the opposite trend is observed beyond the intersection. A similar conclusion was found for the reaction forces acting within the range x ∈ [0,1.6 mm] and x ∈ [1.6, 2.54 mm] (Fig. 15). The maximum force then results from the competition between two opposite trends and depends on the stress distribution on the coating bottom. The minimal hardness observed in this work is thus produced by the boundary conditions influence on hardness measurement.

Conclusion

A finite element model of the HR15Y hardness test was developed to identify the plastic parameters of abradable materials. To provide a general approach applicable to any kind of material, including the abradables, the parameters identification was performed by connecting optimization algorithms to FE analysis. The first order optimization method was used in this work. This method is time consuming but has been found accurate and ef- ficient to provide the global minimum of the objective function. In the considered range of parameters variation, non-uniqueness has not been observed (i.e. identical hardness's obtained from different combination of plastic parameters). Since two different load-indentation depth curves corresponding to two different thicknesses are considered for the reverse analysis, the proposed approach likely works as the dual method used for sharp indentation to provide a single combination of appropriate material parameters.

The plastic parameters were used to study the variation of hardness with the coating thickness. The results show a fairly good agreement between computed and measured hardness. However, a minimum hardness value was found while it was expected that hardness should be always decreasing with thickness. By focussing the study on the central part of the coating, which is the most significant region to account for plastic effects, we have demonstrated that this minimum hardness is produced by the boundary conditions influence on hardness measurement.

The good agreement found between numerical and experimental results demonstrates that a bilinear hardening model is sufficient to simulate the loading cycle of the HR15Y hardness test applied to the abradable layers studied in this work. This simple model may be proposed as a starting assumption for the simulation of other experiments, to be validated afterwards. However, in order to investigate more accurately a widest range of materials and loading conditions, the proposed method could be easily extended to multilinear plastic materials. To include such materials in the identification process, the number of design variables used for the optimization stage should be increased and a multilinear plastic law used rather than a bilinear one for the finite element analysis.

Fig. 1 .

 1 Fig. 1. Micrograph of a typical abradable coating structure; left: surface view -right: cut view.

Fig. 2 .

 2 Fig. 2. Elastic-plastic behaviour law.

Fig. 3 .

 3 Fig. 3. Schematic of HR15Y test rig.

Fig. 4 .

 4 Fig. 4. Multi-step loading.

Fig. 5 .

 5 Fig. 5. Axisymmetrical finite element model.

Fig. 6 .

 6 Fig. 6. Computation time.

Fig. 7 .

 7 Fig. 7. Objective function versus design variables.

Fig. 8 .

 8 Fig. 8. Zoom on the minimum of the objective function.

Fig. 9 .

 9 Fig. 9. Numerical load-indentation depth curves.

Fig. 10 .

 10 Fig. 10. Comparison between measured and computed hardness.

Fig. 11 .

 11 Fig. 11. Computed hardness.

Fig. 13 .

 13 Fig. 13. Reaction force applied to the coating bottom within the range x ∈ [0,2.54] .

Fig. 14 .

 14 Fig. 14. σyy stress component distributed on the central part of the bottom coating.

Fig. 15 .

 15 Fig. 15. Reaction force applied to the coating bottom within the range: (1) x ∈ [0, 1.6] -(2) x ∈ [1.6, 2.54].

Table 1 .

 1 Comparison between measured and computed hardness.

	Coating thickness [mm] Measured hardness Computed hardness Difference (%)
	1,2	24,1	24	0,4
	1,5	19	17,72	6,7
	1,55	16	16,85	-5,3
	1,6	17	16,2	4,7
	1,8	14,5	13,62	6,1
	2,3	9,5	9,6	-1,1

Fig. 12. Von Mises plastic strain component for various thicknesses.

Acknowledgements

This research work was completed from May 2002 to May 2006 within the European Seal-Coat project funded by the European Commission under the FP5 Growth Program. The Seal-Coat project has involved UTBM (Peter Chandler, in charge of the Project Management, and the authors of this paper) and the following partners: ESIL, Dublin/IE (Barry O'Reilly) Rolls-Royce plc, Derby/GB (Noel Hopkins and Chris Sellars) MTU Aero Engines GmbH, München/DE (Manuel Hertter and Ulrike Hain) Euromat GmbH, Heinsberg/DE (Ulrich Morkramer) RWTH Aachen University, Aachen/DE (Jochen Zwick) Neomet Limited, Manchester/GB (Dieter Sporer) Institute of Plasma Physics, Praha/CZ (Jiri Matejicek) RALSA, Langreo/ES (Ramiro Rilla)

This research work was completed within the Seal-Coat project funded by the European Commission under the FP5 Growth Program.