Yunjiao Lu

Pierre Hodara

Charles Kervrann

Alain Trubuil

Probabilistic reconstruction of truncated particle trajectories on a closed surface

Investigation of dynamic processes in cell biology very often relies on the observation in two dimensions of 3D biological processes. Consequently, the data are partial and statistical methods and models are required to recover the parameters describing the dynamical processes. In the case of molecules moving over the 3D surface, such as proteins on walls of bacteria cell, a large portion of the 3D surface is not observed in 2D-time microscopy. It follows that biomolecules may disappear for a period of time in a region of interest, and then reappear later. Assuming Brownian motion with drift, we address the mathematical problem of the reconstruction of biomolecules trajectories on a cylindrical surface. A subregion of the cylinder is typically recorded during the observation period, and biomolecules may appear or disappear in any place of the 3D surface. The performance of the method is demonstrated on simulated particle trajectories that mimic MreB protein dynamics observed in 2D time-lapse fluorescence microscopy in rod-shaped bacteria.

INTRODUCTION

In 2D and 3D live-cell imaging, spatiotemporal events and biomolecule dynamics are frequently observed with an incomplete field of view. Very often these observations are related to regions of observation (ROO) inside a tissue, a cell, or in the neighborhood of membranes. Nevertheless, it is quite unusual to analyze 3D dynamics of biomolecules or events occurring on a closed surface and observed on a 2D plane. Our work is motivated by the study of dynamics of MreB proteins, moving close to the inner membrane during cell wall construction in rod-shaped bacteria ( [START_REF] Billaudeau | Contrasting mechanisms of growth in two model rod-shaped bacteria[END_REF], [START_REF] Van Teeffelen | Recent advances in understanding how rod-like bacteria stably maintain their cell shapes[END_REF]). Its dynamics can only be observed in a small region and are recorded as 2D time-lapse movies (Fig. 1a). As for 3D image acquisition, even it can solve the problem of partial observation, is not always appropriate, especially if the objective is to capture fast and temporally short events as described in [START_REF] Billaudeau | Contrasting mechanisms of growth in two model rod-shaped bacteria[END_REF]. The frame rate adapted to the scale of dynamics may be too high when compared to the period of time to acquire temporal series of 3D volume ( [START_REF] Boulanger | Fast high-resolution 3d total internal reflection fluorescence microscopy by incidence angle scanning and azimuthal averaging[END_REF] and [START_REF] Cornilleau | Methods for studying membrane-associated bacterial cytoskeleton proteins in vivo by tirf microscopy[END_REF]).

To our knowledge, identifying re-entrance events of the same entities inside the ROO is not addressed in the literature. In experimental data, when the unobserved region represents a significant part of the entire surface, a complete description of the dynamics on these closed surfaces becomes of paramount importance for deciphering the mechanisms of some processes. In our study of the regulation of the dynamics of MreB protein, as inputs, we consider a set of trajectories estimated by tracking algorithms (e.g. [START_REF] Jaqaman | Robust single-particle tracking in live-cell time-lapse sequences[END_REF], [START_REF] Chenouard | Multiple hypothesis tracking for cluttered biological image sequences[END_REF]). These tracking algorithms are very sophisticated and allow us to handle large sets of particles, different stochastic dynamical models [START_REF] Blom | The interacting multiple model algorithm for systems with markovian switching coefficients[END_REF], [START_REF] Bressloff | Stochastic models of intracellular transport[END_REF], and observation models [START_REF] Genovesio | Multiple particle tracking in 3-d+ t microscopy: method and application to the tracking of endocytosed quantum dots[END_REF], [START_REF] Li | Target perceivability and its applications[END_REF]. They take into account birth/death events, and/or split/merge events. Particles may be unobserved or undetected for short periods of time, especially in 2D+time microscopy. However, none computational or statistical method manages the situation corresponding to a large hidden region inside the region of interest. Also, the identification of particles leaving the ROO through one border of the domain and re-entering from a far border is not addressed. Our objective is then to provide a generic approach to tackle the problem of the reconstruction of particle trajectories observed on a small part of a closed surface as illustrated in Fig. 1b.

In this paper, we focus on the design and evaluation of a self-contained mathematical framework to tackle the reconstruction of particle trajectories on cylindrical surfaces, given the tracklets observed in a small window sampled on the surface. In our study, the particles are assumed to obey a stochastic Brownian motion with drift and may appear or disappear during the observation period. Split or merge events are not considered in the modeling framework. The trajectory reconstruction problem is defined as the maximization of the likelihood function given tracklets inside the ROO. The optimization problem to be solved is formulated as an integer linear programming problem. The final algorithm is a data-driven algorithm with no hidden parameter to be set by the user. We demonstrate the performance and robustness of our computational method on simulation data, by varying the ratio of observed to unobserved region, the drift and variance of particles, as well as the rates of birth and death of particles.

The remainder of the paper is organized as follows. In Section 2, we present the problem formally and introduce notations. In Section 3, we describe the probabilistic framework, including Poisson processes used to describe birth and death events, and Brownian motion with drift to represent particle motion. We also describe the computational procedure aiming at connecting tracklets belonging to the same trajectory, and then recovering the dynamics of particles on the whole closed surface. Note that we suppose that the curvature of the cylinder is known and so that the movements are represented on a 2D unwrapped surface. In Section 4, the performance of our algorithm is evaluated on simulated data. Finally, we conclude and propose some future work. A summary of notations useful for the evaluation of the likelihood is given in Supplementary Materials [START_REF] Blom | The interacting multiple model algorithm for systems with markovian switching coefficients[END_REF].

PROBLEM STATEMENT AND NOTATIONS

We consider a probabilistic model to represent particles that are born, move and die on a cylindric membrane. Formally, let us denote H and L the height and perimeter of the cylinder respectively (see Fig. 1). We associate 2D coordinates (x, y) ∈ [-L, 0] × [0, H] to each point of the underlying cylindric manifold. The particles are "born" with a constant rate λ and appear uniformly at random on the membrane surface. We consider a Poisson process with intensity λ to statistically represent the birth events. Each particle is assumed to have the same constant rate of death τ d such that life duration T d of a particle follows an exponential law of parameter τ d . During its lifetime, a particle k born at time t 0 and located at Z k 0 = (X k 0 , Y k 0 ), moves according to Brownian motion with drift. On the set ] -L, 0[×[0, H], the position of the particle at time t ≥ t 0 prior to its death time is given by

Z k t = Z k 0 + v(t -t 0 ) + ΣB k t-t 0 (1)
where

Z k t = (X k t , Y k t ), v = (v x , v y ), Σ = σ x 0 0 σ y , B k t is a two-dimensional Wiener process.
In order to model the topology of the cylinder as illustrated in Fig. 1, we impose deterministic jumps when the process reaches one of the two borders {-L} × [0, H] or {0} × [0, H]. For any y ∈ [0, H], the process reaching position (-L, y) jumps to position (0, y) and vice versa.

In y direction the initial position of a particle lies between [0, H]. When a particle hits the vertical borders, its following trajectory is no longer considered. Finally, we assume that each particle behaves independently from the others and that there is no fission or fusion of particles.

In the sequel, we observe the dynamics at discrete times ∆t, 2∆t, 3∆t . . . We denote ∆t the time step on the subset [-l, 0] × [0, H] with l < L. The observations are recorded during a time interval [0, T S ]. As we suppose that a particle does not change its drift direction along its trajectory, we assume that v x > 0, even though particles can actually move in both directions, which requires a classification to separate them into two groups. We consider that an observed tracklet of a given trajectory is an output if the last observed point of the segment is within a neighborhood of {0} × [0, H]. Meanwhile, we consider that it is an input if the first observed point is within a neighborhood of {-l} × [0, H]. Our main objective is then to associate the set of tracklets exiting the observed set [-l, 0] × [0, H] with the set of tracklets entering this observation set. The challenge is to correctly match the outputs and the inputs associated to particles (see Fig. 1).

PROBABILISTIC MODELS AND METHODS

Let us consider a given sample S, the observation set of all the trajectories. Define the sets O S = {o 1 , ..., o p } and I S = {i 1 , ..., i q } of p outputs and q inputs. Each output o = (t o , y o ) ∈ O S is characterized by its output time t o and its position y o ∈ [0, H] where the particle left the observed region. Similarly each input i = (t i , y i ) ∈ I S is characterized by its input time t i and its position y i ∈ [0, H] where it entered the observed region. A particle "involved" in an output o ∈ O S either died after time t o in the unobserved region, or is "involved" in a given input i ∈ I S with t i > t o . We will denote this event by {o → i}. Similarly a particle "involved" in an input i ∈ I S was either born before time t i in the unobserved region, or is "involved" in a given output o ∈ O S with t i > t o , which corresponds to the event {o → i}. 6) and ( 7))

Define c = (D c , B c , b c ) with D c ⊂ O S , B c ⊂ I S
Computation of the probability P (Dc) that the outputs in Dc died in the unobserved region. (section 3.1 Eq. ( 5))

Computation of the probability P ( c) that a set of outputs are connected to a set of inputs. (section 3.1 Eqs. ( 8) and ( 9) Likelihood (section 3.1 Eqs. ( 4), ( 10)-( 12)) Maximization of likelihood with CPLEX, (section 3.2 Eqs. [START_REF] Hussain | Mreb filaments align along greatest principal membrane curvature to orient cell wall synthesis[END_REF] and( 14))

Tracks reconstruction and estimation of connection accuracy, (section 4.2)

(+) Lifetime of a particle T d ∼ E(τ d ) (+)
The first passage time on l of a Brownian motion T l follows an Inverse Gaussian distribution (Prop 1).

(+) To emphasize, τα is a novel estimator Eq. [START_REF] Tweedie | Inverse statistical variates[END_REF] (+) The possibility to find the n th optimal solution, second part of section (3.2)

Figure 2: An outline of the connection procedure: from the estimation of the parameters to connection accuracy measurement, including likelihood formulation. All notations are defined in the corresponding sections.

Likelihood of a configuration

In this section, our objective is to derive an analytic expression of the likelihood Q(c) of a configuration c. The aim is to find, for a given sample S, the configuration ĉ such that P (ĉ/S) is maximal. It is difficult to calculate directly P (ĉ/S). Since c ⊂ S ⊂ O S , we can compute P (ĉ/S) working conditionally on O S .

However, since the model is in continuous time and involves random variables with continuous densities with respect to the Lebesgue measure, the conditional probability P (c/O S ) is equal to 0. This prevents to compute directly P (ĉ/S) with the classical conditional formula Therefore, for each input i = (t i , y i ) ∈ I S , we consider a spatiotemporal neighborhood

V i = T i × H i with T i = [t i -2 , t i + 2 ] and H i = [y i -2 , y i + 2 ] for some > 0.
The idea is to replace a given configuration c by a set C c of configurations where each element c * ∈ C c is similar to c but each input i ∈ I S is replaced by an input in V i . Formally, for each configuration c leading to the input set I S , C c is the set of configurations defined as follows:

c * = (D c * , B c * , b c * ) ∈ C c if and only if for each i ∈ I S , there exist i * ∈ V i satisfying:      D c * = D c , B c * = {i * , i ∈ B c }, For each i ∈ I S \ B c , b c * (b -1 c (i)) = i * .
With this definition, we have

P (c/S) = lim →0 P (C c /S) = lim →0 P (C c /O S ) c ∈C S P (C c /O S ) . (2)
In what follows, we study the behavior of P (C c /O S ) when goes to 0. We will always work conditionally on the realization of the output set O S but we will keep this conditioning implicit and write P (C c ) instead of P (C c /O S ) in order to simplify the notations. The study of P (C c ) will involve the probability for a particle to die in the unobserved region but also the probability that a particle born in this unobserved region enters the observed one in a given spatiotemporal neighborhood V i .

Furthermore, we assume that the particles born in the unobserved region, enter the observed one with a constant rate τ α and with a uniform distribution on {-l} × [0, H]. This is consistent with the fact that the particles are born with constant rate λ and appear uniformly at random on the membrane surface. Therefore, denote by N α the Poisson process of intensity τ α counting the number of inputs involved by particles born in the unobserved region.

Consider an output o ∈ O S and the possibility for the particle involved in o to die in the unobserved region. We have the following proposition (see [START_REF] Schrodinger | Zur theorie der fall-und steigversuche an teilchen mit brownscher bewegung[END_REF], [START_REF] Tweedie | Inverse statistical variates[END_REF], [START_REF] Wald | Sequential analysis[END_REF]). Proposition 1 Given the particle motion model as Brownian motion with drift as described in equation 1, the first passage time noted as T l on the entrance line {-l}×[0, H] of a particle starting at position z 0 = (0, y 0 ) for some y 0 ∈ [0, H], follows a law of inverse Gaussian, that is T l ∼ IG lu vx , lu σx 2 where l u := L -l is the length of the unobserved region.

Recall that if X ∼ IG(µ, λ), then X ≥ 0 almost surely, and for each x ≥ 0,

P (X ≤ x) = x 0 λ 2πy 3 exp - λ(y -µ) 2 2µ 2 y dy. (3)
In our framework, the event corresponding to the death of a particle with life duration T d following an exponential law of parameter τ d in the unobserved region is precisely {T d < T l }. Hence, we can derive an explicit expression of P (C c ).

Assume small enough so that for each i, i ∈ I S , T i ∩ T i = ∅. For a given configuration c and a given > 0, we will write

C c = (D c , B c , c ) with B c = {B c * , c * ∈ C c } and c = {b c * , c * ∈ C c }.
Due to the independent behavior of the particles, we have the following decomposition:

P (C c ) = P (D c )P (B c )P ( c ). (4)
We can then compute separately the probabilities of events D c , B c and c . First, note that we can assume without loss of generality that each output o ∈ D c starts at time t o = 0 and that only the position y o ∈ [0, H] fluctuates with o, but with no influence on T d or T l . Moreover, the loss of memory property of the exponential law ensures that the life duration T d of the particle after the output o still follows an exponential law of parameter τ d .

Since all outputs behave identically and independently, we have P (D c ) = P (T d < T l ) |Dc| ,where |D c | stands for the cardinal of D c . According to proposition 1, and since T d and T l are independent, we have

P (T d < T l ) = +∞ 0 t l 0 f T d (t d )f T l (t l )dt d dt l , (5)
= +∞ 0 t l 0 τ d e -τ d t d l u σ x 2πt 3 l exp - (v x t l -l u ) 2 2σ 2 x t l dt d dt l , = +∞ 0 l u (1 -e -τ d t l ) σ x 2πt 3 l exp - (v x t l -l u ) 2 2σ 2 x t l dt l ,
where f T d and f T l respectively stand for the density functions of T d and T l .

Now, consider the event B c . We call "spontaneous input" an input related to a particle born in the unobserved region that has never been observed. The set B c is defined so that, for each input i ∈ B c , we have exactly one "spontaneous input" appearing during the time interval T i , with a position in H i . Moreover, outside ∪ i∈Bc T i , there is no "spontaneous input". Formally, we have

B c = N α [0, T S ] \ i∈Bc T i = 0 i∈Bc {N α (T i ) = 1} ∩ H i , (6)
where N α is a Poisson process of intensity τ α associated to the counting of inputs involved by particles born in the unobserved region on the time interval [0, T S ]. In order to simplify the notations, H i denotes also the event of "spontaneous" appearance of an input i in H i . This event is independent of the process N α , and since the "spontaneous inputs" appear uniformly on [0, H], we have

P (H i ) = H .
Meanwhile, for any time interval I, N α (I) follows a Poisson law of parameter τ α |I| where |I| denotes the length of the interval I. Since is small enough so that for each i, i ∈ I S , T i ∩ T i = ∅, N α (T i ) and N α (T i ) are independent. Consequently, we can compute P (B c ) as follows:

P (B c ) = e -τα(T S -|Bc| ) τ α e -τα H |Bc| = 2 τα H |Bc| e -ταT S . (7)
Finally, consider the event c . For each input i ∈ I S \ B c , we denote by {o i c → V i } the survival event of the particle involved in the output o i c = b -1 c (i) in the unobserved region which appears in the spatiotemporal neighborhood V i . Since the particles behave independently, we have

P ( c ) = i∈I S \Bc P {o i c → V i } . (8)
In the sequel, we consider a given input i ∈ I S \ B c and its related output o = b -1 c (i). Defining s i = t i -t o and h i = y i -y o allows us to center the situation around the output o in the following way. A particle born at time 0 in position z 0 = (0, 0) has a life duration T d following an exponential law of parameter τ d . During its lifetime, the position of the particle is driven by a Brownian motion with drift Z t = (X t , Y t ): Z t = vt + ΣB t , where B t is a two-dimensional Wiener process and v and Σ are given in Equation ( 1). Define T l the first reaching time of l u = L -l of the process X t . The event {o → V i } can now be written as follows:

{o i c → V i } = {T d > T l } T l ∈ s i - 2 , s i + 2 Y T l ∈ h i - 2 , h i + 2 .
This expression corresponds exactly to the fact that in order to realize {o i c → V i } the particle needs to have a life duration longer than its first reaching time of l u and to appear in the spatiotemporal neighborhood s i -2 , s i + 2 × h i -2 , h i + 2 . Furthermore, T d follows an exponential law of parameter τ d , Y t follows a Gaussian law of parameters v y t and σ 2 y t and T l ∼ IG lu vx , lu σx 2 . Moreover, due to the fact that Σ is diagonal, the process Y t is not only independent of T d but also of T l . This allows us to write

P {o i c → V i } = s i + 2 s i -2 f T l (t l ) +∞ t l f T d (t d ) h i + 2 h i -2 f Yt l (y)dy dt d dt l .
As the two integrals involve a small domain of size , P

({o i c → V i }) ∼ O( 2 )
, and

lim →0 P ({o i c → V i }) 2 = f T l (s i )f Ys i (h i ) +∞ s i f T d (u)du (9) = l u σ x 2πs 3 i exp - (v x s i -l u ) 2 2σ 2 x s i 1 σ y √ 2πs i exp - (h i -v y s i ) 2 2σ 2 y s i e -τ d s i = l u 2πσ x σ y s 2 i exp - (v x s i -l u ) 2 2σ 2 x s i - (h i -v y s i ) 2 2σ 2 y s i -τ d s i .
For each configuration c, we calculate the likelihood Q(c) of the configuration c as follows:

Q(c) := lim →0 P (C c ) 2|I S | .
From (4) and Equations (5, 7, 8 and 9), we finally obtain the likelihood

Q(c) = τ α H |Bc| e -ταT S +∞ 0 l u (1 -e -τ d t l ) σ x 2πt 3 l exp - (v x t l -l u ) 2 2σ 2 x t l dt l |Dc| × i∈I S \Bc l u 2πσ x σ y s 2 i exp - (v x s i -l u ) 2 2σ 2 x s i - (h i -v y s i ) 2 2σ 2 y s i -τ d s i . (10)
Note that the limit when goes to 0 of P (C c ) 2|I S | is well defined, strictly positive, and that the exponent 2|I S | does not depend on the configuration c.

Recalling (2), this allows us to write

P (c/S) = Q(c) c ∈C S Q(c ) (11)
and as a consequence, we have

ĉ = argmax c∈C S {Q(c)}. (12)

Maximum likelihood and optimal configuration

The aim of this section is to identify the configuration c corresponding to the maximal likelihood Q(c) (see Equation [START_REF] Domínguez-Escobar | Processive movement of mreb-associated cell wall biosynthetic complexes in bacteria[END_REF]). Define

β := -log τ α H , δ := -log +∞ 0 l u (1 -e -τ d t l ) σ x 2πt 3 l exp - (v x t l -l u ) 2 2σ 2 x t l dt l
and for each configuration c and each i ∈

I S \ B c γ i c := -log l u 2πσ x σ y s 2 i exp - (v x s i -l u ) 2 2σ 2 x s i - (h i -v y s i ) 2 2σ 2 y s i -τ d s i . It follows that ĉ = argmax c∈C Q(c) = argmin c∈C -log (Q(c)) (13) = argmin c∈C   β|B c | + δ|D c | + i∈I S \Bc γ i c   .
This decomposition allows us to consider a linear optimization problem where β represents the cost of the spontaneous birth of an input, δ the cost of the death of an output and γ i c the cost of the connection between the output b -1 c (i) and the input i. The cost of connection can be defined for any couple (o, i) ∈ O S × I S as

γ i o := -log l u 2πσ x σ y s 2 o,i exp - (v x s o,i -l u ) 2 2σ 2 x s o,i - (h o,i -v y s o,i ) 2 2σ 2 y s o,i -τ d s o,i , where s o,i := t i -t o , h o,i = y i -y o and the convention γ i o = +∞ if t i ≤ t o .
In order to write in a canonical way this linear optimization problem, we associate to each configuration c a family of coefficients (c

o,i ) (o,i)∈O S ×I S such that c o,i = 1 if b c (o) = i and c o,i = 0 if b c (o) = i.
Since an output can be connected to at most one input, for each o ∈ O S , i∈I S c o,i ∈ {0, 1} and i∈I S c o,i = 0 corresponds to the death of the output o. Similarly, for each i ∈ I S , o∈O S c o,i ∈ {0, 1} and o∈O S c o,i = 0 corresponds to the fact that the input i is a "spontaneous input".

Our optimization problem is then equivalent to finding the family of coefficients (c o,i ) (o,i)∈O S ×I S that minimizes the quantity

β i∈I S 1 - o∈O S c o,i + δ o∈O S 1 - i∈I S c o,i + o∈O S i∈I S γ i o c o,i
or equivalently

K(c) := o∈O S i∈I S γ i o -β -δ c o,i s.t.      ∀o ∈ O S , ∀i ∈ I S , c o,i ∈ {0, 1}, ∀o ∈ O S , i∈I S c o,i ∈ {0, 1}, ∀i ∈ I S , o∈O S c o,i ∈ {0, 1}. (14)
In order to avoid to have infinite costs γ i o when t i ≤ t o ,, we can also impose c o,i = 0 if t i ≤ t o . Actually the problem ( 14) is a conventional linear optimization problem which can be solved by applying the CPLEX Linear Programming solver (https://www.ibm.com/analytics/cplexoptimizer).

The configuration ĉ is then the solution of the optimization problem [START_REF] Jaqaman | Robust single-particle tracking in live-cell time-lapse sequences[END_REF] and corresponds to the most likely configuration given the sample S. In order to complete the study, we propose to compute the following most likely configurations in a reccurent way by solving [START_REF] Jaqaman | Robust single-particle tracking in live-cell time-lapse sequences[END_REF] with additional constraints ensuring that the solution is different from the previous ones. In other words we define recursively the sequence (c n ) n∈N in the following way:

c 1 := ĉ ∀n ≥ 2, c n solves (14) with the n -1 additional constraints ∀k ∈ {1, . . . , n -1}, o∈O S i∈I S c o,i n (1 -c o,i k ) + (1 -c o,i n )c o,i k ≥ 1. ( 15
)
With this definition, c n is then the n-th most likely configuration. When n is greater than the number n S of configurations compatible with the sample S, the constraints are impossible to satisfy. In other words this sequence is well defined up to n S .

Estimation of parameters

Several parameters are involved in our computational approach. In this section, we propose clues to set these parameters. First, the parameters v and Σ can be estimated with classical maximum likelihood estimation procedures.

Second, we propose an estimator τd of τ d as explained below. The sample S can be considered as a set of points p = (t p , Z p ) observed at time t p and position Z p = (X p , Y p ) grouped in clusters s corresponding to tracklets of trajectories. The death of a particle in the observed region is detected in S for each point p ∈ S for which the associated tracklet s p has no successor point at time t p + ∆t. In order to be sure that the absence of successor is effectively due to the death of a particle and not to a particle leaving the observed region, we restrict the analysis to a region excluding a neighborhood of the border. However, we can check in this neighborhood the existence of successors for points in the restricted region.

We denote by S r ⊂ S the sample of points in the restricted region. For each point p ∈ S r , we denote by D p the event corresponding to the absence of successor for p. This corresponds to the fact that the particle involved in p died during the time interval [t p , t p + ∆t]. Since the life duration T d of a particle follows an exponential law of parameter τ d , and the absence of memory property of the exponential law, we have

P (D p ) = P (T d ∈ [0, ∆t]) = 1 -e -τ d ∆t . (16)
Hence, we define our estimator τd as

τd = 1 ∆t|S r | p∈Sr 1[D p ], (17)
where

Our estimator τd is then consistent as ∆t is small enough. Moreover, since the variables 1[D p ] are i.i.d Bernoulli random variables, we can calculate the related confidence interval. If q α denotes the α-quantile of the standard normal distribution, we have the following confidence interval of level α for 1-e -τ d ∆t ∆t :

CI α =   τd -q α τd 1 ∆t -τd |S r | , τd + q α τd 1 ∆t -τd |S r |   . ( 19
)
If ∆t is small enough, we get a good approximation of a confidence interval of level α for τ d since

lim ∆t→0 1 -e -τ d ∆t ∆t = τ d .
Now, we describe the estimation procedure for the rate τ α of "spontaneous inputs" induced by particles born in the unobserved region [-L, -l] × [0, H] and reached the border {-l} × [0, H]. We assume here that the parameters v, Σ and τ d are known, keeping in mind that in practice estimators are used instead. As introduced earlier, L is the perimeter of the cylinder, l is the length of the observed region, and l u = L -l is the length of the unobserved region. For a given length x, we denote by N x the number of tracklets born in the region In Fig. 3, N l = 2 correspond to tracks #1 and #4, S * l = {2, 5}, and the event B le 2 is realized while B le 5 is not.

Note that since the particles have the same independent dynamics, P (B x s ) does not depend on s. For x < l, this probability can easily be estimated as follows:

px = N x S o ,
where S o is the set of tracklets having an output in {0} × [0, H]. The strong law of large numbers yields a consistent estimator and allows us, in the case where l e < l, to define our estimator τα as follows:

τα = N l + ple |S * l | T S . (20)
Intuitively, this estimator amounts to counting the number of particles reaching {0} × [0, H] with weight 1 for each tracklet that we actually saw being born in the observed region and with weight ple for each spontaneous input that appeared in {-l} × [0, H]. Note that, as

N lu = N l + s∈S * l 1[B le s ],
τα is an unbiased estimator of τ α . Moreover, if we assume that the number of observed tracklets grows linearly with the observation time T S , this estimator is consistent when T S goes to +∞. Now, we consider the case l < l e < 2l which can easily be extended to the general case l < l e . Consider s ∈ S * l and denote for each interval J ⊂ [-L, 0] the event B J s where the tracklet s is born in the region J × [0, H]. The event B le s can be decomposed as follows:

B le s = B [-2l,-l] s B [-2l,-l] s ∩ B [-lu,-2l] s .
The loss of memory and homogeneity properties of the dynamics lead to the following estimator ple : ple := pl + (1 -pl ) ple-l .

Limits of the model

The main assumptions in this work are homogeneity in time and space, induced by the constant death and birth rates, as well as constant speed and noise. While these assumptions lead to a simple model and allows a reasonably technical study, it is natural to question it. The main reason of this choice is that it corresponds to uniform laws when we have no reason prioritize one specific behavior in particular.

Note that a similar study can be made with different speeds among trajectories. This can be done by classifying the trajectories according to their speeds and applying the present procedure to each class. This would lead to the same estimation procedure with smaller datasets but theoretical results will still hold.

We then discuss the homogeneity in time, for which the most questionable assumption is the constant death rate that could possibly depend on the position or on the age of the particle. Concerning the dependence in space, this modification would lead to the estimation of a function of the position instead of the simple constant τ d . From a practical point of view, this would increase the dimension of the parameter to estimate, with the same size of dataset. From a theoretical point of view a more technical study can be made as long as we assume the death function rate (depending on the position) constant on each tracklet {y}×] -L, 0] in order to overcome the issue of partial observation.

Concerning the dependence in time, the assumption that the death rate depends on the age of the particle prevents to propose a similar study. Indeed, due to partial observation, the age of each particle entering the observed region is unknown and can not be estimated.

Modeling hypothesis and MreB dynamics

The study of the dynamics of MreB patches or assemblies in the vicinity of the internal membrane of Bacillus subtilis bacteria reveals several subpopulations undergoing constrained, randomly or directionally moving [START_REF] Billaudeau | Contrasting mechanisms of growth in two model rod-shaped bacteria[END_REF]. Herein we are interested in the directionally moving subpopulation dynamics. This subpopulation moves possessively around the cell diameter [START_REF] Garner | Coupled, circumferential motions of the cell wall synthesis machinery and mreb filaments in b. subtilis[END_REF][START_REF] Domínguez-Escobar | Processive movement of mreb-associated cell wall biosynthetic complexes in bacteria[END_REF]. Following Hussain et al [START_REF] Hussain | Mreb filaments align along greatest principal membrane curvature to orient cell wall synthesis[END_REF], Billaudeau et al [START_REF] Billaudeau | Mreb forms subdiffraction nanofilaments during active growth in bacillus subtilis[END_REF] confirmed that directionally moving filaments travel in a direction close to their main axis, perpendicularly to the long axis of the cell (angle γ = 89.9 • ± 37.0 • ). Hence, for some filaments, the speed vector may have a component in the main direction of the bacteria.

According to Wong et al. [START_REF] Wong | Mechanics and dynamics of translocating mreb filaments on curved membranes[END_REF] a motion model (named "biased random walk") reproduces the dynamics patterns of MreB filaments. In their simulations, the speed is constant, the noise variance between several time steps depends on the duration and, possibly on the local curvature of the surface. These properties are shared with the Brownian motion model with constant drift we consider.

SIMULATION STUDY

In this section, we present a series of experiments performed on synthetic datasets. These experiments aim to evaluate and analyze the sensitivity of the reconstruction procedure when the characteristics of the dynamics as well as the spatio-temporal sampling resolution of observations vary. In addition to demonstrate the potential of our procedure, these experiments might also be useful for the design of the experimental setting for images acquisition. The reconstruction procedure has been implemented in MATLAB ver. R2018b. The codes are available on Github https://github.com/atrubuil/ReconstructionOfTruncatedTrajectories.

Generation of trajectories

Trajectories are generated on a rectangular unwrapped cylindrical surface of size [0, L] × [0, H] (Fig. 4). In our experiments, we set L = 50, H = 30. The initial position of each trajectory is drawn from uniform distribution on the surface. Time duration T between two births follows an exponential law with birth rate parameter λ. At each birth, the intrinsic properties of a trajectory i are given, such as velocity v i , variance Σ i , and lifetime T i d . The lifetime T d follows an exponential law, with the same death rate τ d for all trajectories in the whole simulated image sequence. The drift v i = (v xi , v yi ) and noise Σ i = σ xi 0 0 σ yi are set to be constant along one given trajectory.

According to the assumptions made on real biological context, unless otherwise stated, it is set by default, θ = 0.01(≈ 0.6 • ) is the angle between the direction of motion of particle and the X direction, v y = tan(θ)v x , σ x = σ y = σ, v x = 0.6, σ = 0.2, λ = 0.03, τ d = 0.005. The time interval between two images ∆t = 0.25. As known, the theoretical depth of the observation field of TIRFM is 200nm, the diameter of the bacteria cell is 1µm, therefore the width of the ROO l is set to 14.76 and that of the unobserved region l u = 35.24 (unit in pixel, note that in TIRF images 1 pixel ≈ 64nm).

As there is no particle on the surface at the beginning, the simulated set of trajectories needs some warm-up time to reach the stationary regime, where the law of the number of trajectories does not depend on time. The assumed dynamic process is a special case of birth and death process. As a known result [START_REF] Karlin | A first course in stochastic processes[END_REF], the expectation of the trajectories number N during stationary regime is E(N ) = λ τ d . To ensure that the dynamics are in a stationary regime, the images sequence is simulated long enough, for around 2 hours (Fig. 5).

The "Adjusted Rand Index" for the evaluation of connection results

Given the true and estimated class assignments, we compute the so called Adjusted Rand Index to evaluate similarity or consensus between the two sets. The Adjusted Rand index is the corrected-for-chance version of the Rand index. It is scored exactly 1 when the two sets are identical, close to 0 for random labeling. It could be negative when the index is lower than the expectation under random labeling. More precisely, let G and K be the true and estimated assignments respectively, let us define a and b as: a the number of pairs of elements that are in the same class in G and in the same class in K, b the number of pairs of elements that are in different classes in G and in different classes in K. The raw (unadjusted) Rand index is then given by:

RI = a + b C M 2 , ( 21
)
where C M 2 is the total number of possible pairs in the dataset (without ordering) of size M . The RI score does not guarantee that random assignments will get a value close to zero. This is especially true if the number of clusters has the same order of magnitude as the number of samples. To overcome this difficulty, we prefer to consider the Adjusted Rand Index defined by [START_REF] Rand | Objective criteria for the evaluation of clustering methods[END_REF]:

ARI = RI -E(RI) 1 -E(RI) . ( 22
)
Here E(RI) denotes the expectation of the Rand Index where the estimated assignment K is replaced by an assignment chosen uniformly at random. This means that the assignment procedure does not do better than random assignment if the ARI score is zero, and that it does worse than random if ARI < 0.

Experimental results

The good performance of the connection procedure relies on the estimation of the characteristics of the dynamics: the speed, v, the diffusion variance, Σ, the arrival rate τ α and the death rate τ d , as these quantities are used in the calculation of the likelihood (Eq. 13). Here we evaluated the impact of spatio-temporal sampling (l/l u , T S ) on the estimators and the impact of parameters of the dynamics (v, Σ, τ α , τ d ) on the accuracy of the reconstruction. The estimator τα is proposed in Eq. 20. Here we test how it performs with different spatiotemporal sampling (l/l u , T S ), and different birth rate λ and death rate τ d . By its definition in section 3.1, τ α , the rate of "spontaneous input" induced by particles born in the unobserved region and reach the border of the ROO, is not a preset parameter. A reference of the "true" value of τ α is given by

N lu
Ts , where N lu denotes the number of tracklets born in the region ] -l u , 0] × [0, H] and reached the border {0} × [0, H], l u is the width of the unobserved region.

Next, we test the robustness of the estimator τα w.r.t. l/l u (Fig. 6). To avoid the influence of T S on the consistency of the estimator, T S is set to be long enough as 30 min. We can conclude that, naturally, the more the observed area is larger, better is the performance of the estimator τα . In the case of the simulation of the real situation, where l/l u = 0.42, the estimator works reasonably good.

Following, we test the robustness of τα w.r.t. T S (Fig. 7). This test is essential because in reality it is impossible to use a 30-min movie, because of technical issues like photobleaching of fluorophores and natual growth in living samples. At this stage, the propotion of observed and unobserved region l/l u is set to 0.42. T S varies from 2.5 min to 30 min. In Fig. 7, it can be noticed that the reference 'ground truth' of τ α (blue boxes) decreases as T S lengthens. Actually, the reference is only a pseudo 'ground truth'. It is sensible to T S when T S is small The absolute value of τ α depends on λ and τ d . Fig. 8 displays for different combinations of λ and τ d , the estimations of τα by 5-min movies (magenta) and 30-min movies (blue). It shows that τ α increases linearly as the birth rate λ increases, and decreases slightly linearly as the death rate τ d increases.

The estimator τd is unbiased and performs reasonably well with 5-min movies

As explained in section 3.3, τd is a rather classical estimator. Fig. 9 shows the estimator with 5-min movies (magenta) and 30-min movies (blue) respectively. It confirms that the estimator is unbiased. Black horizontal lines represent the true value of τ d . Naturally, the variance is bigger with shorter movies.

The choice of T S

According to Figs. 8 and9, when T S = 30 min, the estimators of τ α and τ d perform well, being converged with small variance. As 30-min movie acquisition is almost infeasible under the situation of fluorescence microscopy, we need to find a compromise with smaller T S and The estimator is unbiased. The variance of the estimator is larger with shorter movies (magenta). For a given τ d , when birth rate λ increases (e.g. the first four boxes), then the number of particles also increases, in consequence, the variance decreases.

reasonably good estimators. We tested especially T S =2.5 min and T S =5 min. Comparing the estimation results with 2.5-min movies, we found that T S = 5 min is a good choice to limit the estimation error of τ d and to ensure a good connection performance (more details about the experiments for the choice of T S in Supplementary Materials 1).

4.3.4

The connection procedure works well, even when true parameters are unknown

In this part, we assess the performance of the connection algorithm with different parameters λ and τ d . We evaluate as well the impact of the error of the estimator, by using in the connection procedure respectively true parameters τ α , τ d , v, σ and their estimators τα , τd , v, σ. The duration of movies T S is set to 5 min. The connection results measured by ARI are presented in Fig. 10.

Each pair of blue and magenta box represents the connection result of a setting of λ and τ d . The black line represents the mean value of the number of tracklets fluctuating with different settings of λ and τ d . The performance of connection is affected by the number of tracklets in each movie to be connected. The higher the density of tracklets is, the more difficult it is to find the right ones.

It can be noticed that the ARI value when we use the estimators τd , τα , v and σ, is almost as good as when we use true values for all the parameters. This is an encouraging result as it means that it is feasible to apply the algorithm in real image sequences. When the number of tracklets is around 20 (e.g. λ = 0.04 and τ d = 0.008), the median values of ARI are higher than 0.9, showing a promising connection performance. Even for the case with the highest particle density, when the average number of tracklets reaches 100 (λ = 0.1 and τ d = 0.004), the median value of ARI is still higher than 0.7.

4.3.5

The connection procedure is robust even when each particle moves at different speed (but with constant speed along a trajectory). However v and σ should be well estimated

In the previous experiments, all trajectories are generated with the same speed v and standard error σ. In this section, we design experiments to test the performance of the connection algorithm, when the drift v varies from particle to particle, v x ∼ U(0.5, 0.9). In one movie, as all particles are in the same environment, there is no obvious reason for different particles to have different σ . Therefore the standard error σ is set to be constant for particles in one movie. However, we test in independent movies, when σ = 0.2, 0.3 or 0.4, the influence of σ on the performance of connection procedure. Other parameters to be specified are the angle between the direction of the motion and the circumferential direction of the cylinder, θ = 0.15(≈ 8.6 • ), v y = tan(θ)v x , σ y = tan(θ)σ x , σ x = σ, birth rate and death rate are fixed, with λ = 0.08 and τ d = 0.02.

The normalized error (NE) of an estimator is defined by the error of the estimator normalized by its ground truth. For example, the NE of v x equals to vx-vx vx . In Fig. 11, the NEs of vx , vy , σx and σy when σ takes different values are presented. It shows that when σ increases, the variance of vx and vy increases.

For tracklets connection, we compare the results when true values of v and σ or when the estimated value v and σ respectively are taken by the connection procedure. The experiments are carried under three situations, when σ = 0.2, 0.3, and 0.4. The results in Fig. 12 shows that whether using true v and σ or estimated value v and σ, the performance measured by ARI degrades when σ increases. When the standard error σ = 0.2, using true v and σ, the median value of ARI reaches to 1. When using the estimated v and σ, the median value of ARI is approximately 0.75. It can be concluded that the estimation of v and σ has an impact on the performance of the algorithm.

Analysis of the connection results

An example of tracklets connection

Figure 13 shows, on the left, trajectories in a movie and on the right, the results of tracklets connection. The path from an output to the matched input is represented by the dashed straight line, as we don't know how exactly the particle went through the hidden zone. The only wrong connection corresponds to the bold line. Compared with the figure on the left, we can find the realization of these two tracklets. In reality, the orange bold tracklet disappeared at the hidden region and the bold purple tracklet appeared nearby and entered into the observed zone.

In fact, not only the optimal configuration can be calculated, but also the most likely alternative configurations in decreasing order of probability (Fig. 14). It should be noticed that the optimization algorithm tries to minimize K(c) = -log Q(c), instead of finding the c * maximizing P (c). It costs too much to obtain the probability P (c * ), as it requires the enumeration of all the possible configurations c ∈ C (Eq. 11). However, The number of configurations can be determined in order to guarantee that the sum of the probability of these configurations will be greater than given threshold (see Supplementary Materials 3). As a result, we can upper bounds for the probability.

For example, we see that the second most likely the realization of trajectories, 0.187 < P (c t ) < 0.205 according to the algorithm (Fig. 14).

Combining with Fig. 13, the optimal configuration found by the algorithm, committing one connection error, does not correspond to the realization. In section 2, we evaluated the connection error caused by randomness.

The number of rotations around the cylinder

Once the connection procedure is achieved, we can address the question of the number of rotations of a particle around the cylinder. In the context of simulation, the death rate τ d and the dynamic velocity v x are known. Accordingly, the value of the number of rotations is known to be equal to

CONCLUSION

In this paper, we proposed a probabilistic framework and a computational approach with no hidden parameter to connect tracklets from 2D partial observations. We provided several consistent estimators of parameters to automatically drive the connection procedure. The performance of our procedure is satisfying if we consider the ARI criterion. Moreover, an ordered set of the best reconstructions could also be proposed. The robustness of the procedure has been tested for different drifts, diffusion of the dynamics, and trajectory densities. Our computational approach can be extended to the case when the drift/speed is not the same for all particles but remains constant along time. In that case, it is straightforward to estimate and classify the drifts before applying our connection procedure to each class of drift since the tracklets with different speeds are not likely to be connected.

After recovering the whole trajectory on the surface of the cylinder, we can have a better understanding of the average duration of a particle, and more accurate statistics about the spatio-temporal organization of particles. The simulation study can also serve as a guideline for the design of experiments.

The connection procedure is tested with a real TIRFM dataset. The experimental results are illustrated in Appendix A. For future works, we plan to investigate more on real TIRFM datasets. Experiments on real data show that the observed region corresponds approximately to one-third of the total surface, which is rather small. However, we have shown that we are able to cope with the hidden region of such size. Nevertheless, several assumptions and approximations need to be further investigated. For instance, we assumed spatial homogeneity, suggesting that the particles are born or die uniformly on the membrane surface. Moreover, we assumed a memoryless lifetime while dependency with respect to particle "age" could be more realistic.

A AN ILLUSTRATION OF THE CONNECTION ALGORITHM APPLIED TO REAL MREB DYNAMICS

Data obtained using TIRF microscopy of MreB aggregates in Bacillus subtilis ( [START_REF] Billaudeau | Contrasting mechanisms of growth in two model rod-shaped bacteria[END_REF]) are considered. A typical movie from this dataset shows several MreB aggregates moving inside one or several cells (see Fig. 1 and Supplementary Materials 3). The pixel size, frame rate and duration are respectively ∆x = ∆y = 64nm, ∆t = 1s, T = 2mn. Hereafter, we selected one cell to illustrate the application of our algorithm. First, tracklets exhibiting directed motion should be extracted from the movie data, then tracklets should be projected back on the cylinder shape of the cell and unwrapped, eventually the connection algorithm is applied and a list of likelihood decreasing ordered configurations of trajectories connections is presented to the user.

A.1 Construction of the local cell referential

Once MreB aggregates pixels are separated from the background inside each image of the movie, a bounding box is drawn around a given cell and a local x-y referential is estimated using Principal Component Analysis (PCA) on the coordinates of pixels belonging to aggregates (Fig. 16). The z coordinate of an aggregate is inferred using as a prior the cylinder shape of the bacteria and its radius, R, so z(x, y) = R -√ R 2 -x 2 . Using U-Track [START_REF] Jaqaman | Robust single-particle tracking in live-cell time-lapse sequences[END_REF], MreB aggregates are tracked and constitute a set of tracklets. The automatic classification of these tracklets in three classes, respectively Brownian, subdiffusive and, directed motion is done using two algorithms: the classical MSD algorithm and a recent algorithm ( [START_REF] Briane | Statistical analysis of particle trajectories in living cells[END_REF]). The tracklets classified as directed motion by either one of the two algorithms are selected for the application of our connection algorithm (Fig. 17). The tracklets were projected back on the cylinder and unwrapped, as explained in the technical part of the paper. As we can see, only a few aggregates crossed the borders of the visible region. Others aggregates, according to our definitions are born or die in the visible region, which is not true. When an aggregate approaches the borders, its intensity becomes weak as it is farther from the support plane, and less excitation light penetrates higher z-position in TIRF microscopy settings. As a result, the detection algorithm fails to detect the aggregates when they approach the borders.

A.3 The connection of tracklets

All the selected tracklets crossing the magenta lines in Fig. 17 First, the speed and diffusion are estimated (Fig. 18) for each tracklet, respectively. Two populations of tracklets evolving in opposite x directions are identified. These two populations are considered one after the other in the connection procedure. Tracklets corresponding to speed lower than 0.4 are filtered out. For the population of tracklets associated with positive (resp. negative) v x , death rate τ d is estimated as 0.0691 (resp. 0.0756). The arrival rate τ α is estimated as 0.0310 (resp. 0.0220). † tracklets of positive speed v x The first, fifth, seventh and eighth optimal configuration suggests one connection. The second suggests that there is no connection. The third, fourth and sixth configurations suggest two connections. Some of these configurations are shown in Fig. 19. † tracklets of negative speed v x

The first optimal configuration suggests no connection. The second one suggests one connection Fig. 20. In Fig. 21 we show a 3D reconstruction of the aggregates and two tracks that could correspond to aggregates more than loop around the cylinder surface of the cell. For the positive (resp. negative) speed set of tracklets, the eighth (resp. second) optimal solution was selected.

Figure 1 :

 1 Figure 1: (a): Several consecutive images from a real TIRFM movie[1]. Tracks are superposed on the images.(b) left: Illustration of trajectories observed during recorded time [0, T S ] on the surface of a cylinder. Only the motions inside the ROO ] -l, 0[×[0, H] can be observed, even though the dynamics happen on the whole surface; right: Representation of the dynamics on a 2D unwrapped surface ] -L, 0[×[0, H]. The objective is to recover the dynamics on the whole surface from the partial observations, by coordinating the inputs through {-l} × [0, H] and the outputs through {0} × [0, H] in a movie during T S , taking into account particles birth and death events.

  and b c a bijection from O S \ D c to I S \ B c in order to describe the configuration for which all outputs in D c died in the unobserved region, all inputs in B c are born in the unobserved region, and the event o∈Os\Dc {o → b c (o)} was realized. Our aim is to determine the maximum likelihood configuration c given the sample S. The outline of the connection procedure is given in Fig. 2, to facilitate the understanding of the modeling steps.Sample S, the outputs and inputs sets O S and I S , (section 3.1) Estimation of parameters τα, τd , v and σ, (sections 3.1, 3.3) Computation of the probability P (Bc) that the inputs in Bc are born in the unobserved region. (section 3.1 Eqs. (

P

  (c/S) = P (c/O S ) P (S/O S ) , because it gives P (S/O S ) = c∈C S P (c/O S ) = 0.

  |S r | stands for the number of points in S r and 1[•] denotes the indicator function. Due to the absence of memory property of the exponential law, the random variables 1[D p ] are i.i.d. As |S r | goes to +∞, the strong law of large numbers yields to lim |Sr|→∞ τd = 1 -e -τ d ∆t ∆t a.s. The justification of this choice for τd relies in the following almost sure convergence: lim ∆t→0 lim |Sr|→∞ τd = τ d a.s.

Figure 3 :

 3 Figure 3: An artificially constructed zone ] -l u , 0] × [0, H] having the same size as the unobserved region ] -L, -l] × [0, H]. The observed region is ] -l, 0] × [0, H] ; as the width of the invisible part is l u , the extended zone has width l e = l u -l.

Figure 4 :

 4 Figure 4: A set of simulated trajectories during 2.5 minutes (in stationary regime). X(resp. Y) axis represents the unfolded circumferential (resp. main) direction of the cylinder. Colors from light to dark green represent time evolution. Shadowed area corresponds to the unobserved region and white area corresponds to the ROO.

Figure 5 :

 5 Figure 5: Fluctuations of the number of trajectories w.r.t. time. At around t = 20 min, the trajectories number fluctuates around the theoretical expectation value 6.

Figure 6 :

 6 Figure6: At l/l u = 0.42, which corresponds to the realistic situation in 2D-TIRF, more than half of the trials presents a relative error smaller than 10%. The proposed estimator τα is unbiased and the variance decreases as l/l u increases.

Figure 7 :

 7 Figure7: The counted 'ground truth' τ α and the estimated τα obtained by movies of different duration, varying from 2.5 min to 30 min. Blue boxes (resp. Magenta boxes) correspond to the 'ground truth', i.e. the counted value (resp. the estimator τα ). The blue horizontal line represents the 'ground truth' value when T S = 30 min.

Figure 9 :

 9 Figure9: τd with different λ and τ d . The estimator is unbiased. The variance of the estimator is larger with shorter movies (magenta). For a given τ d , when birth rate λ increases (e.g. the first four boxes), then the number of particles also increases, in consequence, the variance decreases.

Figure 10 :

 10 Figure 10: Connection performance comparison for different λ and τ d , when T S = 5 min. Blue (resp. magenta) boxes represents ARI values obtained with true parameters τ d , τ α , v and σ (resp. with estimators τd , τα , v and σ). The black line represents the mean number of tracklets in a movie.

4 Figure 11 :

 411 Figure 11: The NEs of v x , v y , σ x and σ y in cases where σ = 0.2, 0.3 and 0.4.

Figure 12 :

 12 Figure 12: Comparison of connection performance measured by ARI, when the procedure takes true v and σ or estimated v and σ, in three experiments where σ = 0.2, 0.3 and 0.4 respectively.

Figure 13 :

 13 Figure 13: Left: Trajectories in one movie; Right: the connection results.

Figure 14 :

 14 Figure 14: The probability of n th optimal configuration and the probability of the realization.

  vxT d L , where T d ∼ E(τ d ) ensures a theoretical expectation value of vx τ d L . By counting the tracklets for each trajectory, we can obtain a proxy of the number of rotations around the cylinder. In Fig. 15, λ is set to 0.04 and different values of τ d between 0.004 and 0.01 are evaluated. Blue bars represent the distribution of the number of rotations of true connections. The magenta bars to display the distribution of the number of rotations estimated by the connection procedure. The corresponding ARI values, indicating the connection accuracy, are given as well. The density of the theoretical values of the number of rotations is presented in green color. The vertical lines represent the median values of the corresponding distribution.Overall, when τ d is small, the median value of number of rotations is higher and the distribution has a heavier tail. In general, the distributions with all three colors are similar to each other.

Figure 15 :

 15 Figure 15: The density distribution of 'number of rotations'. Each row represents the result for a different τ d value.

Figure 16 :A. 2

 162 Figure 16: The estimation of the local x -y referential for a cell.

  are considered.

Figure 17 :

 17 Figure 17: The tracklets classification. (a) MSD classification. (b) STP classification. Brownian tracklets (blue), Directed tracklets (red). Blue lines represent the border of the visible region. Magenta lines represent the 0.1 quantile and the 0.9 quantile of x coordinate values.

Figure 18 :

 18 Figure 18: The distribution of drift and variance in the selected tracklets population. (a) without filtering. (b) after filtering.

Figure 19 :

 19 Figure19: Some optimal configurations for the population of positive speed v x . From left to right, first, third and sixth better configurations. Connected tracklets are drawn with the same color.

Figure 20 :

 20 Figure 20: Second optimal configuration for the population of negative speed v x .

Figure 21 :

 21 Figure 21: The 3D reconstruction of tracks. The centroids of aggregates are represented as blue squares.The arrows indicates the direction of motion. The full blue lines represent tracklets crossing the observed region. The dotted red lines represent the simplified extrapolation of aggregate motion in the invisible region.

  Figure8: The estimation of arrival rate τ α w.r.t. different λ and τ d . For example, when λ = 0.04, τ d = 0.004, the median value of τα is around 0.025, which means that at each moment, the probability that a particle born in the invisible zone arrives at {-l} × [0, H] is around 0.025.

	0.16
		Estimated by a long film of around 30 min
	0.14	Estimated by a short film of 5 min
	0.12
	0.1
	0.08
	0.06
	0.04
	0.02
	0
		4/0.4 6/.4 8/.4 10/.4 4/0.6 6/.6 8/.6 10/.6 4/0.8 6/.8 8/.8 10/.8 4/1 6/1 8/1 10/1
		birth rate / death rate	10 -2
	0.03
		Estimated by a long film of around 30 min
		Estimated by a short film of 5 min
	0.025	True value
	0.02
	0.015
	0.01
	0.005
	0
		4/0.4 6/.4 8/.4 10/.4 4/0.6 6/.6 8/.6 10/.6 4/0.8 6/.8 8/.8 10/.8 4/1 6/1 8/1 10/1
		birth rate / death rate	10 -2

ACKNOWLEDGEMENTS

The authors thank L. Tournier for fruitful discussions about the optimization procedure, R. Carballido-López and C. Billaudeau for their inspiring work in MreB studies which triggered this research. This work was partially supported by ANR DALLISH, Programme CES232016.


Andrey Yu Zubarev

Maxim Raboisson-Michel

Gregory Verger-Dubois

Pavel Kuzhir

To the theory of ferrohydrodynamic circulating flow induced by running magnetic field

Keywords: Magnetic fluid, alternating magnetic field, field-induced flow

We present results of theoretical modelling of macroscopic circulating flow induced in a drop of ferrofluid by oscillating running magnetic field. The drop is placed in a narrow flat channel filled by a non-magnetic liquid. The aim of this work is development of a scientific basis for a progressive method of address drug delivery to thrombus clots in blood vessels with the help of the magnetically induced circulation flow. Our results show that the oscillating running field allows inducing the carrier fluid flow with velocity amplitude 1-10 cm/s. This is the range of values, presenting interest from the point of view of the drug delivery.

Introduction

Most common method of the treatment of brain stokes consists in injection of tissue plasminogen activator to dissolve the blood clots and restore the blood flux. However, the thrombolytic drug diffuses very slowly toward blood clots through blocked vessels. A smart solution of this problem has been proposed in [START_REF] Creighton | Magnetic-based systems for treating occluded vessels[END_REF]. The key idea of this method is injection of a ferrofluid with nano-sized magnetic particles in the blocked vessel and generation of recirculating flows under applied external magnetic field. Under the action of the applied field, the magnetic nanoparticles or their needle-shaped aggregates become to rotate; this rotation provokes recirculating flows in the blocked vessels which in turn considerably enhance convective transport of the drug towards the clots. Only a few works have been published on this topic [START_REF] Clements | A mathematical model for magnetically-assisted delivery of thrombolytics in occluded blood vessels for ischemic stroke treatment[END_REF][START_REF] Gabayno | Controlled manipulation of Fe 3 O 4 nanoparticles in an oscillating magnetic field for fast ablation of microchannel occlusion[END_REF][START_REF] Li | Thrombolysis Enhancing by Magnetic Manipulation of Fe3O4 Nanoparticles[END_REF]. However, the physical understanding of the origin of the recirculating or oscillatory flows is necessary for development and practical application of this method.

A theoretical analysis of the circulation flow in a liquid in a flat channel, a few moments after injection of a ferrofluid drop , under the action of a rotating field has been proposed in ref. [START_REF] Musickhin | Field-induced circulation flow in magnetic fluids[END_REF]. The considered ferrofluid was supposed to consist of individual spherical particles with permanent magnetic moments . It was shown that the flow appears only under the field with spatially nonuniform amplitude and only in the presence of nanoparticle concentration gradients (naturally induced by injection of a ferrofluid drop). Possible nanoparticle aggregation induced by applied magnetic fields observed in experiments and considerably enhancing recirculation flows [START_REF] Gabayno | Controlled manipulation of Fe 3 O 4 nanoparticles in an oscillating magnetic field for fast ablation of microchannel occlusion[END_REF][START_REF] Li | Thrombolysis Enhancing by Magnetic Manipulation of Fe3O4 Nanoparticles[END_REF] was not taken into account.

In this paper, we propose a modified theoretical model of the circulation flow induced in a drop of magnetic fluid, situated in a thin channel, modelling a blood vessel. The ferrofluid consists of elongated magnetizable particles mimicking field induced aggregates. Furthermore, we consider not rotating, but oscillating running magnetic fields. They can be generated using only one pair of Helmholtz coils and would allow a much more compact design of the final medical device than the rotating magnetic fields requiring two pairs of Helmholtz coils or mechanically rotating permanent magnets [START_REF] Creighton | Magnetic-based systems for treating occluded vessels[END_REF]. The obtained results demonstrate that the velocity of the field-induced recirculating flows in a closed channel is well above the minimal value required for efficient drug delivery in real blood vessels. This is quite encouraging result, and we believe that the brain stroke treatment can be realized with low magnetic fields/ electric currents, while the developed theoretical model could be useful for parametric optimization of this medical technique.

Mathematical model

We consider a flat channel of an infinite size in the y direction and with two opposite solenoids placed at the extremities along the horizontal axis x, as illustrated in Fig. 1. The distance 2Lx between the solenoids is supposed to be much larger than the channel thickness Lz. For simplicity, we suppose that inside the ferrofluid drop the viscosity 𝜂 is homogeneous and is the same as outside the drop. This means that at any moment of time, the volume concentration Φ of the particles inside the drop is relatively low, about a few volume per cent. A drop of the ferrofluid is situated in the middle, between the solenoids. The initial volume concentration Φ 0 of the particles is supposed to be known and depends only on the coordinate x. For the sake of simplicity, we consider here the two-dimensional model of the system, and suppose that all physical events take place in the plane x,z . Then, we consider the linear dependence of the ferrofluid magnetization M on the local magnetic field. Note that this approximation is not principal and cannot affect qualitative effects revealed by the present model. The "non-linear" generalization will make the calculations more cumbersome.

We suppose that magnetic fields H1 and H2, created by the solenoids 1 and 2, take the form:

𝐻 1𝑥 = ℎ 1𝑥 (𝑥, 𝑧)𝑐𝑜𝑠𝜔𝑡; 𝐻 1𝑧 = ℎ 1𝑧 (𝑥, 𝑧)𝑐𝑜𝑠𝜔𝑡 (1)
𝐻 2𝑥 = ℎ 2𝑥 (𝑥, 𝑧)𝑠𝑖𝑛𝜔𝑡; 𝐻 2𝑧 = ℎ 2𝑧 (𝑥, 𝑧)𝑠𝑖𝑛𝜔𝑡
where  is the angular frequency of the alternating magnetic field, t is the time and the spatial parts h(x,z) of the fields can be calculated by using the standard formulas for the fields of the cylindrical solenoids with given values of the current intensity, radius and length. Exact expressions for h(x,z) are provided in our previous work [START_REF] Musickhin | Field-induced circulation flow in magnetic fluids[END_REF]. Below we will use the notation 𝑇 = 2𝜋/𝜔 for the period of the field oscillation The results of calculations of the local absolute value H(x,z,t) of the total field H=H1 + H2 are illustrated in Fig. 2 for the realistic solenoid parameters indicated in the figure caption. We suppose that the ferrofluid consists of highly elongated magnetizable particles (mimicking field-induced aggregates composed of spherical ferrofluid nanoparticles). At this stage, we ignore the possible aggregation of the elongated particles as it was not reported in experimental studies [START_REF] Gabayno | Controlled manipulation of Fe 3 O 4 nanoparticles in an oscillating magnetic field for fast ablation of microchannel occlusion[END_REF][START_REF] Li | Thrombolysis Enhancing by Magnetic Manipulation of Fe3O4 Nanoparticles[END_REF].Estimates show that the time of the particle reorientation in the alternating field is much smaller than the period of the field oscillation with the typical angular frequency 𝜔~10 -20𝑠 -1 We can therefore suppose that the elongated particles are always oriented along the local magnetic field H. At such condition, the torque term (proportional to 𝐌 × 𝐇) disappears from the equation of ferrofluid motion. One can also show that for the considered channel dimensions and for oscillating flows with the typical frequency 𝜔~10 -20𝑠 -1 , the inertia terms can also be ruled out from the equation of motion. Under these approximations, equations of motion of the whole ferrofluid and of its solid phase (composed of elongated particles) take the following form [START_REF] Rosensweig | Ferrohydrodynamics[END_REF]:

𝟎 = -∇𝑝 + 𝜂Δ𝒗 + Φ𝑭 𝑚 (1a) 𝑑𝑖𝑣 𝒗 = 𝟎 (1b) 𝜕Φ 𝜕𝑡 + 𝑑𝑖𝑣 (Φ𝑭 𝑚 𝑉 𝑎 𝛽 ) + 𝑑𝑖𝑣(Φ𝒗) = 𝟎 (1c)
Here 𝒗 is the flow velocity, p is the pressure, Φ(𝑥, 𝑧, 𝑡) is volume concentration of the particles, Fm is the magnetic force acting by unit volume of each elongated particle by the magnetic field gradient, Va is the particle volume, 1/ 𝛽 is the component of the hydrodynamic mobility tensor of the elongated particle along the force 𝑭 𝑚 .

In the frames of the linear magnetization approximation, the magnetic force by unit volume of elongated particles can be estimated as follows: [START_REF] Landau | Electrodynamics of Continuum Media[END_REF] 𝐹 𝑚 = -∇𝑈, 𝑈 = -𝜇 0 2 𝜒 𝑎 𝐻 2 (2)

Here U is the potential energy of elongated particles oriented at each point along the local the field H(x,z,t) and evaluated under realistic approximation of high length-to-diameter ratio of particles, 𝜒 𝑎 is magnetic susceptibility of these particles, and 𝜇 0 = 4𝜋 × 10 -7 H/m is magnetic permeability of vacuum. . In the frames of these approximations, applying rot operator to the both parts of the first equation (1a), one gets:

∆rot 𝒗 = 1 𝜂 rot(𝛷𝛻𝑈) (3)
This is convenient to introduce the standard stream function Ψ :

𝑣 𝑥 = 𝜕Ψ 𝜕𝑧 , 𝑣 𝑧 = - 𝜕Ψ 𝜕𝑥 (4) By using (3), (4) in (2)
, we come to the following relation:

Δ 2 Ψ = - 1 𝜂 [ 𝜕𝑈 𝜕𝑧 𝜕Φ 𝜕𝑥 - 𝜕𝑈 𝜕𝑥 𝜕Φ 𝜕𝑧 ] (5)
The boundary conditions for the function Ψ are:

Ψ = 0, 𝑎𝑡 𝑧 = 0, 𝐿 𝑧 (6a)
∂Ψ 𝜕𝑧 = 0, 𝑎𝑡 𝑧 = 0, 𝐿 𝑧 (6b)

Ψ → 0, 𝑎𝑡 𝑥 → ±∞ (6c)
We should solve eq. ( 5), ( 6) together with the third equation of (1). In can be done only numerically.

To get some transparent analytical results, let us consider the initial stage of the oscillations, a few moments after the injection of a ferrofluid droplet, when the concentration Φ is practically equal to its initial profile Φ 0 (𝑥).

Since Lz<<Lx, all derivates of Ψ over z are much larger than the derivates over x. Taking it into account, the equation (5) reduces to:

𝜕 4 𝜕𝑧 4 Ψ = - 1 𝜂 𝜕𝑈 𝜕𝑧 𝜕Φ 0 (𝑥) 𝜕𝑥 (7)
Solution of the problem (4), ( 7) is:

Ψ = - 1 𝜂 𝜕Φ 0 (𝑥) 𝜕𝑥 ∫ 𝑑𝑧′ 𝑧 0 ∫ 𝑑𝑧′′ 𝑧′ 0 ∫ 𝑈(𝑥, 𝑧′′′, 𝑡)𝑑𝑧′′′ 𝑧′′ 0 + 1 6 𝐴𝑧 3 + 1 2 𝐵𝑧 2 (8)
where the integration constants A,B are determined from the boundary condition Ψ = 0, 𝜕 𝜕𝑧 Ψ = 0 at 𝑧 = 𝐿 𝑧 . We omit their explicit expressions for brevity. Thus the calculation steps can be summarized as follows. Using a known concentration profile Φ 0 (𝑥), we first calculate the stream function with the help of Eq. ( 8). Then, replacing Eq. ( 8) in Eq. ( 4), we obtain both components 𝑣 𝑥 and 𝑣 𝑧 of the velocity profile. The numerical results allow us to compare the obtained velocity values with those required for successful treatment of brain strokes by the considered method.

Results and discussion

The calculations of the velocity profile in the flat channel were done for a gaussian-like profile of initial concentration, rather typical for fluid injection processes:

Φ 0 (𝑥) = Φ 0 exp(- 𝑥 2 𝜎 2 ) , Φ 0 = 0.05; 𝜎 = 1с𝑚 (9)
This means that we consider a cloud of ferroparticles with the typical length 1cm. Another important parameter intervening into calculations is the magnetic susceptibility 𝜒 𝑎 of these particles. We take the value 𝜒 𝑎 = 25, close to the one (𝜒 𝑎 = 24) reported in experiments with ferrofluid microdroplets [START_REF] Bacri | Behavior of a magnetic fluid microdrop in a rotating magnetic field[END_REF]. Some results of calculations of the longitudinal component of the velocity 𝑣 𝑥 = 𝜕Ψ 𝜕𝑧 in middle plane of the channel (z=Lz/2) are shown in Fig. 3. The results show that for the chosen parameters of the system the velocity, close to the drop edges, can achieve several centimeters per second . The characteristic value (10 -3 -10 -2 m/s) obtained in our simulations is three-to-four orders of magnitude larger than the value 𝑣 𝑥 ≈ 3 × 10 -6 𝑚/𝑠 reported in experiments of [START_REF] Chang | Magnetic control of Fe3O4 nanomaterial for fat ablation in microchannel[END_REF] with fat ablation by generated microcirculation flows of magnetic suspension. Nevertheless, we warrant the reader that theoretical velocity values obtained in the present paper are valid only at the very onset after the magnetic field application, and the temporal evolution of the speed caused by spreading of the ferrofluid cloud along the channel has to be checked in order to confirm the effectiveness of the proposed method during a few minutesa typical timescale for the medical intervention previewed for the brain stroke treatment by this method [START_REF] Creighton | Magnetic-based systems for treating occluded vessels[END_REF].

Conclusion.

Our results show that the running oscillating field with a quite realistic amplitude can induce in a ferrofluid drop, therefore in a liquid, the drop is situated in, circulating flows with the velocity amplitude about several centimeters per second. This effect can be used for the intensification of drug delivery in blood vessels.

Fig. 1 .

 1 Fig.1. Sketch of the flat channel with a ferrofluid drop (gray rectangle) injected in the middle of the channel.

Fig. 2 .

 2 Fig.2. Absolute value H of the instantaneous local field H vs. the horizontal coordinate x for z=Lz/2 (middle of the channel thickness) for the time moments: t0=kT/2; t1=kT/2+T/8; t2=kT/2+T/4; k=0, 1, 2, …. The solenoid parameters are as follows: radius 10cm; length 5cm; the product I0N= 18 kA, where I0 is current amplitude and N is the number of turns; the current angular frequency 𝜔 = 15𝑠 -1 . The distance between the solenoids 2Lx=20cm; the channel thickness Lz=5mm

Fig. 3 .

 3 Fig.3. The longitudinal velocity 𝑣 𝑥 at the middle plane of the flat channel vs. coordinate x for several moments of time t0=kT/2; t1=kT/2+T/8; t2=kT/2+T/4; k=0, 1, 2, …. The other parameters of the system are described in the figure caption of Fig.2.

Acknowledgements

PK acknowledges the French 'Agence Nationale de la Recherche', Project Future Investments UCA JEDI, No. ANR-15-IDEX-01 (projects ImmunoMag and MagFilter) and the private company Axlepios Biomedicals for financial support,. AZ thanks the Russian Science Foundation, project 20-12-00031, for the financial support.