Vincent Briane

Antoine Salomon

Myriam Vimond

Charles Kervrann
email: charles.kervrann@inria.fr

A computational approach for detecting microdomains and confinement domains in cells: a simulation study

In this paper, we aim at detecting trapping areas (equivalently microdomains or confinement areas) within cells, corresponding to regions where molecules are trapped and thereby undergo subdiffusion. We propose an original computational approach that takes as input a set of molecule trajectories estimated by appropriate tracking methods. The core of the algorithm is based on a combination of clustering algorithms with trajectory classification procedures able to distinguish subdiffusion, superdiffusion and Brownian motion. The idea is to automatically identify trapping areas where we observe a high concentration of subdiffusive particles. We evaluate our proof of concept on artificial sequences obtained with a biophysicsbased simulator (Fluosim), and we illustrate its potential on real TIRF microscopy data.

Introduction

As the interior of a living cell is a fluctuating environment, the trajectories of intracellular particles can be modeled with stochastic processes having continuous paths. Diffusions belong to this class of processes and can model a large range of intracellular movements. They are widely used in the biophysical literature [Qian et al., 1991, Saxton and[START_REF] Saxton | Single-particle tracking: applications to membrane dynamics[END_REF]. Biophysicists distinguish four main types of diffusion, namely Brownian motion (also referred to as free diffusion), superdiffusion, confined diffusion and anomalous diffusion. These different diffusion types correspond to specific biological scenarios as briefly described in [Briane et al., 2019a]. A particle evolving freely inside the cytosol or along the plasma membrane is modeled by free diffusion. Its motion is due to constant collisions with smaller particles undergoing thermal fluctuations. Thus, 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 suggests that particle motion is direct. The particles (called in this context cargos) are carried by molecular motors along microtubules and actin filament networks. Superdiffusion models the motion of molecular motors and their cargos. Confined or restricted diffusion [START_REF] Metzler | The random walk's guide to anomalous diffusion: a fractional dynamics approach[END_REF]Klafter, 2000, Hoze et al., 2012] 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 particles that encounters dynamics or fixed obstacles [Saxton, 1994, Berry and[START_REF] Berry | Anomalous diffusion due to hindering by mobile obstacles undergoing brownian motion or ornstein-ulhenbeck processes[END_REF], or particles slowed down by the contrary current due to the viscoelastic properties of the cytoplasm. Very recently, [START_REF] Smith | Macromolecular crowding directs the motion of small molecules inside cells[END_REF] studied the effect of heterogeneous crowding on the motion of intracellular particles. They assumed that the motion of the particle is hindered by spherical objects, and then derived an advection-diffusion equation to describe the motion of a small particle traveling across an arbitrary distribution of large crowder molecules. Finally, we emphasize that particles can switch from one motion type to another over time; for instance, [START_REF] Fusco | Single mRNA molecules demonstrate probabilistic movement in living mammalian cells[END_REF] showed experimentally that mRNA molecules undergo cycles between anchoring, free diffusion and superdiffusion (active transport).

An amount of work has been dedicated to the description of motions of individual particles in the cell. The second step consists in understanding the spatial distribution of these motions within the cell and interpreting the underlying biological mechanisms. As a pioneer work, [START_REF] Hummer | Position-dependent diffusion coefficients and free energies from bayesian analysis of equilibrium and replica molecular dynamics simulations[END_REF] studied the spatial mapping of diffusion coefficients using Bayesian inference. This paper inspired other mapping methods (e.g. see [START_REF] Masson | Mapping the energy and diffusion landscapes of membrane proteins at the cell surface using high-density single-molecule imaging and bayesian inference: application to the multiscale dynamics of glycine receptors in the neuronal membrane[END_REF] and [START_REF] Hoze | Heterogeneity of ampa receptor trafficking and molecular interactions revealed by superresolution analysis of live cell imaging[END_REF]). In this paper, we are interested in localizing areas where subdiffusion occurs as in [START_REF] Hoze | Heterogeneity of ampa receptor trafficking and molecular interactions revealed by superresolution analysis of live cell imaging[END_REF]. Regulation of cellular physiological processes relies on molecular interactions (binding and unbinding) at specific places and involves trafficking in confined local microdomains. Our objective is to automatically detect these microdomains from groups of particle trajectories. [START_REF] Hoze | Heterogeneity of ampa receptor trafficking and molecular interactions revealed by superresolution analysis of live cell imaging[END_REF] model these microdomains as potential wells that attract intra-cellular particles. The authors use an Eulerian method based on a non-parametric estimation of the drift parameter of the underlying diffusion process. This method assumes a high concentration of particles over the spatial domain of interest to be meaningful. Moreover, this Eulerian approach is not designed to handle a mixture of different diffusion types observed in a small region. All the trajectories in a small region are uniformly averaged, potentially leading to wrong conclusions. To overcome this problem, we propose an original procedure to detect microdomains.

Our method is based on the combination of the clustering algorithm Dbscan [START_REF] Ester | A density-based algorithm for discovering clusters in large spatial databases with noise[END_REF] and a three-label classification procedure [START_REF] Briane | Statistical analysis of particle trajectories in living cells[END_REF] which assigns a label (Brownian motion, subdiffusion and superdiffusion) to each trajectory. Our results on simulations mainly demonstrate that:

(i) our method detects satisfyingly the microdomain shapes;

(ii) our method is able to handle a mixture of several diffusion types.

The remainder of this paper is organized as follows. In Section 2 we present the mathematical models and the notations used in the paper. In Section 3, we present related works and focus on the Eulerian method proposed by [START_REF] Hoze | Heterogeneity of ampa receptor trafficking and molecular interactions revealed by superresolution analysis of live cell imaging[END_REF]; we emphasize the underlying assumptions, the key parameters, and the properties of this approach used to detect potential wells within cells. In Section 4, we describe the simulation framework corresponding to the Fluosim simulator [START_REF] Lagardère | FLUOSIM: Simulator of single molecule dynamics for fluorescence live-cell and super-resolution imaging[END_REF] used in our experiments. In Section 5, we describe the procedure that combines the Dbscan clustering algorithm [START_REF] Ester | A density-based algorithm for discovering clusters in large spatial databases with noise[END_REF] with the test procedure described in [START_REF] Briane | Statistical analysis of particle trajectories in living cells[END_REF]. In Section 6, we evaluate the performance of the computational method on simulated sequences, and we illustrate its potential on real TIRF microscopy image data depicting small zones with a high concentration of subdiffusive vesicles. In Section 7, we discuss the main features of our approach.

Notations and mathematical framework

In the sequel, we assume a population of m independent trajectories, with m of order of magnitude of 10 3 particles. We denote X (i) t ∈ D the position of the i th particle at time t, where D ⊂ R 2 denotes a spatial bounded domain. Let us define P subregions S 1 , . . . , S P such that S p ⊂ D in which the particles can be trapped (see Fig. 1(a)). The union of subdomains is denoted S = ∪ P p=1 S p , and S = D \ S is defined as the complementary domain where the particles undergo a Brownian motion with diffusion coefficient σ S . We impose that the particles are normally reflected at the boundaries of D = S ∪ S.

Inside a subregion S p ⊂ D, a particle can switch between two states: trapped or non-trapped (also denoted as free):

(i) the particle is trapped in S p : the particle motion is a confined Brownian motion normally reflected at the boundary ∂S p of S p , with a diffusion coefficient σ T ≤ σ S ;

(ii) the particle is not trapped in S p : the particle motion is a Brownian motion with diffusion coefficient σ S .

Formally, we introduce the following indicator variable:

φ(X (i) t ) = 1 if the particle i is trapped at time t, 0 if the particle i is not trapped at time t. (1)
Given the binding rate k + and the unbinding rate k -, the probabilities of a particle belonging to S to be trapped are defined as follows:

P (φ(X t+h ) = 1| φ(X t ) = 0, X t+h ∈ S) = k + h + o(h), (2)
P (φ(X t+h ) = 0| φ(X t ) = 1, X t+h ∈ S) = k -h + o(h),
where h > 0 and o(h)/h → 0 as h → 0. In the sequel, the binding and unbinding rates k + and k -are two constants which do not depend on space. Hence, a trapped particle close to the boundaries ∂S p has the same probability to unbind to the trapping region S p than a trapped particle in the middle of the region S p . For the sake of simplicity, we assume that all the subregions S p share the same common parameters k + , k -, σ T and σ S . In our estimation framework, we observe a collection of m two-dimensional trajectories X m . We denote X (k) n i the vector of n i spatial positions along the trajectory associated to the i th particle:

X (i) n i = X (i) t 0 , . . . , X (i)
tn i , i = 1, . . . , m, (3)
X m = X (i) n i , i = 1, .
. . , m . We assume that each discrete trajectory is generated by a stochastic process (X (i) t ) with continuous path defined on the spatial domain D ⊂ R 2 and which is a solution of the stochastic differential equation (SDE),

dX (i) t = µ(X (i) t )dt + σ (i) dB (i),h (i) t , t ∈ [t 0 , t n i ], (4)
where

B (i),h t
is a two-dimensional fractional Brownian motion of unknown Hurst parameter h (i) ∈ (0, 1), σ (i) > 0 is the diffusion coefficient, and µ (i) : D → R 2 is the drift term. The different trajectories do not need to have the same drift parameter and the same diffusion coefficients, provided that the drift functions and diffusion coefficients do not depend on time.

Most related work

In this section, we focus on the two most related methods developed to detect 2D microdomains. First, we describe the Eulerian method proposed by [START_REF] Hoze | Heterogeneity of ampa receptor trafficking and molecular interactions revealed by superresolution analysis of live cell imaging[END_REF], based on the two following steps:

(i) non-parametric estimation of the drift vector field from multiple trajectories;

(ii) fitting to a parametric drift function.

As starting point, [START_REF] Hoze | Heterogeneity of ampa receptor trafficking and molecular interactions revealed by superresolution analysis of live cell imaging[END_REF] model microdomains as potential wells that attract intra-cellular particles. A high density of particles within a local neighborhood is assumed. Formally, let us consider a collection of m independent trajectories X m such that all the trajectories X (i)

n i = X (i) t 0 , . . . , X (i) tn i
are generated from the common diffusion process (X t ) solution of the overdamped Langevin equation defined as:

dX t = -∇U (X t )dt + σdB t , ( 5
)
where ∇ is the gradient operator, U : R 2 → R is the potential function, and σ > 0 the diffusion coefficient. In the sequel, the two-dimensional drift function (or vector field) is defined as µ(x) = -∇U (x). The the drift function is defined through the limit as:

µ(x) = lim ∆t→0 1 ∆t E(X t+∆t -X t |X t = x). (6)
In [START_REF] Hoze | Heterogeneity of ampa receptor trafficking and molecular interactions revealed by superresolution analysis of live cell imaging[END_REF], an empirical estimator μ(x) of the drift is computed over a square r × r window W r (x) centered at point x ∈ D ⊂ R 2 , in order to reflect the local conditioning X t = x in Equation ( 6). The estimation of μ(x) amounts to computing the average of all displacements of particles passing through the window W r (x). The estimation of μ(x) is spatially varying and its accuracy depends on the number of trajectory points in W r (x). Fortunately, we generally observe a high density of particles in trapping areas. In these microdomains, the drift vector field is thus reliably estimated but it is assumed that, locally in the r × r neighborhood, all the trajectories undergo the same diffusion process with the same drift function. The parameter r is a free parameter of the algorithm, which can be interpreted as a bandwidth parameter in density estimation. In Fig. 2, we illustrate the estimation of drift from trajectories on a artificial sequence depicting two circular trapping areas. A general view of trapping areas with variable sizes generated with Fluosim [START_REF] Lagardère | FLUOSIM: Simulator of single molecule dynamics for fluorescence live-cell and super-resolution imaging[END_REF], is shown in Fig. 1(b)).

In the second step, the authors assume that the potential U (x) has a parametric form and can be approximated by a truncated polynomial of second order, defined as:

U (x) = A( x -x 0 2 a,b -1) if x -x 0 2 a,b < 1, 0 otherwise, (7)
where x 0 is the attractor, and A > 0 is the depth of the potential modeling the strength of the attractive force toward x 0 , restricted to the underlying 2D ellipse. The parameters (a, b) are the axis lengths of the ellipse of center x 0 in which the particle are under the action of an attractive force. In the computational scheme, the parameters a, b, x 0 and A are assumed partially known and/or estimated by minimizing a least-square error. To conclude, the approach of [START_REF] Hoze | Heterogeneity of ampa receptor trafficking and molecular interactions revealed by superresolution analysis of live cell imaging[END_REF] is very appropriate to analyze a large number of tracks as computed in single-particle tracking PALM imaging. The motivation was to compute the residence time of receptors in dendritic spines using short receptor trajectories. The main advantage is the processing of very short tracks. Nevertheless, it cannot handle mixtures or changes of motion occurring in a local neighborhood. Another Eulerian method was described in [START_REF] Masson | Mapping the energy and diffusion landscapes of membrane proteins at the cell surface using high-density single-molecule imaging and bayesian inference: application to the multiscale dynamics of glycine receptors in the neuronal membrane[END_REF]. The authors propose a Bayesian framework to infer the drift vector field and diffusion coefficients. The idea is first to decompose the image into non overlapping blocks. All the subtrajectories inside a block (or mesh element) are supposed to be driven by the same 1). We compute μ(x) on each point of this 200 × 200 lattice from trajectories (a) (r = 5 pixels). We can clearly see the two microdomains S 1 and S 2 (b): the vector field converges to the centers of these domains delineated by two red circles (b). We can see that the vector field is orthogonal to the boundaries of the square D. It is due the normal boundary condition on ∂D. Elsewhere, the μ(x) is almost null, reflecting Brownian motion.

stochastic differential equations (SDE). Instead of using non-parametric estimates of drift and diffusion coefficients as in [START_REF] Hoze | Heterogeneity of ampa receptor trafficking and molecular interactions revealed by superresolution analysis of live cell imaging[END_REF], the authors consider a Gaussian approximation of the SDE to model the likelihood function. Meanwhile, a Jeffreys prior is used as a default prior distribution for the drift and the diffusion coefficient, supposed to be constant in each block. It follows that the a-posteriori distribution on the whole space is defined as the product of the a-posteriori distributions of each block. An alternative prior is proposed to penalize gradients of the drift field (or the coefficient diffusion field) [START_REF] Beheiry | A primer on the bayesian approach to high-density single-molecule trajectories analysis[END_REF]; in that case, the maximum a-posteriori (MAP) estimator is more computationally demanding. Finally, once the local drifts (and diffusion coefficients) are estimated, [START_REF] Masson | Mapping the energy and diffusion landscapes of membrane proteins at the cell surface using high-density single-molecule imaging and bayesian inference: application to the multiscale dynamics of glycine receptors in the neuronal membrane[END_REF] propose to fit the estimated vector field to a parametric vector field as in [START_REF] Hoze | Heterogeneity of ampa receptor trafficking and molecular interactions revealed by superresolution analysis of live cell imaging[END_REF]. They consider a penalized least square estimator of the vector field to impose a smooth solution. The main advantage of the Bayesian method is that prior information about the drift and diffusion coefficient can be easily incorporated, including errors due to misconnections of particles during tracking.

In this paper, we are interested in the variation over space of particle motion. More specifically, different areas within the cell induce different dynamics. We assume that we can describe these dynamics with a finite collection of motions. We must note that there also exists a vast literature analyzing the switch of motion over time. Again, the motions belong to a finite collection of transport states. Then, [START_REF] Monnier | Inferring transient particle transport dynamics in live cells[END_REF] propose a Bayesian model selection to infer the transient transport states along a trajectory over time. The switches between the different transient states are modeled with a Markov process. In their case, the different states correspond to Brownian motion with different diffusion coefficients and Brownian motion with drift with different values for both the drift and the diffusion coefficients.

Modeling and simulation of trapped particles

In this section, we propose a simulation framework based on general concepts in biophysics to mimic trapping areas. The artificial data sets will serve to evaluate the performance of detection and estimation methods described in Section 4. As said earlier in Section 1, we assume that the particles can switch between the trapped and the free states. When the particles are trapped, they undergo Brownian motion in confined regions. When the particles are free, they are driven by Brownian motion, only constrained by the boundary conditions on the limits of the domain D. We model the proportions of trapped and free particles in the whole domain and in the trapping regions through two systems of differential equations.

Modeling of particle dynamics

In the aforementioned situation, the dynamics of trapped and free particles are described by two connected processes:

(i) the dynamics of individual particles corresponding to normally reflected Brownian motion in D and in S when trapped;

(ii) the trapping process.

The trapping process induces a different motion for the trapped particles; inversely the particles can be trapped only if their motion drive them inside S. Formally, we assume m particles independently drawn from the uniform distribution over D:

X (i) 0 ∼ U(D), i = 1, . . . , m. Note that n i = n if
the trajectory length is the same for all the particles. Instead of considering the exact dynamics of a particle, we simply model its probability to be in any trapping region S p (equivalently to be in S = ∪ P p=1 S p ) by:

p S := P (X (i) t ∈ S) = |S| |D| , i = 1, . . . , m. (8)
Consequently, it is assumed that at every time t, the spatial point process (X

(1) t , . . . , X (m) t
) is a binomial point process over D of parameter m and which density function is the uniform density [START_REF] Baddeley | Spatial point processes and their applications[END_REF]. Furthermore, when the number m of particles is large enough, the evolution of the two populations can be modeled by a system of differential equations depending on parameters p S , k + and k -.

Let us denote b(t) = P (φ(X t ) = 1) and f (t) = P (φ(X t ) = 0) the proportions of trapped particles and free particles respectively. The functions t → b(t) and t → f (t)

are solutions of    db dt = γ 1 f (t) -γ 2 b(t), b(t) + f (t) = 1, (9)
where γ 1 and γ 2 are respectively the global binding rate and global unbinding rate, such that:

P (φ(X t+h ) = 1| φ(X t ) = 0) = γ 1 h + o(h), (10)
P (φ(X t+h ) = 0| φ(X t ) = 1) = γ 2 h + o(h),
with h > 0 and h → 0. Consequently, (φ(X t )) is a continuous-time homogeneous Markov chain with states {0, 1} with infinitesimal generator parameters k + and k - [START_REF] Brémaud | Markov chains: Gibbs fields, Monte Carlo simulation, and queues[END_REF][Chapter 8, Section 2.2]. Finally, we can show that (see Appendix A):

γ 1 = k + p S , (11)
γ 2 = k -.
The solution of the system ( 9) is as follows:

   b(t) = b(0) - k + p S k + p S + k -e -(k + p S +k -)t + k + p S k + p S + k -, b(t) + f (t) = 1. (12)
As a uniform distribution is assumed, we have the constraint that b(0) ≤ p S . In other words, only particles inside S can be trapped. The situation b(0) = p S matches with the situation where initially all the particles inside S are trapped. From Equation ( 12), the transitory regime is exponential and it converges to an equilibrium proportion:

lim t→∞ b(t) = k + p S k + p S + k -. ( 13
) If b(0) < k + p S /(k + p S + k -), function b(t) is decreasing toward its equilibrium point k + p S /(k + p S + k -). If b(0) > k + p S /(k + p S + k -), function b(t) is increasing toward its equilibrium point. If b(0) = k + p S /(k + p S + k -), function b(t)
is constant and equal to its equilibrium point. The characteristic time is defined as τ c = 1/(k + p S + k -). We can consider that the stationary regime is reached when t > 5τ c .

In the same way, we can model the proportions of trapped particles and free particles inside S. If t → b S (t) denotes the proportion of trapped particles inside S and t → f S (t) denotes the proportion of free particle inside S, we have:

   db S dt = k + f S (t) -k -b S (t), b S (t) + f S (t) = 1. ( 14
)
As before, the solution has an exponential form. The equilibrium proportion of trapped particles inside any trapping regions S p (same for all regions) is as follows

lim t→∞ b s (t) = k + k + + k -. (15)
Below, we describe the simulation of microdomains generated by running the

Fluosim software [START_REF] Lagardère | FLUOSIM: Simulator of single molecule dynamics for fluorescence live-cell and super-resolution imaging[END_REF] based on the set of aforementioned equations.

Table 1. Parameters of the first Fluosim simulation.

Type Parameters Value

Biophysics parameters

σ 2 1 µm 2 sec -1 k + 0.2 sec -1 k - 0.05 sec -1
Microscopy parameters ∆t 0.1 sec ∆x 0.025 µm First, we evaluate the model ( 12) describing the evolution of the proportion of trapped and free particles in the whole domain D. Second, we assess the model ( 14) describing the evolution of the proportion of trapped and free particles in the set of trapping areas S.

Simulations of trajectories and microdomains

In our first experiment, we designed two trapping regions S 1 and S 2 in a a square domain D ⊂ R 2 with radius 5 µm (origin (0, 0) at the bottom left corner of D):

(i) S 1 is a circle with radius r 1 = 0.65 µm and center (2.5, 2.5).

(ii) S 2 is a circle with radius r 2 = 0.39 µm and center (4, 4).

This situation typically corresponds to the analysis of a region of interest as illustrated in Fig. 1(b) (yellow square). The parameters of simulation are given in Table 1. Here, 2 500 particles are uniformly distributed over D. Fluosim allows us to identify which particle is trapped at time t. In that case φ(X (i) t ) is known, meaning that if the particle i is trapped or not trapped at time t. In Fig. 3, we plot the positions of the particles labelled as free or trapped at time t = 10 sec (transitory regime) and t = 100 sec (stationary regime). 12); green: asymptote of the model ( 12). We computed the red curve over m = 2 500 trajectories simulated with Fluosim (see Table 1).

Evaluation of proportions of trapped and free particles in the whole domain

We compare the proportion of trapped particles t → b(t) corresponding to Equation ( 12) to the true proportion of trapped particles. The curves are shown in Fig.

(4)(a). The fit of the transitory regime (t < 5τ c = 77.5 sec) is rather good even if the stationary regime is not exactly the same as the one predicted by the model. The predicted equilibrium is 22.53% while the mean true proportion of trapped particles computed on the last 100 steps of time is 21.35%. The relative error of the model is 5.24%, mainly due to the fact that the model oversimplifies the underlying trapping process; the dynamics of the particles are only controlled by the parameter p S , and we assume that the particles are always in the stationary regime of a normally reflected Brownian motion on D.

This assumption holds at the beginning, as the particles are initially drawn from the uniform distribution on D. Nevertheless, as t increases, the trapping process makes this assumption invalid. This explains the relative good fit of the transitory phase and the relative lack of fit of the asymptotic phase. However, we can consider that it is a rather satisfying modeling, regarding the flexibility and the number of parameters.

Evaluation of proportions of trapped and free particles in trapping areas

We compare the proportion of trapped particles t → b S (t) given by Equation ( 14) to the true proportion of trapped particles inside S. Interestingly, we observe the opposite behavior compared to the case of b(t) (proportion of trapped particles in the whole domain D). There is a clear lack of fit during the transitory regime but a perfect fit during the stationary regime (see Fig.

(4)(b). We can see that the transitory regime has the same duration for b(t) and b S (t). During this phase, the number of particles inside S increases due to the process of Brownian particles entering in S and getting trapped. However, it is explicitly assumed in Equation ( 14) that the number of particles is constant (and large). This explains the lack of fit during the transitory phase which matches with a period during which the number of particles increases in S. On the contrary, once b(t) has reached the stationary regime, the number of particles inside S is approximately constant (even if it is not the same particles that remain in S from one time to another). In the latter case, the model ( 12) is relevant and the stationary regime is well predicted by the model.

Clustering of trajectories for microdomain detection

In this section, we describe the computational approach used to automatically detecting microdomains where particles undergo subdiffusion. The proposed approach is not based on the biophysical modeling and the Fluosim simulator presented earlier. Fluosim serves here to quantitatively demonstrate the performance of our approach on synthetic sequences (ground truth).

Our four-step detection procedure is based on the combination of the clustering algorithm Dbscan [START_REF] Ester | A density-based algorithm for discovering clusters in large spatial databases with noise[END_REF] and the labeling of trajectories as explained in [START_REF] Briane | Statistical analysis of particle trajectories in living cells[END_REF]. The procedure can be used for a large range of application problems in live cell imaging, and is summarized as follows (see Fig. 5):

Step #1: Classification of trajectories: we run our test procedure [START_REF] Briane | Statistical analysis of particle trajectories in living cells[END_REF] on the collection of trajectories X m (Fig. 5(b)) to label the tracks into three diffusion groups: free diffusion (Brownian motion), superdiffusion, and subdiffusion. We denote Ω 1 (X m ) the set of trajectory indexes corresponding to the acceptance of hypothesis H 1 (subdiffusion).

Step #2: Representative point selection: we choose a unique point xi to represent each trajectory X (i)

n i (Fig. 5(c)). We remove points associated to superdiffuse trajectories and Brownian trajectories (Fig. 5)(d).

Step #3: Point clustering: We partition the set

Ω = {x i | i ∈ Ω 1 (X m )} into clusters (Fig. 5(e)).
Step #4: Detection of trapping areas: We use the clusters to delineate confinement areas and microdomains (Fig. 5(f)).

In the sequel, we successively present Steps #2, #3 and #4 of the method. As already explained, we will use the simulation framework presented in Section 3 to illustrate the different steps of our method.

Classification of trajectories

Step #1 is straightforward as it corresponds to the test procedure described in [START_REF] Briane | Statistical analysis of particle trajectories in living cells[END_REF]. This method is based on a non-parametric three-decision test to distinguish the three types of diffusion (superdiffusion, subdiffusion, free diffusion). The null hypothesis is that the observed trajectory made of n points is generated from a Brownian motion and the two other alternatives are subdiffusion and superdiffusion. The test statistic T n is the standardized largest distance covered by the particle from its starting point. It can be interpreted as follows: (i) if the value of T n is low, it means that the process stayed close to its initial position and can be trapped in a microdomain (subdiffusion); (ii) if the value of T n is high, the particle went far from its initial position and the particle is actively transported (e.g. via a molecular motor) (superdiffusion).

The test is appropriately normalized, asymptotically converges under the null hypothesis and under parametric alternatives which are currently considered in the biophysics literature. The test procedure, applied to subdiffusive trajectories, is easily controlled by a p-value.

Our procedure is applied to detect areas where subdiffusion occurs, but it can also be used to detect areas in which superdiffusion or Brownian motion are the main dynamics. In that case, we just need to replace the set Ω 1 (X m ) by respectively Ω 0 (X m ) (set of trajectory indexes corresponding to Brownian motion) or Ω 2 (X m ) (set of trajectory indexes corresponding to superdiffusion).

Moreover, if we observe long trajectories, dynamics can vary, switching from one diffusion model to another. In this case, labelling a trajectory with a single diffusion class would be misleading. To overcome this problem, we can use the algorithm described in [Briane et al., 2019b] (also based on the classification proposed in [START_REF] Briane | Statistical analysis of particle trajectories in living cells[END_REF]) to detect the change points where the diffusion switches occur and segment the trajectories accordingly. Then we can use these sub-trajectories as inputs of the proposed algorithm.

Finally, it is possible that, due to the spatial extensions of the trapping regions, a trapped particle never hits the edges of the region during the period of observation. In this case, these particles have no chance to be detected as subdiffusive by the procedure described in [START_REF] Briane | Statistical analysis of particle trajectories in living cells[END_REF]. If we assume that the trapped particles have a significantly lower diffusion coefficient than the untrapped particles -a realistic biological assumption -we can use this supplementary assumption to detect the trapped particles even in this case. Consequently, we proceed as follows. We keep all the particles classified as subdiffusive as before. For all the particles classified as Brownian we can compute an estimator of the diffusion coefficient as follows:

σ2 n = 1/(2n∆) n j=1 X t j -X t j-1 2 2 .
In the next step, we classify the trajectories first detected as Brownian by using the method [START_REF] Briane | Statistical analysis of particle trajectories in living cells[END_REF] based on the estimated diffusion coefficient, and we re-label as subdiffusive all the trajectories detected as Brownian by the method [START_REF] Briane | Statistical analysis of particle trajectories in living cells[END_REF] for which the estimated diffusion coefficient is low enough. A threshold must be set by either a statistical procedure -for example based on the value of the diffusion coefficients of the trajectories detected as subdiffusive in the first place -or by an expert. These trajectories would correspond to the subdiffusive trajectories that never hit the edges of the subregions. We can then run the next steps of the procedure on this new set of subdiffusive trajectories.

Representative trajectory points and spatial distribution

In Step #2, a single point is required to represent a trajectory made of n points. In what follows, we consider mapping functions f : R 2×n → R 2 of the following form:

• Mean point : f (x 1 , . . . , x n ) = 1 n n i=1
x i ,

• j-th point located on the trajectory: f (x 1 , . . . , x n ) = x(j), with j ∈ {1; . . . , n}.

As already mentioned, the trapped particles are modeled by confined Brownian motion normally reflected at the boundaries of the trapping regions. In that case, we can see that different representative trajectory points yield very different spatial patterns. As illustrated in Fig. 6, we represent the trajectory X n = (X t 0 , . . . , X tn ) by the mid-point X n/2 (related function f (x 1 , . . . , x n ) = x( n/2 )). It turns out that the spatial distribution of the representative points

X (i) n/2 of the trapped trajectories X (i)
n is uniform with this choice, over the trapped region S 1 and S 2 (Fig. 6(b)). This is due to Proposition 1, as the trapped particles undergo a normally reflected confined Brownian motion in S 1 and S 2 .

Proposition 1 [START_REF] Pinsky | Asymptotics of the principal eigenvalue and expected hitting time for positive recurrent elliptic operators in a domain with a small puncture[END_REF]. Let (X t ) be a normally reflected Brownian motion on a finite volume domain D. Assume that the process (X t ) has a stationary uniform distribution over D: U(D). If the process (X t ) has reached its stationary distribution, then for any t > 0 and B ⊂ D it follows that

P (X t ∈ B) = |B| |D| , ( 16
)
where |B| denotes the area of domain |B|.

Actually, the average position of the trapped trajectories are concentrated in the center of the trapped regions (Fig. 6(a)). Therefore, in the purpose of accurately estimating the boundaries of S 1 and S 2 , it seems more appropriate to represent the trajectory i with X (i) n/2 rather than by the average position, not necessary located on the trajectory.

A clustering algorithm: Dbscan

In Step #3, we apply the popular clustering algorithm Dbscan [START_REF] Ester | A density-based algorithm for discovering clusters in large spatial databases with noise[END_REF] to robustly estimate point clusters from noisy points due to presence of false detections.

Dbscan adds a point into a cluster if it is surrounded by at least κ neighbouring points in a disk of radius -this type of point is called a core point -or if it is within a disk of radius from a core point -this type of point is called a boarder point -. On the contrary, the point does not belong to any cluster and is considered as noise. Other more sophisticated clustering algorithms could be applied. For instance, the so-called a contrario method allows a finer analysis of the clusters than Dbscan through the concept of meaningful clusters [START_REF] Desolneux | Maximal meaningful events and applications to image analysis[END_REF][START_REF] Briane | A sequential algorithm to detect diffusion switching along intracellular particle trajectories[END_REF] organized in a hierarchical structure. We focus on the Dbscan algorithm since it is very popular spatial statistics. It is mainly controlled by two parameters ( , κ) as explained below.

Selection of the Dbscan parameters

The set of parameters ( , κ) should be adapted to each cluster, as each cluster does not have the same density (or concentration) of points. However this information is available once the clusters are estimated. In practice, a unique set of parameters ( , κ) is used for all clusters. A good choice is to set the parameters ( , κ) in order to detect the least dense cluster. Indeed, the clusters with more points are more easily detected than the least dense clusters. [START_REF] Ester | A density-based algorithm for discovering clusters in large spatial databases with noise[END_REF] argue that, in two-dimensional problems, we can set κ = 4 as a rule of thumb. Once the parameter κ is set, we can choose by considering two different strategies:

(i) a data driven method based on the observed distribution of κ-nearest neighbors;

(ii) a parametric method to approximate the distribution of noisy points.

In the following, we present two data driven methods and two parametric methods for selecting . Then, we compare the outcomes of the Dbscan algorithm with different methods for selecting , applied to Fluosim simulations.

Data driven methods

The data driven methods enable to determine from the sample d 1 (κ), . . . , d m (κ) where d i (κ) denotes the distance of the κ-nearest neighbors of point x i ∈ Ω. Let us denote d (1) (κ), . . . , d (m) (κ) the increasing-ordered sample. Note that if = d (i) (κ), the i points corresponding to the values d (1) (κ), . . . , d (i) (κ) will be core points while the other points will be either boarder points or noisy points. The choice < d 1 (κ) corresponds to the case when all the points are noisy points (no core point), whereas the choice ≥ d m (κ) corresponds to the case where all the points are core points. In the latter case, all the points belong to one single cluster. If the number of noisy points ℵ is known, a natural choice proposed by [START_REF] Ester | A density-based algorithm for discovering clusters in large spatial databases with noise[END_REF] is = d (m-ℵ+1) (κ). A alternative graphical approach is also described in [START_REF] Ester | A density-based algorithm for discovering clusters in large spatial databases with noise[END_REF] to determine . The authors plot the sequence d (m) (κ), . . . , d (1) (κ) and is determined as the first d (i) in the first valley of the sequence.

Parametric methods

In [START_REF] Daszykowski | Looking for natural patterns in data: Part 1. density-based approach[END_REF], it is implicitly assumed that the noisy points are uniformly distributed. The parameter is defined as the α quantile of the distribution (e.g. α = 0.05) of the κ-nearest neighbors of |Ω| points drawn from a uniform distribution over D. If the observed noisy points are really drawn from a uniform distribution, none of them will be chosen as a core point with approximately 1 -α probability.

As the number ℵ of noisy points and the clusters C i are unknown, it is not possible to propose a method which exactly controls the probability of detecting a noisy point as Table 2. Values of the parameter of Dbscan obtained with different methods. The parameter is estimated assuming κ = 4 as recommended in [START_REF] Ester | A density-based algorithm for discovering clusters in large spatial databases with noise[END_REF]. The parameter is expressed in pixel units. The image size is 200 × 200 pixels.

Type

Method Data driven Graphical method 5.64 Parametric [START_REF] Daszykowski | Looking for natural patterns in data: Part 1. density-based approach[END_REF] 5.94 [START_REF] Friedman | An algorithm for finding nearest neighbors[END_REF] 5.86 a core point. In [START_REF] Daszykowski | Looking for natural patterns in data: Part 1. density-based approach[END_REF], is estimated with Monte-Carlo simulations, which is time consuming if |Ω| is large.

It is worth noting that [START_REF] Friedman | An algorithm for finding nearest neighbors[END_REF] also studied the distribution of the κ-nearest neighbors of m points uniformly distributed over a finite domain D. If |Ω| is large enough to neglect boundary effects, the ratio of the volume of a two-dimensional sphere centered at a point containing κ neighbors and the volume of the whole space D is governed by a Beta distribution f with parameters (κ, |Ω| -κ):

f (x) = m! (κ -1)!(m -κ)! x κ-1 x m-κ , 0 ≤ x ≤ 1. (17)
Finally, we can define as explained in [START_REF] Daszykowski | Looking for natural patterns in data: Part 1. density-based approach[END_REF]. If F -1 (α) denotes the quantile of order α of the distribution f , we get:

= |D|F -1 (α) π 1/d , ( 18
)
where |D| is the volume (or area in the two-dimensional case) of D.

Comparison of strategies on Fluosim simulations

We apply the Dbscan algorithm on the set Ω obtained with the test procedure [START_REF] Briane | Statistical analysis of particle trajectories in living cells[END_REF] (see Fig.

5(b))

. The parameter is estimated with the different methods. All the methods provided similar estimations of in our experiments (e.g. see Table 2). Consequently, we got the same set of clusters when we run Dbscan as shown in Fig. 7 by setting κ = 4. We found two clusters as expected, corresponding to the two trapping regions S 1 and S 2 . We got similar results in additional experiments, suggesting that selecting one method or another one would not affect the detection of trapping areas.

Detection of trapping areas

In

Step #4, we simply use the convex hull to define the boarders of the clusters. By definition of the convex hull, all the points of the cluster will be either inside the convex hull (or convex polygon) or be some vertexes of the polygon. We note that several extensions of the convex hull have been developed, including α-shapes [START_REF] Edelsbrunner | On the shape of a set of points in the plane[END_REF]. In that case, an extra parameter α > 0 is used to control the smoothness of the estimated shape.

Experimental results

Evaluation on Fluosim simulated sequences

In this section, we use the Fluosim simulation scheme to illustrate the potential of our approach to detect trapping areas. We consider temporal series depicting 4 500 trajectories made of n = 350 points. The trajectories are analyzed once the equilibrium regime is reached. We recall that the trajectories can switch between two diffusion modes:

(i) confined Brownian in small microdomains S p with normal reflection on the boundaries ∂S p ;

(ii) Brownian motion in the whole domain D with normal reflection on the boundaries ∂D.

We designed a cell in with a crossbow shape and seventeen trapping regions depicted in Figure 8:

• one circular nucleus surrounded by one crescent-shaped zone,

• five circular zones with decreasing radius,

• seven tiny circular zones,

• three tiny circular zones really close to each other.

The parameters are given in Table 3. In order to assess the efficiency of the method, we compute two criteria for each true (ground truth) trapping region. First, we compute the number of detected regions overlapping with the real trapping regions. Secondly, we compute the proportion of overlap between the closest detected region to the trapping region (useful if multiple detected regions overlap with the true region). We define the proportion of overlap as the area of intersection of the true region and the detected region over the area of the true region. Then, ideally the two criteria should be equal to one, meaning that one detected region corresponds to one real trapping region and that the overlap is perfect. Also, for each detected cluster we compute the proportion of subdiffusive trajectories which should be compared to Equation ( 15). In this case, the theoretical proportion is k + /(k + + k -) = 0.8333. Table 5. Proportion of subdiffusive trajectories inside the detected regions from the crossbow-shaped cell simulation.

Cluster index Proportion of subdiffusive trajectories 1 3.3 % 2 5.6 % 3 18.5 % 4 8.3 %

From Table 4, we can see that the algorithm detects satisfyingly the regions #1 to #5 while it does not detect at all the regions #6 to #17. The best detected region is region #1 with 90% of overlap. The algorithm does not detect the small regions as there is either no detected subdiffusive trajectories or less than κ = 4 detected subdiffusive trajectories (see Figure 9 (a)). In the latter case, Dbscan can not build a cluster as the minimum number of points in a cluster must be κ = 4. We also note that the proportion of subdiffusive trajectories in the five detected regions is much lower as reported in Table 5, than the expected value k + /(k + + k -) = 0.8333. Then, in the setting of this simulation, the algorithm under-detects subdiffusive trajectories. The parameters that influence the most the outcome of the classification step (Step #1) are the trajectory length, the microdomain size, and the diffusion coefficient of trapped particles. An optimal setting is a long trajectory with a high diffusion coefficient in a small domain. In this case, the particles bounces a lot against the edges of the trapping regions and its subdiffusive behavior is obvious. A more detailed simulation scheme would give the interval of parameters for which the classification step performs the best.

Another issue is the proximity between two trapping regions. We can see in Figure 9 (b) that a detected region overlaps both with region #1 (large disk) and region #2 (crescent shape). We note that the minimum distance to resolve two close regions will depend on the concentration of subdiffusive trajectories, that is on the binding parameters k + and k -. In fact, the outcome of Dbscan depends on the density of points inside the clusters. In Appendix B, we provide another simulation controlled by the parameters given in Table 6. The results are shown in Figure 12 and the quantitative results are reported in Tables 7 and8. In conclusion, we did not choose the optimal simulation parameters here to assess our method. Still, we are able to detect some trapping regions even if the algorithm under-detects subdiffusive trajectories and merges close trapping regions. A more systematic simulation study should be performed to determine the range of parameters for which the method gives satisfying results.

Experiments on real TIRF images: dynamics of Rab11 involved in exocytosis

In this section, we demonstrate the potential of our approach on Rab11-mCherry proteins involved in exocytosis mechanisms and observed in 2D TIRF (Total Internal Reflection Fluorescence) microscopy (Nikon TE2000 inverted microscope equipped with a 100× TIRF objective: numerical aperture = 1.49, exposure time ∆t = 100 msec, dx = dy = 160 nm). Figure 10 shows three 256 × 353 images at different time points, depicting fluorescent vesicles close to the plasma membrane. The image sequences are typically composed of 600 images, and show the late steps of exocytosis driven by Rab11 at 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].

Given several hundreds of Rab11 trajectories estimated with the Icy software [START_REF] Chenouard | Multiple hypothesis tracking for cluttered biological image sequences[END_REF][START_REF] Chenouard | Objective comparison of particle tracking methods[END_REF], we first applied our test procedure [START_REF] Briane | Statistical analysis of particle trajectories in living cells[END_REF] to classify the trajectories. In Fig. 11(a), we displayed all trajectory points onto a 2D three-color map: Brownian motion (blue), subdiffusion (green), and superdiffusion (red). In the particular case shown in Fig. 10, the M10 cells were treated with Latrunculin A, which inhibits actin polymerization and then reduces active transport of vesicles [START_REF] Pécot | A quantitative approach for analyzing the spatio-temporal distribution of 3d intracellular events in fluorescence microscopy[END_REF]. This effect is confirmed on Fig. 11(a) since no track was labeled a superdiffusive trajectory.

Second, we applied our computational approach on subdiffusion trajectories (green pixels, left) represented by points in Fig. 11(a) (right). In the sequel, we assume that the tethering-docking process under study is associated to subdiffusion because of multiple interactions with actin filaments and interactions with the plasma membrane. All points reported in Fig. 11(a) (right) correspond to vesicles that fuse at the surface of the plasma membrane. First, we observe an approximately uniform spatial distribution of points (black dots). Nevertheless, we detected a few additional clusters corresponding to higher densities of fusion events in a local spatial neighborhood. These clusters can be interpreted as local domains related to more intensive interactions with the plasma membrane and actin filaments. Several clusters are actually localized in regions with a high concentration of spots (e.g. see top right of the cell). Trajectory analysis is then helpful to focus the attention on specific areas which cannot be easily detected by examining fluctuations of fluorescence intensities shown in Fig. 10. Finally, we illustrate in Fig. 11(b) the potential of the [START_REF] Hoze | Heterogeneity of ampa receptor trafficking and molecular interactions revealed by superresolution analysis of live cell imaging[END_REF] method on the same sequence. In particular, we reported the spatially varying diffusion coefficients on a 2D map since the drift vectors are very noisy and not interpretable. In our study, the diffusion map is actually not really helpful to draw conclusions if the drift values are noisy. The number of estimated trajectories (about one thousand) is too small to safely apply the methods described in [Hoze et al., 2012, Salomon and[START_REF] Salomon | Precise mapping of intracellular diffusion and drift from SPT data analysis[END_REF].

Discussion

In this paper, we presented a method to detect the trapping regions where the particles are confined and thereby undergo subdiffusion. We combined the clustering algorithm Dbscan [START_REF] Ester | A density-based algorithm for discovering clusters in large spatial databases with noise[END_REF] and a test procedure [START_REF] Briane | Statistical analysis of particle trajectories in living cells[END_REF] in order to identify the areas with a high concentration of subdiffusive particles, assumed to correspond to trapping regions. With this approach, we were also able to estimate the proportion of particles inside the trapping regions which were indeed confined in the domain. Our approach is more flexible than the Eulerian approach of [START_REF] Hoze | Heterogeneity of ampa receptor trafficking and molecular interactions revealed by superresolution analysis of live cell imaging[END_REF], which is not able to handle mixtures of motions, or to capture the motion of particles not confined in the trapping regions but still going through these regions.

We illustrated the potential of our procedure on simulations generated with the Fluosim software, and on real TIRF images. However, the simulation results provided here do not intend to assess systematically the procedure for a large range of parameters values (k + , k -, σ S , σ T , the shapes of the regions). For future works, it would typically involve:

• determining jointly the region sizes and trajectory lengths for which the classification procedure can detect subdiffusive trajectories (Step #1) and further defining a criterion to select the region size for different shapes;

• determining the minimum distance between two regions for the clustering algorithm to be able to detect two distinct regions and not merge them. This distance is expected to depend on the concentration of subdiffusive trajectories in the trapping regions.

In this paper, we focused on the modeling and estimation framework and demonstrated its potential with preliminary experiments. Then, we emphasize the fact that it is possible to use different procedures at each step of the proposed pipe-line to adapt to different scenarios. For instance, in the case of trajectories switching over time, we can use sub-trajectories estimated with the method described in [Briane et al., 2019b], instead of the classification procedure [START_REF] Briane | Statistical analysis of particle trajectories in living cells[END_REF]

used in

Step #1 of our algorithm. The same way, another choice of representative point and other strategies of clustering or shape estimation can also be used respectively in Steps #2 #3 and #4.

We proposed a default choice of procedures for each step of the pipeline but the user is free to select the appropriate procedures depending on the application. Table 7. Evaluation of the procedure on the second simulated crossbow-shaped cell.

True regions Region index (see Fig.

Figure 1 .

 1 Figure 1. Illustrations of configurations with trapping areas. (a) sketch with three trapping areas (S = S 1 ∪ S 2 ∪ S 3 ); (b) example of a simulated image at a given time t, depicting a crossbow-shaped cell and green particles in small trapping (blue) areas and in the cytosol S (orange).

Figure 2 .

 2 Figure2. Drift field computed with the method of[START_REF] Hoze | Heterogeneity of ampa receptor trafficking and molecular interactions revealed by superresolution analysis of live cell imaging[END_REF]. The Fluosim simulation comprises 30 frames of 200 × 200 pixels as the domain D is a square of 5µm and the spatial resolution is ∆x = 0.025 µm (see Table1). We compute μ(x) on each point of this 200 × 200 lattice from trajectories (a) (r = 5 pixels). We can clearly see the two microdomains S 1 and S 2 (b): the vector field converges to the centers of these domains delineated by two red circles (b). We can see that the vector field is orthogonal to the boundaries of the square D. It is due the normal boundary condition on ∂D. Elsewhere, the μ(x) is almost null, reflecting Brownian motion.

  Figure 3. Positions of the particles simulated with Fluosim at time t = 10 sec and t = 100 sec (blue: free particles; red: trapped particles). The red circles delineate the regions S 1 (large disk) and S 2 (small disk).

Figure 4 .

 4 Figure 4. Evolution of the proportions of trapped particles over time. (a) proportions of trapped particles in D; (b) proportions of trapped particles in the trapping region S. red: true proportion of trapped particles computed from the data; blue: proportions computed with model (12); green: asymptote of the model (12). We computed the red curve over m = 2 500 trajectories simulated with Fluosim (see Table1).

Figure 5 .

 5 Figure 5. Overview of our computational approach for detecting trapping areas. (a) input trajectories; (b) classification of trajectories into three diffusion groups (superdiffusion (red), subdiffusion (green), free diffusion (blue)); (c) representative points on trajectories; (d) display of subdiffusive representative points (green); (e) estimation of clusters associated to trapping areas; (f) detection of trapping areas.

Figure 6 .

 6 Figure 6. Spatial distribution of the trapped particles in the Fluosim simulation; (a) the trajectories are represented by their mean point; (b) the trajectories are represented by X n/2 . A trajectory is considered to be trapped if it is trapped during the whole period of observation. The red circles represent the boundaries of the trapping regions S 1 (large circle) and S 2 (small circle).

Figure 7 .

 7 Figure7. Clusters detected by Dbscan (κ = 4) on the set of particles detected as subdiffusive with the test procedure[START_REF] Briane | Statistical analysis of particle trajectories in living cells[END_REF] (single test with a significance level of 5%). We choose = 5.86 derived with the Beta distribution assumption. The noisy points are in red, and points belonging to the same cluster are labelled with the same color. The red circles represent the boundaries of the trapping regions S 1 (large circle) and S 2 (small circle). The black square delineates the borders of the whole domain D.

Figure 8 .

 8 Figure 8. Simulated crossbow-shaped cell with seventeen trapping regions of different sizes and shapes. The white line at the bottom accounts for 16 µm.

Figure 9 .

 9 Figure 9. Results of our procedure on the first crossbow-shaped cell simulation; (a) representative points of the trajectories detected as subdiffusive by the test procedure [Briane et al., 2018] with a significance level of 5%; (b) trapping regions (yellow background, red boundaries) estimated with the convex hull from the clusters detected by Dbscan. The true trapping regions are delineated with blue curves.

Figure 10 .

 10 Figure 10. Rab11-mCherry proteins within a M10 cell treated with Latrunculin and observed at several time points. The images are extracted from a 2D TIRFM sequence composed of 600 images (exposure time ∆t = 100 msec, dx = dy = 160 nm).

Figure 11 .

 11 Figure 11. Analysis and clustering of trajectories. (a) Trajectory classification into three diffusion groups [Briane et al., 2018] (left) and selection/clustering of subdiffusion trajectories (right). The point clusters are labeled with colors while individual subdifffusion trajectory are represented with black dots. (b) Adaptation of the Hoze et al. [2012] method to estimate diffusion coefficients on a sparse set of unclassified trajectories [Salomon and Kervrann, 2020].

Figure 12 .

 12 Figure 12. Results of our procedure on the second simulated crossbow-shaped; (a) representative points of the trajectories detected as subdiffusive by the test procedure of [Briane et al., 2018] at level 5%; (b) trapping regions (yellow background, red boundaries) estimated with the convex hull from the clusters detected by Dbscan. The true trapping regions are depicted in blue.

  Proportion of subdiffusive trajectories inside the detected regions from the crossbow-shaped cell simulation (2) to be compared (k + /(k + + k -) = 0.9091). ) 42.0 19.0 5.0 15.0 8.5 11.0 11.0 8.0

Table 3 .

 3 Fluosim parameters to simulate the first crossbow-shaped cell.

	Type	Parameter	Value
		σ 2 S	0.2 µm 2 .s -1
	Biophysics parameters	σ 2 T k +	0.02 µm 2 .sec -1 10 sec -1
		k -	2 sec -1
	Microscopy parameters	∆t ∆x	0.02 sec 0.160 µm

Table 4 .

 4 Evaluation of the procedure on the crossbow-shaped cell simulation.

	True regions

Table 6 .

 6 Fluosim parameters to simulate the second crossbow-shaped cell.

	Type	Parameter	Value
		σ 2 S	0.2 µm 2 .s -1
	Biophysics parameters	σ 2 T k +	0.1 µm 2 .s -1 10 sec -1
		k -	1 sec -1
	Microscopy parameters	∆t ∆x	0.02 sec 0.160 µm

Acknowledgement

We thank M. Lagardère (Interdisciplinary Institute for Neuroscience, Bordeaux, France) for assistance with the Fluosim software and for helpful discussion.

Funding

Funding was provided by Inria Rennes and CREST-Ensai-Université Bretagne Loire. This work was also supported by the French National Research Agency (France-BioImaging infrastructure-ANR-10-INBS-04, DALLISH-ANR-16-CE23-0005).

Appendix A: Derivation of the global binding and unbinding rates

We demonstrate that the equalities (12) hold. Let us consider the equality involving γ 1 :

The particle can not be trapped outside S. Consequently, the second probability of the sum is zero. As the probability to be in region S at t + h is independent from the event {φ(X t ) = 0} to be a free at t, we finally get

Normally reflected Brownian motion in a domain D with finite volume has a stationnary distribution. This distribution is the uniform distribution over D [START_REF] Pinsky | Asymptotics of the principal eigenvalue and expected hitting time for positive recurrent elliptic operators in a domain with a small puncture[END_REF]. Then, under the uniform distribution assumption, we get:

It follows that γ 1 = ' k + |S|/|D|. Now we prove the second equality of (12) involving γ 2 .

As {φ(X t ) = 1} ⊂ {X t+h ∈ S}, we get the result:

Appendix B: Supplementary evaluation on a second crossbow-shaped cell

In this section, we provide another simulation controlled by the parameters given in Table 6.

The results are shown in Figure 12 and the quantitative results are reported in Tables 7 and8 (see Section 6.1 for comments).