---
author:
- |
    Yunjiao Lu$^{1,2,\$}$, Pierre Hodara$^{1,\$}$, Charles Kervrann$^{2}$, Alain Trubuil$^{1}$\*\
    $^{1}$INRA, UR 1404, MaIAGE, Université Paris-Saclay, Jouy-en-Josas, France\
    $^{2}$Inria, Centre de Recherche Bretagne-Atlantique, EPC SERPICO, Rennes, France\
    ^\$^these authors contributed equally to this work\
    ^\*^corresponding author
bibliography:
- 'biblio.bib'
title: An inverse problem approach to the probabilistic reconstruction of particle tracks on a censored and closed surface
---

abstract {#abstract .unnumbered}
========

Investigation of dynamics processes in cell biology very often relies on the observation of sampled regions without considering re-entrance events.  In the case of image-based observations of bacteria cell wall processes, a large amount of the cylinder-shaped wall is not observed.  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 connection of particle trajectories on a cylindrical surface.  A subregion of the cylinder is typically observed during the observation period, and biomolecules may appear or disappear in any place of the 3D surface.  The performance of the method is mainly demonstrated on simulation data that mimic MreB dynamics observed in 2D time-lapse fluorescence microscopy. 

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, in the neighborhood of membranes.  Nevertheless, it is quite uncommon to analyze 3D dynamics occurring on a closed surface and observed on a 2D plane.  To our knowledge, identifying re-entrance events of the same entities inside the ROO is not addressed in the literature.  Yet, as soon as the unobserved region represents a significant part of the entire surface, a complete description of the dynamics on these closed surfaces is of paramount importance to decipher the mechanisms involved in cell wall construction of bacteria [@billaudeau2017contrasting].  Our objective is to provide a generic approach to tackle the problem of the reconstruction of particle tracks observed on a small part of a closed surface as illustrated on Fig.  \[fig:illus\_partialobservation\].  As inputs, we consider a set of trajectories estimated by tracking algorithms (e.g. [@jaqaman2008robust], [@chenouard2013multiple]).  These algorithms are very sophisticated and allows us to handle large sets of particles, different stochastic dynamical models [@blom1988interacting], [@bressloff2013stochastic], and observation models [@genovesio2006multiple], [@li2001target].  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 large hidden or blind region inside the region of interest, or to the identification of particles leaving the ROO through one border of the domain and re-entering from a far border. 

In this paper, we focus on the design and evaluation of a self-contained mathematical framework to tackle the reconstruction of particle tracks on cylindrical surfaces, given the observations inside a restricted window sampled on the surface.  3D image acquisition is not always appropriate, especially if the objective is to capture fast and temporally short events as described in [@billaudeau2017contrasting].  The frame rate adapted to the scale of dynamics may be too high when compared to the period of time to acquire one 3D volume.  In our study, the particles obey a stochastic Brownian motion with drift and may appear or disappear during the observation period.  In our modeling framework, we do not consider split or merge events.  The track reconstruction problem is defined as the maximization of the tracks likelihood given tracks segments reconstructed inside the region of observation.  The optimization problem to be solved is formulated as a integer linear programming problem.  We derive a data-driven algorithm with no hidden parameter to 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 speed and variance of particles as well as the birth and death rates. 

The remainder of the paper is organized as follows.  In Section 2, we present the problem more formally.  In Section 3, we describe a probabilistic framework, including Poisson processes used to describe birth and death events, and Brownian motion with drift to represent particle motion.  This mathematical framework allows us to derive a computational procedure aiming at connecting segments belonging to the same trajectory, and then recovering particles dynamics on the whole closed surface.  Note that we ignore the curvature of the cylinder and suppose that the movements are observed on a 2D unwrapped surface.  In Section 4, the performance of the segments connection algorithm is evaluated on simulated data.  Finally we conclude and propose some future work. 

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.  \[fig:illus\_partialobservation\]).  We associate 2D coordinates $(x,y) \in [-L,0] \times [0,H]$ to each point of the underlying cylindric manifold.  The particles are “born” with a constant rate $\lambda$ and appear uniformly at random on the membrane surface.  We consider a Poisson process with intensity $\lambda$ to statistically represent the birth events.  Each particle is assumed to have the same constant rate of death $\tau_d$ such that life duration $T_d$ of a particle follows an exponential law of parameter $\tau_d$. During its lifetime, a particle $k$ born at time $t_0$ and located at $\mathbf{Z}_0^k=(X_0^k,Y_0^k)$, moves according to Brownian motion with drift.  On the set $]-L,0[  \times [0,H],$ the position of the particle at time $t\geq t_0$ prior to its death time is given by $$\mathbf{Z}_t^k = \mathbf{Z}_0^k + \mathbf{v} (t-t_0) + \mathbf{\Sigma} \mathbf{B}_{t-t_0}^k
    \label{eq:bmd}$$ where $\mathbf{Z}_t^k=(X_t^k,Y_t^k)$, $\mathbf{v}=(v_x,v_y)$, $\mathbf{\Sigma} =
  \left[ {\begin{array}{cc}
   \sigma_x & 0 \\
   0 & \sigma_y \\
  \end{array} } \right]$, $\mathbf{B}_t^k$ is a two-dimensional Wiener process.  In order to model the topology of the cylinder as illustrated in Fig.  \[fig:illus\_partialobservation\], we impose deterministic jumps when the process reaches one of the two borders $\{-L\} \times [0,H]$ or $\{0\} \times [0,H]$. For any $y \in [0,H],$ the process reaching position $(-L,y)$ jumps to position $(0,y)$ and vice versa.  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 $\Delta t, 2\Delta t, 3\Delta t \dots$ We denote $\Delta t$ the time step on the subset $[-l,0] \times [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 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 segment of a given trajectory is an output if the last observed point of the segment is within a neighborhood of $\{0\} \times [0,H]$. Meanwhile, we consider that it is an input if the first observed point is within a neighborhood of $\{-l\}\times [0,H].$ Our main objective is then to associate the set of trajectories exiting the observed set $[-l,0] \times [0,H]$ with the set of trajectories entering this observation set.  The challenge is to correctly match the outputs and the inputs associated to particles (see Fig.  \[fig:illus\_partialobservation\]). 

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) \in O_S$ is characterized by its output time $t_o$ and its position $y_o \in [0,H]$ where the particle left the observed region.  Similarly each input $i=(t_i,y_i) \in I_S$ is characterized by its input time $t_i$ and its position $y_i \in [0,H]$ where it entered the observed region.  A particle “involved” in an output $o \in O_S$ either died after time $t_o$ in the unobserved region, or is “involved” in a given input $i \in I_S$ with $t_i>t_o.$ We will denote this event by $\{o \to i\}.$ Similarly a particle “involved” in an input $i \in I_S$ was either born before time $t_i$ in the unobserved region, or is “involved” in a given output $o \in O_S$ with $t_i>t_o$, which corresponds to the event $\{o \to i\}.$

Define $c=(D_c,B_c,b_c)$ with $D_c \subset O_S, \, B_c \subset I_S$ and $b_c$ a bijection from $O_S \setminus D_c$ to $I_S \setminus 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 $$\bigcap_{o \in O_s \setminus D_c} \{o \to 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.  \[fig:procedure\], to facilitate the understanding of the modeling steps. 

\[node distance = 1cm, auto,font=, every node/.style=[node distance=3cm]{}, comment/.style=[rectangle, inner sep= 5pt, text width=5cm, node distance=0.25cm, font=]{}, force/.style=[rectangle, draw, fill=black!10, inner sep=5pt, text width=5cm, text badly centered, minimum height=1.2cm, font=****]{}\]

(datas) [Sample $S$, the outputs and inputs sets $O_S$ and $I_S$]{}; (estim) [Estimation of parameters $\hat{\tau}_\alpha, \hat{\tau}_d, \hat{v}$ and $\hat{\sigma}$, (Section 3.3)]{}; (Proba) [Computation of the probability $P(\mathcal{B}_c)$ that the inputs in $B_c$ are born in the unobserved region.\
(section 3.1 Eqs. ( \[eq:BcEps\]) and (\[eq:ProbaBcEps\]))]{}; (Probd) [Computation of the probability $P(D_c)$ that the outputs in $D_c$ died in the unobserved region.\
(section 3.1 Eq. ( \[eq:death\]))]{}; (Probc) [Computation of the probability $P(\flat_c)$ that a set of outputs are connected to a set of inputs.\
(section 3.1 Eqs. ( \[eq:ProbabcEps\]) and (\[eq:probaOtoI\])]{}; (MaxiLikeli) [Likelihood\
(section 3.1 Eqs. ( \[eq:decomp1\]), (\[eq=probaCepsfinal\])-(\[eq:hatc\]))]{}; (cplex) [Maximization of likelihood with CPLEX, (section 3.2 Eqs. ( \[eq:DecompArgMax\]) and(\[eq:optim\]))]{}; (criteriaAri) [Tracks reconstruction and estimation of connection accuracy, (section 4.2)]{};

(comment-Proba) [(+) Lifetime of a particle $T_d \sim \mathcal{E} (\tau_d)$\
(+) The first passage time on $l$ of a Brownian motion $T_l$ follows an Inverse Gaussian distribution (Prop 1).]{}; ( comment-estim) [(+) To emphasize, $\hat{\tau}_\alpha$ is a novel estimator Eq.(19)]{}; (comment-cplex) [(+) The possibility to find the n^th^ optimal solution, second part of section (3.2)]{};

(datas) edge (estim) (comment-Proba) edge (Proba) (comment-Proba) edge (Probd) (comment-Proba) edge (Probc) (Proba) edge (MaxiLikeli) (Probd) edge (MaxiLikeli) (Probc) edge (MaxiLikeli) (MaxiLikeli) edge (cplex) (cplex) edge (criteriaAri);

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 $\hat c$ such that $P(\hat c / S)$ is maximal.  It is difficult to calculate directly $P(\hat c / S)$. Since $c \subset S \subset O_S,$ we can compute $P(\hat c / 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(\hat c / S)$ with the classical conditional formula $$P\left( c/S \right)=\frac{P\left( c / O_S\right)}{P(S/O_S)},$$ because it gives $P(S/O_S)=\sum_{c \in \mathcal{C}_S} P(c / O_S) =0.$

Therefore, for each input $i=(t_i,y_i) \in I_S$, we consider a spatiotemporal neighborhood $V_i^\epsilon = T_i^\epsilon \times H_i^\epsilon$ with $T_i^\epsilon= [t_i- \frac{\epsilon}{2}, t_i+ \frac{\epsilon}{2}]$ and $H_i^\epsilon= [y_i- \frac{\epsilon}{2}, y_i+ \frac{\epsilon}{2}]$ for some $\epsilon > 0.$

The idea is to replace a given configuration $c$ by a set $\mathcal{C}^\epsilon_c$ of configurations where each element $c^* \in \mathcal{C}^\epsilon_c$ is similar to $c$ but each input $i \in I_S$ is replaced by an input in $V_i^\epsilon.$ Formally, for each configuration $c$ leading to the input set $I_S,$ $\mathcal{C}^\epsilon_c$ is the set of configurations defined as follows: $c^*=(D_c^*,B_c^*,b_c^*) \in \mathcal{C}^\epsilon_c$ if and only if for each $i \in I_S,$ there exist $i_\epsilon^* \in V_i^\epsilon$ satisfying: $$\begin{aligned}
\left\{
\begin{array}{l}
D_c^*=D_c,\\
B_c^*= \{i_\epsilon^*, i \in B_c \},\\
\mbox{For each}~ i \in I_S \setminus B_c, \, b_{c^*}\left( b_c^{-1} (i) \right) = i_\epsilon^*. 
\end{array}
\right.\end{aligned}$$ With this definition, we have

$$\label{eq:LimEps}
P(c/S)= \lim_{\epsilon \to 0} P \left( \mathcal{C}^\epsilon_c / S \right) = \lim_{\epsilon \to 0} \frac{P \left( \mathcal{C}^\epsilon_c / O_S \right)}{\sum_{c'\in \mathcal{C}_S} P \left( \mathcal{C}^\epsilon_{c'} / O_S \right)}.$$

In what follows, we study the behavior of $P(\mathcal{C}^\epsilon_c / O_S)$ when $\epsilon$ 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(\mathcal{C}^\epsilon_c)$ instead of $P(\mathcal{C}^\epsilon_c / O_S)$ in order to simplify the notations.  The study of $P(\mathcal{C}^\epsilon_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^\epsilon.$

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

Consider an output $o \in O_S$ and the possibility for the particle involved in $o$ to die in the unobserved region.  We have the following proposition (see [@schrodinger1915], [@tweedie1945inverse], [@wald1973sequential]). 

\[prop:IG\] Given the particle motion model as Brownian motion with drift as described in equation \[eq:bmd\], the first passage time noted as $T_l$ on the entrance line $\{-l\} \times [0,H]$ of a particle starting at position $z_0=(0,y_0)$ for some $y_0 \in [0,H],$ follows a law of inverse Gaussian, that is $T_l \sim IG \left( \frac{l_u}{\mathbf{v}_x}, \big(\frac{l_u}{\sigma_x}\big)^2 \right)$ where $l_u:=L-l$ is the length of the unobserved region. 

Recall that if $X\sim IG(\mu,\lambda)$, then $X \geq 0$ almost surely, and for each $x\geq 0,$ $$P(X\leq x)=\int_0^x\sqrt{\frac{\lambda}{2\pi y^3}}\exp \big(-\frac{\lambda(y-\mu)^2}{2\mu^2y}\big) dy.$$

In our framework, the event corresponding to the death of a particle with life duration $T_d$ following an exponential law of parameter $\tau_d$ in the unobserved region is precisely $\{T_d < T_l\}$. Hence, We can derive an explicit expression of $P(\mathcal{C}^\epsilon_c).$

Assume $\epsilon$ small enough so that for each $i,i' \in I_S, \, T_i^\epsilon \cap T_{i'}^\epsilon = \emptyset.$ For a given configuration $c$ and a given $\epsilon >0,$ we will write $\mathcal{C}^\epsilon_c = (D_c, \mathcal{B}_c^\epsilon, \flat_c^\epsilon)$ with $\mathcal{B}_c^\epsilon = \{B_c^*, c^* \in \mathcal{C}^\epsilon_c \}$ and $\flat_c^\epsilon = \{b_{c^*}, c^* \in \mathcal{C}^\epsilon_c \}.$

Due to the independent behavior of the particles, we have the following decomposition: $$\label{eq:decomp1}
P(\mathcal{C}^\epsilon_c)=P(D_c) P(\mathcal{B}_c^\epsilon) P(\flat_c^\epsilon).$$ We can then compute separately the probabilities of events $D_c, \, \mathcal{B}_c^\epsilon$ and $\flat_c^\epsilon$. First, note that we can assume without loss of generality that each output $o \in D_c$ starts at time $t_o=0$ and that only the position $y_o \in [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 us that the life duration $T_d$ of the particle after the the output $o$ still follows an exponential law of parameter $\tau_d.$

Since all outputs behave identically and independently, we have $P(D_c) = P( T_d < T_l)^{|D_c|},$where $|D_c|$ stands for the cardinal of $D_c.$ According to proposition \[prop:IG\], and since $T_d$ and $T_l$ are independent, we have $$\begin{aligned}
\label{eq:death}
P( T_d < T_l) &=& \int_0^{+\infty} \int_0^{t_l} f_{T_d}(t_d) f_{T_l}(t_l) dt_d \; dt_l, \\
&=& \int_0^{+\infty} \int_0^{t_l} \tau_d e^{-\tau_d t_d} \frac{l_u}{\sigma_x\sqrt{2\pi t_l^3}}\exp \left(-\frac{(\mathbf{v}_xt_l-l_u)^2}{2\sigma_x^2t_l}\right) dt_d \;dt_l, \nonumber \\
&=& \int_0^{+\infty}  \frac{l_u \left( 1 - e^{-\tau_d t_l} \right)}{\sigma_x\sqrt{2\pi t_l^3}}\exp \left(-\frac{(\mathbf{v}_xt_l-l_u)^2}{2\sigma_x^2t_l}\right) dt_l, \nonumber \end{aligned}$$ where $f_{T_d}$ and $f_{T_l}$ respectively stand for the density functions of $T_d$ and $T_l.$

Now, consider the event $\mathcal{B}_c^\epsilon.$ We call “spontaneous input” an input related to a particle born in the unobserved region that has never been observed.  The set $\mathcal{B}_c^\epsilon$ is defined so that, for each input $i \in B_c$, we have exactly one “spontaneous input” appearing during the time interval $T_i^\epsilon,$ with a position in $H_i^\epsilon.$ Moreover, outside $\cup_{i \in B_c} T_i^\epsilon,$ there is no “spontaneous input”.  Formally, we have $$\label{eq:BcEps}
\mathcal{B}_c^\epsilon = \left\{ N_\alpha \left( [0,T_S] \setminus \bigcup_{i \in B_c} T_i^\epsilon \right) =0 \right\} \bigcap \left( \bigcap_{i \in B_c} \Big( \{N_\alpha (T_i^\epsilon)=1\} \cap H_i^\epsilon \Big) \right),$$ where $N_\alpha$ is a Poisson process of intensity $\tau_\alpha$ 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^\epsilon$ denotes also the event of “spontaneous” appearance of an input $i$ in $H_i^\epsilon$. This event is independent of the process $N_\alpha$, and since the “spontaneous inputs” appear uniformly on $[0,H],$ we have $P(H_i^\epsilon)=\frac{\epsilon}{H}.$

Meanwhile, for any time interval $I, \, N_\alpha(I)$ follows a Poisson law of parameter $\tau_\alpha |I|$ where $|I|$ denotes the length of the interval $I.$ Since $\epsilon$ is small enough so that for each $i,i' \in I_S, \, T_i^\epsilon \cap T_{i'}^\epsilon = \emptyset, \, N_\alpha( T_i^\epsilon)$ and $N_\alpha( T_{i'}^\epsilon)$ are independent.  Consequently, we can compute $P(\mathcal{B}_c^\epsilon)$ as follows: $$\begin{aligned}
\label{eq:ProbaBcEps}
P(\mathcal{B}_c^\epsilon)  = e^{-\tau_\alpha (T_S - |B_c|\epsilon)} \left( \epsilon \tau_\alpha e^{-\epsilon \tau_\alpha} \frac{\epsilon}{H} \right)^{|B_c|} 
& = \left( \frac{\epsilon^2 \tau_\alpha}{H} \right)^{|B_c|} e^{-\tau_\alpha T_S}.\end{aligned}$$

Finally, consider the event $\flat_c^\epsilon$. For each input $i \in I_S \setminus B_c,$ we denote by $\{ o_c^i \to V_i^\epsilon \}$ the survival event of the particle involved in the output $o_c^i=b_c^{-1}(i)$ in the unobserved region which appears in the spatiotemporal neighborhood $V_i^\epsilon$. Since the particles behave independently, we have $$\label{eq:ProbabcEps}
P(\flat_c^\epsilon)=\prod_{i \in I_S \setminus B_c} P \Big(\{ o_c^i \to V_i^\epsilon \} \Big).$$

In the sequel, we consider a given input $i \in I_S \setminus B_c$ and its related output $o=b_c^{-1}(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 $\tau_d.$ During its lifetime, the position of the particle is driven by a brownian motion with drift $\mathbf{Z}_t=(X_t,Y_t)$: $\mathbf{Z}_t = \mathbf{v} t + \mathbf{\Sigma} \mathbf{B}_t,$ where $\mathbf{B}_t$ is a two-dimensional Wiener process and $\mathbf{v}$ and $\mathbf{\Sigma}$ are given in Equation (\[eq:bmd\]).  Define $T_{l}$ the first reaching time of $l_u=L-l$ of the process $X_t$. The event $\{ o \to V_i^\epsilon \}$ can now be written as follows: $$\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\{ o_c^i \to V_i^\epsilon \}= \left\{T_d > T_l \right\} \bigcap \left\{T_l \in \left[s_i- \frac{\epsilon}{2},s_i + \frac{\epsilon}{2}\right] \right\} \bigcap \left\{Y_{T_l} \in \left[h_i- \frac{\epsilon}{2},h_i + \frac{\epsilon}{2}\right] \right\}.$$ This expression corresponds exactly to the fact that in order to realize $\{ o_c^i \to V_i^\epsilon \}$ the particle needs to have a life duration longer than its first reaching time of $l_u$ and to appear in the spatiotemporal neighborhood $ \left[s_i- \frac{\epsilon}{2},s_i + \frac{\epsilon}{2}\right] \times \left[h_i- \frac{\epsilon}{2},h_i + \frac{\epsilon}{2}\right].$ Furthermore, $T_d$ follows an exponential law of parameter $\tau_d,$ $Y_t$ follows a Gaussian law of parameters $\mathbf{v}_yt$ and $\sigma_y^2t$ and $T_l \sim IG \left( \frac{l_u}{\mathbf{v}_x}, \big(\frac{l_u}{\sigma_x}\big)^2 \right)$. Moreover, due to the fact that $\mathbf{\Sigma}$ is diagonal, the process $Y_t$ is not only independent of $T_d$ but also of $T_l.$ This allows us to write $$P \left( \{ o_c^i \to V_i^\epsilon \} \right)= \int_{s_i-\frac{\epsilon}{2}}^{s_i + \frac{\epsilon}{2}} f_{T_l}(t_l) \left( \int_{t_l}^{+ \infty} f_{T_d}(t_d) \left( \int_{h_i-\frac{\epsilon}{2}}^{h_i + \frac{\epsilon}{2}} f_{Y_{t_l}}(y) dy \right) dt_d \right) dt_l.$$ As the two integrals involve a small domain of size $\epsilon, \, P \left( \{ o_c^i \to V_i^\epsilon \} \right) \sim \mathcal{O} (\epsilon^2)$, and $$\begin{aligned}
\label{eq:probaOtoI}
\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!  \lim_{\epsilon \to 0} \frac{P \left( \{ o_c^i \to V_i^\epsilon \} \right)}{\epsilon^2} = f_{T_l}(s_i) f_{Y_{s_i}}(h_i) \int_{s_i}^{+ \infty} f_{T_d}(u) du  \\
=  \frac{l_u}{\sigma_x\sqrt{2\pi s_i^3}}\exp \left(-\frac{(\mathbf{v}_xs_i-l_u)^2}{2\sigma_x^2s_i}\right) \frac{1}{\sigma_y \sqrt{ 2 \pi s_i}} \exp \left( - \frac{\left(h_i -  \mathbf{v}_y s_i\right)^2 }{2 \sigma_y^2 s_i} \right) e^{-\tau_d s_i} \nonumber \\
= \frac{l_u}{2 \pi \sigma_x \sigma_y s_i^2} \exp \left( -\frac{(\mathbf{v}_xs_i-l_u)^2}{2\sigma_x^2s_i} - \frac{\left(h_i -  \mathbf{v}_y s_i\right)^2 }{2 \sigma_y^2 s_i} -\tau_d s_i \right).  \nonumber\end{aligned}$$

For each configuration $c$, we derive the likelihood $Q(c)$ of the configuration $c$ as follows: $$Q(c)~:= ~\lim_{\epsilon \to 0} \frac{P(\mathcal{C}^\epsilon_c)}{\epsilon^{2|I_S|}}.$$ From (\[eq:decomp1\]) and Equations (\[eq:death\], \[eq:ProbaBcEps\], \[eq:ProbabcEps\] and \[eq:probaOtoI\]), we finally obtain the likelihood $$\begin{aligned}
\label{eq=probaCepsfinal}
\! \! \! \! \! \!\! \! \! \! \! \! \! \! \!\! \! \!\! \! \! \!   Q(c) ~=~\left( \frac{\tau_\alpha}{H} \right)^{|B_c|} e^{-\tau_\alpha T_S} \left(\int_0^{+\infty}  \frac{l_u \left( 1 - e^{-\tau_d t_l} \right)}{\sigma_x\sqrt{2\pi t_l^3}}\exp \left(-\frac{(\mathbf{v}_xt_l-l_u)^2}{2\sigma_x^2t_l}\right) dt_l \right)^{|D_c|} \nonumber \\
~~~~\times \prod_{i \in I_S \setminus B_c} \left[ \frac{l_u}{2 \pi \sigma_x \sigma_y s_i^2} \exp \left( -\frac{(\mathbf{v}_xs_i-l_u)^2}{2\sigma_x^2s_i} - \frac{\left(h_i -  \mathbf{v}_y s_i\right)^2 }{2 \sigma_y^2 s_i} -\tau_d s_i \right) \right].\end{aligned}$$

Note that the limit when $\epsilon$ goes to $0$ of $\frac{P(\mathcal{C}^\epsilon_c)}{\epsilon^{2|I_S|}}$ is well defined, strictly positive, and that the exponent $2 |I_S|$ does not depend on the configuration $c.$

Recall (\[eq:LimEps\]), this allows us to write

$$\label{eq:LimEpsQ}
P(c/S)=  \frac{Q(c)}{\sum_{c'\in \mathcal{C}_S} Q(c')}$$

and as a consequence, we have

$$\label{eq:hatc}
\hat c= \arg\!\max_{c \in \mathcal{C}_S} \{ Q(c) \}.$$

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 (\[eq=probaCepsfinal\])).  Define $$\beta:= - \log \left(\frac{\tau_\alpha}{H}\right),$$ $$\delta := - \log \left(\int_0^{+\infty}  \frac{l_u \left( 1 - e^{-\tau_d t_l} \right)}{\sigma_x\sqrt{2\pi t_l^3}}\exp \left(-\frac{(\mathbf{v}_xt_l-l_u)^2}{2\sigma_x^2t_l}\right) dt_l \right)$$ and for each configuration $c$ and each $i \in I_S \setminus B_c$ $$\gamma_c^i:= - \log \left[ \frac{l_u}{2 \pi \sigma_x \sigma_y s_i^2} \exp \left( -\frac{(\mathbf{v}_xs_i-l_u)^2}{2\sigma_x^2s_i} - \frac{\left(h_i -  \mathbf{v}_y s_i\right)^2 }{2 \sigma_y^2 s_i} -\tau_d s_i \right) \right].$$ It follows that $$\begin{aligned}
\label{eq:DecompArgMax}
\hat c ~=~  \arg\!\max_{c\in C} Q(c) &=& \arg\!\min_{c\in C} - \log \left(Q(c) \right) \\\nonumber
 &=& \arg\!\min_{c\in C} \left( \beta |B_c| + \delta |D_c| + \sum_{i \in I_S \setminus B_c} \gamma_c^i \right).\end{aligned}$$ This decomposition allows us to consider a linear optimization problem where $\beta$ represents the cost of the spontaneous birth of an input, $\delta$ the cost of the death of an output and $\gamma_c^i$ the cost of the connection between the output $b_c^{-1}(i)$ and the input $i.$ The cost of connection can be defined for any couple $(o,i) \in O_S \times I_S$ as $$\gamma_o^i:= - \log \left[ \frac{l_u}{2 \pi \sigma_x \sigma_y s_{o , i}^2} \exp \left( -\frac{(\mathbf{v}_xs_{o , i}-l_u)^2}{2\sigma_x^2s_{o , i}} - \frac{\left(h_{o , i} -  \mathbf{v}_y s_{o , i}\right)^2 }{2 \sigma_y^2 s_{o , i}} -\tau_d s_{o , i} \right) \right],$$ where $s_{o,i}:= t_i-t_o, \, h_{o,i}=y_i-y_o$ and the convention $\gamma_o^i = + \infty$ if $t_i \leq 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) \in O_S \times I_S}$ such that $c^{o,i}=1$ if $b_c(o)=i$ and $c^{o,i}=0$ if $b_c(o) \neq i.$ Since an output can be connected to at most one input, for each $o \in O_S$, $\sum_{i \in I_S}c^{o,i} \in \{0,1\}$ and $\sum_{i \in I_S}c^{o,i} = 0$ corresponds to the death of the output $o.$ Similarly, for each $i \in I_S,$ $\sum_{o \in O_S}c^{o,i} \in \{0,1\}$ and $\sum_{o \in 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) \in O_S \times I_S}$ that minimizes the quantity $$\beta \left(\sum_{i \in I_S} \left(1- \sum_{o \in O_S} c^{o,i} \right) \right) + \delta \left( \sum_{o \in O_S} \left(1- \sum_{i \in I_S} c^{o,i} \right) \right) + \sum_{o \in O_S} \sum_{i \in I_S} \gamma_o^i c^{o,i}$$ or equivalently $$\begin{aligned}
\!\! \! \! \! \! \! \! \! \! \!\!\! \! \! \! \! \! \! \! \!  \!K(c):=\sum_{o \in O_S} \sum_{i \in I_S} \left(\gamma_o^i -\beta -\delta \right) c^{o,i}   ~~\mbox{subject to}~
\left\{
\begin{array}{l}
\forall o \in O_S, \, \forall i \in I_S, \, c^{o,i} \in \{0,1\}, \\ 
\forall o \in O_S, \, \sum_{i \in I_S} c^{o,i} \in \{0,1\},\\ 
\forall i \in I_S, \, \sum_{o \in O_S} c^{o,i} \in \{0,1\}.  
\end{array}
\right. 
\label{eq:optim}\end{aligned}$$ In order to avoid to have infinite costs $\gamma_o^i$ when $t_i \leq t_o,$, we can also impose $c^{o,i}=0$ if $t_i \leq t_o$. Actually the problem (\[eq:optim\]) is a conventional linear optimization problem which can be solved by applying the CPLEX Linear Programming solver (https://www.ibm.com/analytics/cplex-optimizer). 

The configuration $\hat c$ is then the solution of the optimization problem (\[eq:optim\]) 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 by solving (\[eq:optim\]) with additional constraints ensuring that the solution is different from the previous ones.  In other words we define recursively the sequence $\left( c_n \right)_{n \in \mathbb N}$ in the following way:

-   $c_1 := \hat c$

-   $\forall n \geq 2, \, c_n$ solves (\[eq:optim\]) with the $n-1$ additional constraints $$\forall k  \in \{ 1, \dots, n-1 \}, \,  \sum_{o \in O_S} \sum_{i \in I_S} \left[ c_n^{o,i} (1-c_k^{o,i}) + (1-c_n^{o,i})c_k^{o,i} \right] \geq 1.$$

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 $\mathbf{v}$ and $\mathbf{\Sigma}$ can be estimated with classical maximum likelihood estimation procedures. 

Second, we propose an estimator $\hat \tau_d$ of $\tau_d$ as explained below.  The sample $S$ can be considered as a set of points $p=(t_p,\mathbf{Z}_p)$ observed at time $t_p$ and position $\mathbf{Z}_p=(X_p,Y_p)$ grouped in clusters $s$ corresponding to segments of trajectories.  The death of a particle in the observed region is detected in $S$ for each point $p \in S$ for which the associated segment $s_p$ has no successor point at time $t_p+\Delta 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 \subset S$ the sample of points in the restricted region.  For each point $p \in S_r,$ we denote by $D_p$ the event corresponding to the absence of successor for $p.$ This correspond to the fact that the particle involved in $p$ died during the time interval $[t_p, t_p + \Delta t].$ Since the life duration $T_d$ of a particle follows an exponential law of parameter $\tau_d$, and the absence of memory property of the exponential law, we have $$\label{eq:probDp}
P(D_p)=P(T_d \in [0,\Delta t])=1-e^{-\tau_d \Delta t}.$$ Hence, we define our estimator $\hat \tau_d$ as $$\label{eq:defHatTaud}
\hat \tau_d= \frac{1}{\Delta t |S_r|} \sum_{p \in S_r} 1[D_p],$$ where $|S_r|$ stands for the number of points in $S_r$ and $1[\cdot]$ 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 $+ \infty,$ the strong law of large numbers yields to $$\lim_{|S_r| \to \infty} \hat \tau_d = \frac{1-e^{-\tau_d \Delta t}}{\Delta t} \quad \quad a.s.$$ The justification of this choice for $\hat \tau_d$ relies in the following almost sure convergence: $$\label{eq:HatTauConsistent}
\lim_{\Delta t \to 0} \lim_{|S_r| \to \infty} \hat \tau_d = \tau_d \quad \quad a.s.$$ Our estimator $\hat \tau_d$ is then consistent as $\Delta 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_\alpha$ denotes the $\alpha$-quantile of the standard normal distribution, we have the following confidence interval of level $\alpha$ for $\frac{1-e^{-\tau_d \Delta t}}{\Delta t}$: $$CI_\alpha = \left[  \hat \tau_d - q_\alpha \sqrt{\frac{\hat \tau_d \left( \frac{1}{\Delta t} - \hat \tau_d \right)}{|S_r|}} ,  \hat \tau_d + q_\alpha \sqrt{\frac{\hat \tau_d \left( \frac{1}{\Delta t} - \hat \tau_d \right)}{|S_r|}} \right].$$ If $\Delta t$ is small enough, we get a good approximation of a confidence interval of level $\alpha$ for $\tau_d$ since $$\lim_{\Delta t \to 0} \frac{1-e^{-\tau_d \Delta t}}{\Delta t} = \tau_d.$$

Now, we describe the estimation procedure for the rate $\tau_\alpha$ of “spontaneous inputs” induced by particles born in the unobserved region $[-L,-l]\times [0,H]$ and reached the border $\{-l\}\times[0,H].$ We assume here that the parameters $\mathbf{v}, \, \mathbf{\Sigma}$ and $\tau_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 segments born in the region $]-x,0]\times[0,H]$ and reached the border $\{0\}\times[0,H].$ Accordingly, $\frac{N_{l_u}}{T_S}$ is a consistent estimator of $\tau_\alpha$ since the dynamics are assumed to be homogeneous on the surface of cylinder.  Our aim is actually to build an estimator for $\tau_\alpha$ in the case where $l_u >l$ which prevents us to compute directly $N_{l_u}$. Therefore, we compute $N_l$ by taking the whole observed region into account, and denote by $S_l^*$ the set of segments having an input in $\{-l\}\times[0,H]$ and an output in $\{0\}\times[0,H].$ For each segment $s \in S_l^*$ and each length $x \in [0,l_u],$ we denote by $B_s^x$ the event corresponding to the birth of $s$ within $]-l-x,-l]\times [0,H].$ Let $l_e=l_u-l$ be the length of the extended zone $[-l_u,-l] \times [0,H],$. We are now interested in the realization of the events $B_s^{l_e}.$

In Fig.  \[fig:illust\_Bs\], $N_l=2$ correspond to tracks \#1 and \#4, $S_l^*= \{2,5\},$ and the event $B_2^{l_e}$ is realized while $B_5^{l_e}$ is not. 

Note that since the particles have the same independent dynamics, $P\left( B_s^x \right)$ does not depend on $s.$ For $x < l,$ this probability can easily be estimated as follows: $$\hat p_x = \frac{N_x}{\big| S_o \big|},$$ where $S_o$ is the set of segments having an output in $\{0\}\times[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 $\hat \tau_\alpha$ as follows: $$\label{eq:defHatTauAlpha}
\hat \tau_\alpha = \frac{N_l + \hat p_{l_e} |S_l^*|}{T_S}.$$

![ An artificially constructed zone $]-l_u,0]\times [0,H]$ having the same size as the unobserved region $]-L,-l]\times [0,H].$ The observed region is $]-l,0]\times [0,H]$ ; as the width of the invisible part is $l_u$, the extended zone has width $l_e=l_u-l.$[]{data-label="fig:illust_Bs"}](images/EstimTauAlpha.png)

Intuitively, this estimator amounts to counting the number of particles reaching $\{0\} \times [0,H]$ with weight 1 for each segment that we actually saw being born in the observed region and with weight $\hat p_{l_e}$ for each spontaneous input that appeared in $\{-l\} \times [0,H].$ Note that, as $N_{l_u}=N_l + \sum_{s \in S_l^*} 1[B_s^{l_e}], \, \hat \tau_\alpha$ is an unbiased estimator of $\tau_\alpha.$ Moreover, if we assume that the number of observed segments grows linearly with the observation time $T_S,$ this estimator is consistent when $T_S$ goes to $+ \infty.$

Now, we consider the case $l< l_e <2l$ which can easily be extended to the general case $l< l_e.$ Consider $s \in S_l^*$ and denote for each interval $J \subset [-L,0]$ the event $B_s^J$ where the segment $s$ is born in the region $J \times [0,H].$ The event $B_s^{l_e}$ can be decomposed as follows: $$B_s^{l_e} = B_s^{[-2l,-l]} \bigcup \left( \overline{B_s^{[-2l,-l]}} \cap B_s^{[-l_u,-2l]} \right).$$ The loss of memory and homogeneity properties of the dynamics leads to the following estimator $\hat p_{l_e}$: $$\hat p_{l_e} := \hat p_l + \left( 1-\hat p_l \right) \hat p_{l_e-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 $\tau_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 segment $\{y\}\times ]-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. 

Simulation study
================

We present next a series of experiments performed on synthetic datasets.  These experiments aim at the evaluation and sensitivity analysis of the reconstruction procedure when varying the characteristics of the dynamics as well as the spatio-temporal support of observations.  In addition to demonstrate the potential of our procedure, these experiments might also be useful for the design of the experimental setting when considering the observation of censored dynamics on closed surfaces.  The reconstruction procedure has been implemented in Matlab. 

Generation of trajectories
--------------------------

Trajectories are generated on a rectangular unwrapped cylindrical surface of size $[0,L]\times[0,H]$ (Fig.  \[track\_gen\]).  In our experiments, we set $L = 50, H = 30$. The initial position of each trajectory is drawn uniformly on the surface.  Time duration $T$ between two births follows an exponential law with birth rate parameter $\lambda$. At each birth, the intrinsic properties of a trajectory $i$ are given, such as velocity $\mathbf{v}_i$, variance $\mathbf{\Sigma}_i$, and lifetime $T_d^i$. The lifetime $T_d$ follows an exponential law, with a constant death rate $\tau_d$ for all trajectories in the whole simulated image sequence.  The speed $\mathbf{v}_i = (v_{xi}, v_{yi})$ and noise $\mathbf{\Sigma}_i =
  \left[ {\begin{array}{cc}
   \sigma_x & 0 \\
   0 & \sigma_y \\
  \end{array} } \right]$ are constant along one given trajectory.  It is set $v_x = 0.01 v_y$ and $ \sigma_x = \sigma_y= \sigma$.

![ A set of simulated trajectories during 2.5 minutes (in stationary regime).  $\lambda = 0.03, \tau_d = 0.005, v_x = 0.6$ and $\sigma = 0.2$. 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.[]{data-label="track_gen"}](images/track_gen)

As there is no particle on the surface at the beginning, the simulated trajectories set needs some warm-up time to reach the stationary regime, where the law of the number of trajectories does not depend on time.  According to the supposed dynamic process, the expectation of the trajectories number $N$ during stationary regime is $E(N)=\frac{\lambda}{\tau_d}$. To ensure that the dynamics are in a stationary regime, the images sequence is simulated long enough, for around 2 hours (Fig.  \[nb\_tracks\]).  The time interval between two images is $\Delta t =0.25s$. Figure \[nb\_tracks\] shows how the number of trajectories on the surface changes along time $t$, when $\lambda = 0.03$ and $\tau_d = 0.005$. It can be noticed that after a while, around $t = 20$ min, the trajectories number fluctuates around the theoretical expectation value 6. 

![ Fluctuations of the number of trajectories w.r.t.  time.  The dynamics parameters are the same as in Fig.  \[track\_gen\]](images/Nb_tracks)

.  \[nb\_tracks\]

Analysis of the connection criterion: The Adjusted Rand Index
-------------------------------------------------------------

Given the true and estimated class assignments, we compute the Adjusted Rand Index to evaluate similarity or consensus between the two sets.  The Rand index ignores permutation and is symmetric: swapping the arguments does not change the score.  Perfect matching is scored 1.0 and 0 or negative on the contrary: 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 = \frac{a+b}{C_2^M},$$ where $C_2^M$ 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 [@rand1971objective]: $$ARI = \frac{RI-E(RI)}{1-E(RI)}.$$

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 connection procedure relies on the estimation of the characteristics of the dynamics: the speed, $\mathbf{v}$, the diffusion variance, $\mathbf{\Sigma}$, the arrival rate $\tau_{\alpha}$ and the death rate $\tau_d$, as these quantities are used in the calculation of the likelihood (Fig.  \[fig:procedure\]).  In this section, we evaluated the impact of several combinations of parameters of the dynamics $(\mathbf{v}, \mathbf{\Sigma}, \tau_{\alpha}, \tau_d)$ and spatio-temporal sampling, $(l/l_u, T_S)$ (Appendix , Table \[TabImpact\]) on the estimators and the accuracy of the reconstruction. 

### Performance of estimator $\hat \tau_\alpha$

In the following section, the proposed estimator $\hat{\tau}_\alpha$ defined in eq.(\[eq:defHatTauAlpha\]) is evaluated by simulation.  The impact of spatial (resp.  temporal) sampling through $l/l_u$ (resp.  $T_S$) as well as of birth rate $\lambda$ and death rate $\tau_d$ is tested. 

To calculate the “true” value of $\tau_\alpha$ in the invisible zone of size $l_u \times H, l_u = L-l$, we propose to count the arrival trajectories during $T_S$ and approximate $\tau_\alpha$ by $N_{l_u}/T_S$. This value can serve as ’ground truth’.  Then, by varying the value of $l<l_u, l_u$ is fixed, we test the robustness of the estimator $\hat{\tau}_\alpha$.

![ The normalized $\hat{\tau}_{\alpha}$ as $l/l_u$ fluctuates, $T_S = 30$ min.  The proportion of $l$ to $l_u$ is set from 0.42 to 1 (the lowest value $l/l_u=0.42$ corresponds to a realistic situation in 2D time-lapse fluorescence images).  The value of $y$ axis is the difference between the estimated value $\hat{\tau}_\alpha$ and the ’ground truth’, normalized w.r.t.  the ’ground truth’.  It can be noticed that, at $l/l_u$=0.42, more than a half of the trials presents a relative error smaller than $10\%$. The proposed estimator $\hat{\tau}_\alpha$ is unbiased and the variance decreases as $l/l_u$ increases.[]{data-label="fig:tau_alpha_portion"}](images/tau_alpha_2nd)

In Fig.  \[fig:tau\_alpha\_portion\], each box is computed from 1000 replications and the movie length used to calculate $\hat{\tau}_\alpha$ is 30 min, which is long enough to have a stable value.  However, in reality it is difficult to obtain a film as long as 30 min, because of photo bleaching and the natural growth in living samples.  Therefore, the robustness of $\hat{\tau}_\alpha$ w.r.t.  $T_S$ is also evaluated.  The proportion $l/l_u$ is set to 0.42, which corresponds to the real experimental data.  The result is shown in Fig.  \[fig:tau\_alpha\_time\].  $T_S$ varies from 2.5 min to 30 min.  It can be noticed that the counted ’ground truth’ (blue boxes) decreases as $T_S$ becomes longer.  Actually the counted value is only a pseudo ’ground truth’ and it is sensible to $T_S$ especially when $T_S$ is small.  The counted value converge to the hidden true value as $T_S \sim \infty$. The distributions of count ’ground truth’ and estimator become close to each other for $T_S \geq 10$ min. 

![ The counted ’ground truth’ $\tau_\alpha$ and the estimated $\hat{\tau}_{\alpha}$ obtained by images movies of different length, varying from 2.5 min to 30 min.  Blue boxes show the “ground truth” when all the information, even the invisible region, is available.  Magenta boxes correspond to the estimator $\hat{\tau}_\alpha$. The blue line is the counted value when $T_S=30$ min. []{ data-label="fig:tau_alpha_time"}](figs/tau_alpha_time)

The absolute value of $\tau_\alpha$ depends on $\lambda$ and $\tau_d$. Fig.  \[fig:hat\_tau\_alpha\_5min\_delta\_t0\_25\] displays for different combinations of $\lambda$ and $\tau_d$, the estimations of $\hat{\tau}_\alpha$ by 5-min movies (magenta) and 30-min movies (blue).  It shows that $\tau_\alpha$ increases linearly as the birth rate $\lambda$ increases, and decreases slightly linearly as the death rate $\tau_d$ increases. 

![ The estimation of arrival rate $\tau_\alpha$ w.r.t.  different $\lambda$ and $\tau_d$. For example, when $\lambda=0.04, \tau_d = 0.004,$ the median value of $\hat{\tau}_\alpha$ is around 0.025, which means that at each moment, the probability that a particle born in the invisible zone arrives at $\{-l\}\times[0,H]$ is around 0.025.[]{data-label="fig:hat_tau_alpha_5min_delta_t0_25"}](figs/tau_alpha_5min_delta_t0_25)

### Performance of the estimator $\hat{\tau}_d$

As explained in previous section, $\hat{\tau}_d$ is a rather classical estimator.  Fig.  \[fig:hat\_tau\_d\_5min\_delta\_t0\_25\] shows the estimator with 5-min movies (magenta) and 30-min movies (blue) respectively.  It confirms that the estimator is unbiased.  Black lines represent the true value of $\tau_d$. Naturally, the variance is bigger with shorter movies.  Some estimators of death rate with movies of 5 min falls to zero, which means that the particle is immortal.  This error will surely cause connection error because the algorithm will force the connection.  We found that $T_S = 5$ min is a good choice to limit the estimation error of $\tau_d$ and to ensure a good connection performance (more details in Appendix). 

![ $\hat{\tau}_d$ with different $\lambda$ and $\tau_d$. The estimator is unbiased.  The variance of the estimator is larger with shorter movies (magenta).  For a given $\tau_d$, when birth rate $\lambda$ increases, and therefore the number of particles also increases, then the variance decreases.[]{data-label="fig:hat_tau_d_5min_delta_t0_25"}](figs/tau_d_5min_delta_t0_25)

According to Figs.  \[fig:hat\_tau\_alpha\_5min\_delta\_t0\_25\] and \[fig:hat\_tau\_d\_5min\_delta\_t0\_25\], when $T_S = 30$ min, the estimators of $\tau_\alpha$ and $\tau_d$ perform well, being converged with small variance.  As 30-min movies acquisition is almost infeasible under the situation of fluorescence microscopy, these estimators are noted as $\tilde{\tau}_\alpha$ and $\tilde{\tau}_d$. In the following, the estimators obtained with tested $T_S$ (ex.  5 min or 2.5 min) are denoted as $\hat{\tau}_\alpha$ and $\hat{\tau}_d$.

### Connection performance with parameters $\tau_\alpha, \tau_d, v, \sigma$ and with their estimators $\hat{\tau}_\alpha, \hat{\tau}_d, \hat{v}, \hat{\sigma}$

In this part, we assess the performance of the connection algorithm with different parameters $\lambda$ and $\tau_d$. The duration of movies $T_S$ is set to 5 min.  This choice of $T_S$ is a result of some experimental tests.  In the Appendix, The estimation results with 2.5-min movies are presented.  It shows that at $T_S=2.5$ min, the estimators are not accurate enough to ensure the good performance of the connection procedure.  Therefore it justifies the choice of $T_S= 5$ min is a good balance between data acquisition and connection performance.  The results of connection measured by ARI is presented in Fig.  \[fig:Comp\_Ari\_5min\_delta\_t0\_25\]. 

For every combination of $\lambda$ and $\tau_d$, 100 replications of segments connection procedure are performed on 100 5-min movies.  The performance of connection is effected by the number of segments in each movie to be connected.  The higher the density of segments is, the more difficult it is to find the right ones.  One point in the black line represents the mean value of the number of segments of 100 movies corresponding to the setting of $\lambda$ and $\tau_d$.

![ Connection performance comparison for different $\lambda$ and $\tau_d$, when $T_S = 5$ min.  Blue boxes represents ARI values obtained with true parameters $\tau_d, \tau_\alpha, v$ and $\sigma$; Magenta color represents ARI values obtained with estimators $\hat{\tau}_d, \hat{\tau}_\alpha, \hat{v}$ and $\hat{\sigma}$. The black line represents the mean number of segments in a movie.[]{data-label="fig:Comp_Ari_5min_delta_t0_25"}](figs/Comp_Ari_5min_delta_t0_25)

It can be noticed that the ARI value when we use the estimators $\hat{\tau}_d, \hat{\tau}_\alpha, \hat{v}$ and $\hat{\sigma}$, 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 images sequences.  For the highest particle density, when the average number of segments reaches 100, the ARI criterion is still around 0.8.  When the number of segments is around 20, for more than half of the replicates, the ARI values are higher than 0.9. 

### The impact of the intrinsic characteristics of dynamics, the speed $v$ and the variance $\sigma$, to the connection algorithm

In the following, birth rate $\lambda = 0.08$ and death rate $\tau_d=0.004$ are fixed for all the experiments.  Several experiments are conducted, varying the distribution of the speed and the variance.  Two situations for the speed is designed.  In the first one, all particles in one movie have the same speed as 0.6/s (note as *vc* or *v const*) and in the second one, the speed of particles $v\sim \mathcal{U}(0.4,0.8)$ (note as *vu* or *v unif*) in one movie.  The variance $\sigma$ is constant in one movie, however, depending on the experiments it varies from 0.2 to 0.5. 

![ Estimation of $v$ and $\sigma$. $v_x$ (resp.  $v_y$) and $\sigma_x$ (resp.  $\sigma_y$) are the speed and the diffusion in the circumferential (resp.  main axis) direction of the cylinder.  By default $v_y = 0.01v_x$. $\sigma_x=\sigma_y=\sigma$ (isotropic diffusion). ’ vc’ means $v$ is constant and ’vu’ means that $v \sim \mathcal{U}(0.4, 0.8)$.[]{data-label="fig:estim_v_sigma"}](figs/v_sigma_2_5min_delta_t0_25)

In Fig.  \[fig:estim\_v\_sigma\] each box shows the distribution of the estimator in one experiment with 100 replication on 100 generated movies.  In all these four plots, the first four boxes correspond to the case where $\mathbf{v}$ is constant and the last four boxes correspond to the case where $v\sim \mathcal{U}(0.4,0.8)$. In the plot of $\hat{v}_x$, it can be noticed that $\hat{v}_x$ are unbiased when $v$ is constant.  The variances increase as $\sigma$ increases.  However when $v\sim \mathcal{U}(0.4,0.8)$, the red lines on the boxes, which represent the median of the estimator is higher than 0.6.  It should be specified that the estimation of $v$ and $\sigma$ is averaged over the entire movie.  This choice makes the estimation of $v$ more stable when all the particles have the same speed.  However, when the speed is drawn from a random law, the distribution of $\hat{v}_x$ is no longer centered on the expected value and is positively biased.  The explanation is that faster particles can be observed more often in the ROO.  The estimations of $\sigma_x$ and $\sigma_y$ are biased.  The bias might be related to the temporal resolution $\Delta t$.

![ Comparison of segments connection when $v$ and $\sigma$ vary.  Blue (resp.  green, dark green and magenta) boxes represents ARI values obtained with parameters ($\tau_d, \tau_\alpha, \mathbf{v}, \sigma$) (resp. ( $\tilde{\tau}_d, \tilde{\tau}_\alpha, \mathbf{v}, \sigma$), ($\hat{\tau}_d, \hat{\tau}_\alpha, \mathbf{v}, \sigma$) and ($\hat{\tau}_d, \hat{\tau}_\alpha, \hat{\mathbf{v}}, \hat{\sigma}$)).[]{data-label="fig:Comp1_Ari_2_5min_delta_t0_25_v_sigma"}](figs/Comp1_Ari_2_5min_delta_t0_25_v_sigma)

![ Comparison of segments connection when $v$ and $\sigma$ vary.  Blue (resp.  violet, cyan and magenta) boxes represents ARI values obtained with parameters ($\tau_d, \tau_\alpha, \mathbf{v}, \sigma$) (resp. ( $\tau_d, \tau_\alpha, \hat{\mathbf{v}}, \hat{\sigma}$), ($\tilde{\tau}_d, \tilde{\tau}_\alpha, \hat{\mathbf{v}}, \hat{\sigma}$), and ($\hat{\tau}_d, \hat{\tau}_\alpha, \hat{\mathbf{v}}, \hat{\sigma}$)).[]{data-label="fig:Comp2_Ari_2_5min_delta_t0_25_v_sigma"}](figs/Comp2_Ari_2_5min_delta_t0_25_v_sigma)

The connection results show that when we use true values for all parameters in the algorithm, which correspond to blue boxes in both Figs.  \[fig:Comp1\_Ari\_2\_5min\_delta\_t0\_25\_v\_sigma\] and \[fig:Comp2\_Ari\_2\_5min\_delta\_t0\_25\_v\_sigma\], the connection performance decreases considerable and reasonably with $\sigma$ increases.  When we replace all parameters by there estimators ($\hat{\tau}_d, \hat{\tau}_\alpha, \hat{\mathbf{v}}, \hat{\sigma}$), corresponding to magenta boxes in both Figs.  \[fig:Comp1\_Ari\_2\_5min\_delta\_t0\_25\_v\_sigma\] and \[fig:Comp2\_Ari\_2\_5min\_delta\_t0\_25\_v\_sigma\], it can be seen that the last four magenta plots, where the speed $v_x $ is not the same for all the particles, have low ARI values.  To investigate in detail which estimators cause the degradation of the connection performance, we plot the intermediate cases to analyze the individual impact of the estimators (Figs \[fig:Comp1\_Ari\_2\_5min\_delta\_t0\_25\_v\_sigma\] and \[fig:Comp2\_Ari\_2\_5min\_delta\_t0\_25\_v\_sigma\]).  Although in Figs \[fig:hat\_tau\_alpha\_5min\_delta\_t0\_25\] and \[fig:hat\_tau\_d\_5min\_delta\_t0\_25\] it shows that $\tilde{\tau}_\alpha$ and $\tilde{\tau}_d$ (estimators with 30-min movies) have better performance than $\hat{\tau}_\alpha$ and $\hat{\tau}_d$ (estimators with 5-min movies), the connection performance using ($\tau_\alpha, \tau_d$), ($\tilde{\tau}_\alpha, \tilde{\tau}_d$), or ($\hat{\tau}_\alpha, \hat{\tau}_d$) are almost equally good (Fig.  \[fig:Comp1\_Ari\_2\_5min\_delta\_t0\_25\_v\_sigma\] in blue, green and dark green).  we now focus on the last four boxes in Fig.  \[fig:Comp2\_Ari\_2\_5min\_delta\_t0\_25\_v\_sigma\], when we use true value of $\tau_\alpha$ and $\tau_d$, and use the estimators of $\hat{\mathbf{v}}, \hat{\sigma}$, the ARI values drop brutally (from blue to violet).  Combing both figures, it can be concluded that good estimation of $\mathbf{v}$ and $\sigma$ is essential to ensure the good performance of the algorithm. 

Analysis of the connection results
----------------------------------

![ The difference between $K(c^t)$ and $K(c^*)$ versus ARI for different values of birth rate and death rate.[]{data-label="Analyse_ari"}](images/Analyse_ari_true_parameters)

In this section, we evaluated the connection error caused by randomness.  We display in Fig.  \[Analyse\_ari\] the scatter plots of ARI value vs $K(c^t)-K(c^*)$, where $c^*$ denotes the optimal configuration calculated by the “Segment Connection Algorithm” while $c^t$ is the true configuration.  Each scatter plot displays the results of 100 simulations for a given combination of $\lambda$ and $\tau_d$.

The difference between $K(c^t)$ and $K(c^*)$ is always positive or null, showing the optimization procedure works correctly to find the optimal solution.  When $K(c^t)>K(c^*)$, it means that the configuration $c^*$ has higher probability than the true realization $c^t$, which can be due to randomness.  We can notice that ARI decreases as soon as $K(c^t)-K(c^*)$ increases.  This error occurs when the realization is significantly different from the optimal configuration.  Overall, ARI values are generally above 0.7. 

### Number of rotations around the cylinder

![ The distribution of ’number of rotations’ under the assumption that the life duration follows an exponential law of parameter $\tau_d$.[]{data-label="Nb_segs"}](images/Nb_segs)

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 $\tau_d$ and the dynamic velocity $v_x$ are known.  Accordingly, the value of the number of loops is known to be equal to $\frac{v_xT_d}{L},$ where $T_d \sim \mathcal{E}(\tau_d)$ which gives a theoretical expectation value of $\frac{v_x}{\tau_d L}.$ By counting the segments for each trajectory, we can obtain a proxy of the number of rotations around the cylinder.  In Fig.  \[Nb\_segs\], the theoretical values of the number of rotations, after being rounded, are presented in blue color for different values of $\tau_d$ between $0.002$ and $0.005$. Cyan color is used to display the distribution of the true number of segments for each trajectory.  Magenta color is used to display the distribution of the number of segments estimated by our connection procedure.  Overall, when $\tau_d$ is small, the distribution has a heavier tail.  In general, the cyan and magenta distributions are similar to each other. 

Conclusion {#conclusion .unnumbered}
==========

In this paper, we proposed a probabilistic framework and a computational approach with no hidden parameter to connect trajectory segments 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 segments 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 to the design of experiments. 

For future works, we plan to investigate real TIRF (Total Internal Reflection Fluorescence Microscopy [@axelrod1984total]) microscopy 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 such hidden region.  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. 

Appendix: Estimators $\hat{\tau}_\alpha$ and $\hat{\tau}_d$ and the connection performance with 2.5-minute movies {#appendix-estimators-hattau_alpha-and-hattau_d-and-the-connection-performance-with-2.5-minute-movies .unnumbered}
=================================================================================================================

Estimation of birth rate $\tau_\alpha$ {#estimation-of-birth-rate-tau_alpha .unnumbered}
--------------------------------------

The performance of segments connection in 2.5-min movies is evaluated and the result is presented in Fig.  \[fig:hat\_tau\_alpha\_2\_5min\_delta\_t0\_25\].  Compared to Fig.  \[fig:hat\_tau\_alpha\_5min\_delta\_t0\_25\] where 5-min movies are used, the estimations with 2.5-min movies have bigger variance.  With 2.5-min movies, there are more estimations equal to zero. 

![ The estimation of arrival rate $\tau_\alpha$ with different $\lambda$ and $\tau_d$. Magenta boxes represents $\hat{\tau}_\alpha$ estimated by 2.5-min movies and blue boxes by 30-min movies.[]{data-label="fig:hat_tau_alpha_2_5min_delta_t0_25"}](figs/tau_alpha_2_5min_delta_t0_25)

Estimation of death rate $\tau_d$ {#estimation-of-death-rate-tau_d .unnumbered}
---------------------------------

![ The estimation of death rate $\tau_d$ with different $\lambda$ and $\tau_d$. Magenta boxes represents $\hat{\tau}_d$ estimated by 2.5-min movies and blue boxes by 30-min movies.  Black lines represent the true value of $\tau_d$.[]{data-label="fig:hat_tau_d_2_5min_delta_t0_25"}](figs/tau_d_2_5min_delta_t0_25)

It is shown that both in Figs \[fig:hat\_tau\_alpha\_2\_5min\_delta\_t0\_25\] and \[fig:hat\_tau\_d\_2\_5min\_delta\_t0\_25\], the estimators $\hat{\tau}_\alpha$ and $\hat{\tau}_d$ have many zero values.  This phenomenon will deny the birth and death of particles, and force the connection of observed segments.  The result of connection is presented in the next section and it shows that these biases of parameter estimation have severe influence to the segments connection performance. 

Evaluation of the connection procedure {#evaluation-of-the-connection-procedure .unnumbered}
--------------------------------------

![ Connection performance comparison.  Color representation is the same as in Fig.  \[fig:Comp1\_Ari\_2\_5min\_delta\_t0\_25\_v\_sigma\].  Light green boxes, representing the result when $\tilde{\tau}_\alpha, \tilde{\tau}_d$ (estimators with 30-min movies) and true value $v, \sigma$ are used, shows performance as good as the blue boxes, where true parameters values are used.  However, the dark green boxes, where $\hat{\tau}_\alpha, \hat{\tau}_d$ (with 2.5-min movies) and $v, \sigma$ are used, show much degraded results.[]{data-label="ARI_taux_v_comp1"}](figs/Comp1_Ari_2_5min_delta_t0_25)

![ Connection performance for different $\lambda$ and $\tau_d$, when $v_x$ = 0.6 and $\sigma = 0.2$. Color representation is the same as in Fig.  \[fig:Comp2\_Ari\_2\_5min\_delta\_t0\_25\_v\_sigma\].  Boxes with blue, violet, and cyan colors are similar, which means that the estimators $\hat{v}, \hat{\sigma}, \tilde{\tau}_\alpha$ and $\tilde{\tau}_d$ do not cause degradation of the connection.  Only when $\hat{\tau}_\alpha$ and $\hat{\tau}_d$ are used, shown in magenta boxes, the connection results degrades.[]{data-label="ARI_taux_v_comp2"}](figs/Comp2_Ari_2_5min_delta_t0_25)

Figs.  \[ARI\_taux\_v\_comp1\] and \[ARI\_taux\_v\_comp2\] show the results of segments connection measured by ARI, according to different settings of birth rate $\lambda$ and death rate $\tau_d$. Concerning the experiment with true parameters (blue boxes in both figures), the more the number of segments, the less precise the connection performance is.  However, when the estimators $\hat{\tau}_d, \hat{\tau}_\alpha$ are used (dark green boxes and magenta boxes), and the number of segments is too small, the errors introduced by the estimator cause the crush of the connection algorithm. 

Through Fig.  \[ARI\_taux\_v\_comp1\] and \[ARI\_taux\_v\_comp2\], it can be concluded that the error of estimators $\hat{\tau}_d, \hat{\tau}_\alpha$ is the main cause of the dramatic decrease of ARI when all the estimators $\hat{\tau}_d, \hat{\tau}_\alpha, \hat{v}$ and $\hat{\sigma}$ are used. 

The accuracy of $\hat{\tau}_\alpha$ increases as the total observed time $T_S$ increases (see Fig.  \[fig:tau\_alpha\_time\]).  Therefore, in the study we chose movies of 5 min to estimate $\tau_\alpha$ and $\tau_d$ and then to evaluate the connection performance. 

Acknowledgements {#acknowledgements .unnumbered}
================

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.