Philippe Roudot
email: philippe.roudot@gmail.com

Charles Kervrann

Cedric Blouin

Francois Waharte

Lifetime estimation on moving sub-cellular objects in frequency domain FLIM imaging

come

Lifetime estimation on moving sub-cellular objects in frequency domain FLIM imaging

INTRODUCTION

The fluorescence lifetime of a given fluorescent species is defined as the average time during which the fluorophore stays in excited state before relaxing to its ground state and possibly emitting a photon. The fluorophore energy can be relaxed to its ground state through various de-excitation pathways. Those pathways are influenced by various biophysical phenomena such as ion concentration, collision or non-radiative energy transfer. Exploiting those properties, Fluorescence Lifetime Imaging Microscopy (FLIM) is now a commonly-used technique for sensing fluorophore environment in a living biological sample (pH, ions, temperature ...). Fluorescence lifetime measurement is also particularly useful to detect Förster resonance energy transfer (FRET) between fluorescent molecules in close proximity [START_REF] Jares-Erijman | Fret imaging[END_REF][START_REF] Lakowicz | Principles of fluorescence spectroscopy[END_REF].

In this area, two main technologies can be considered: time domain methods and frequency domain methods. Time domain methods such as TCSPC (Time-Correlated Single-Photon Counting) measure the mean delay between a photon emission and a pulsed excitation. This approach is now well established but mainly recommended for fixed sample with no motion during acquisition [START_REF] Birch | Time-domain fluorescence spectroscopy using time-correlated single-photon counting[END_REF][START_REF] Chessel | A maximum likelihood method for lifetime estimation in photon counting-based fluorescence lifetime imaging microscopy[END_REF]. Unlike TCSPC, a FLIM measurement in the frequency domain is defined as a short sequence of images presenting a sinusoidal intensity footprint that reflects the sinusoidal excitation of the sample followed by the phase modulation of the emitted signal. A fluorescence lifetime variation results in a phase variation in the detection signal. The acquisition process is 10 to 500 times faster than time domain methods [START_REF] Lajevardipour | The effect of translational motion on flim measurements-single particle phasor-flim[END_REF][START_REF] Schneider | Rapid acquisition, analysis, and display of fluorescence lifetime-resolved images for real-time applications[END_REF][START_REF] Holub | Fluorescence lifetime imaging (fli) in real-time -a new technique in photosynthesis research[END_REF]. FD FLIM is then intrinsically more suitable for dynamical process analysis and for the study of living cells.

In FD FLIM, Fourier transform or equivalently least mean squares fitting of a sine function is used traditionally to estimate lifetime maps from FD FLIM measurements [START_REF] Lakowicz | Principles of fluorescence spectroscopy[END_REF][START_REF] Spring | Image analysis for denoising full-field frequency-domain fluorescence lifetime images[END_REF]. Two problems arise from those approaches. First, it is implicitly assumed that the noise induced during acquisition is temporally stationnary and Gaussian. This traditional assumption does not hold because of the photonic nature of the emitted signal and the high intensity gradient present in FD FLIM measurement. Secondly, the specimen is assumed to remain immobile during acquisition. Even if frequency domain methods are faster than time domain methods, transient motions of intracellular structures can still happen during FD FLIM acquisitions and thus inducing artifacts in the estimated lifetime map [START_REF] Lajevardipour | The effect of translational motion on flim measurements-single particle phasor-flim[END_REF][START_REF] Hanley | Fluorescence lifetime imaging: multi-point calibration, minimum resolvable differences, and artifact suppression[END_REF]. In this paper, we propose to overcome those two difficulties to address important biological issues in dynamic intra-cellular traffic experiments. Our line of work consists in exploiting the usual FD FLIM intensity model combined with the original noise modeling described in [START_REF] Roudot | Noise modeling for intensified camera in fluorescence imaging: application to image denoising[END_REF] to estimate accurately fluorescence lifetime on static background and on moving vesicles that underwent non-negligible motions during the acquisition. In what follows, we propose a novel and original statistical framework based on robust statistics to estimate alternately the trajectory and the fluorescence lifetime of moving vesicles. To our knowledge, this issue is not addressed in the literature and this paper describes an innovative solution to solve this problem. The remainder of the paper is organized as follows. Section 2 presents the principles of fluorescence lifetime imaging techniques in the frequency domain (FD FLIM). Section 3 describes our estimators and correction terms for fluorescence lifetime on the static background. Section 4 describes our lifetime estimation approach on dynamical structures. In Section 5, the proposed methodology is evaluated on experimental image series with controlled fluorophores and simulated measurements. On experimental data, it is extremely difficult to obtain a ground truth for the fluorescence lifetime of moving vesicles. Accordingly, the accuracy of our method is tested via the motion estimation results that are compared against competitive particle trackers. Finally, the complete lifetime map reconstruction framework is applied to the study of a GFP-tagged donorreceptor expressed in RPE1 cells.

S θ (k) = G ccd (F (t)M k (t)) = C S + A S cos 2πk K + φ τ S θ (k) = G ccd (F (t)M k (t)) = C S + A S cos 2πk K + φ τ

FD FLIM AND SIGNAL MODELING

A. Basics in frequency domain FLIM measurement

Instead of the pulsed laser excitation that characterizes timedomain methods, the FD FLIM excitation source is a sinusoidally modulated signal E(t) (blue ray Figure 19) defined as:

E(t) = C E + A E sin(ω E t + φ E ) (1)
where ω E denotes the radial frequency of the excitation signal, C E the offset, A E the amplitude, φ E is the phase delay and t the temporal variable. The excitation E(t) can be operated with various illumination methods, such as wide-field or confocal microscopy. The signal emitted by the sample F(t) (green ray in Figure 19) is defined as :

F(t) = E(t) * R(t) (2) = C F + A F sin(ω E t + φ E -atan(ω E τ))
where * denotes the convolution operator, R(t) = R 0 e -t/τ is the response of a fluorescent sample at time t with a given lifetime τ and R 0 is the response at time t = 0. Accordingly, the phase delay φ Eatan(ω E τ) allows us to recover the fluorescence lifetime τ.

Experimentally, for measuring lifetime in the order of nanoseconds (ns) the frequency ω E has to be > 1 MHz, too high to be compatible with the frame rate of a CCD sensor. Therefore the signal F(t) is phase-modulated with a finite set of K sinusoidal signals M k (t), k ∈ [1, K] (see black signal in Figure 19) at the same frequency ω E (homodyne detection) defined as:

M k (t) = C M + A M sin ω E t + φ M 0 + 2πk K ( 3
)
where C M , A M and φ M 0 denote respectively the offset and the amplitude of M k and the phase of signal M k for k = 0. This modulation is operated by applying a time-varying voltage to the photocathode terminal in the intensified CCD. The intensified CCD is in turn required to acquire the lowest intensity images resulting from modulation. Thanks to the low pass effect of the CCD detector, the higher frequencies of the phase modulated-signal F(t)M k (t) are attenuated and become negligible. After phase-modulation, the K time-independent signals are described as a function of k ∈ [1, K] :

S θ (k) = G CCD (F(t)M k (t)) = C S + A S cos 2πk K + φ τ ( 4
)
where G CCD represents the CCD function, In what follows, we will denote θ = (C S , A S , φ τ ) T ∈ R 3 the FD FLIM parameters at a given pixel. FD FLIM thus consists in K images presenting a sinusoidal footprint (see Figure 2.A). For each pixel in the image domain, the phase φ τ combined with proper calibration yields the lifetime estimate τ. A detailed derivation for the fluorescent signal expression and the phase modulation can be found in Appendix A.

C S = C F C M , A S = A F A M and φ τ = (φ M 0 -φ E + atan(ω E τ)) is the phase.

B. Conventional estimation of FD FLIM parameters

Estimating fluorescence lifetime from FD FLIM measurements amounts to estimating the parameters that control the sinusoidal footprint (see Figure 2

.B).

Let us denote a measured image sequence of K images I : 2 is the image domain. Let φ τ : Ω → R + and τ : Ω → R + be respectively the unknown phase and lifetime map. The usual approach to estimate the phase φ τ (x) is based on the Fourier transform and is defined as:

Ω × [1, K] → R + where Ω ⊂ R
φ τ (x) = -atan ∑ K k=1 I(x, k) cos( 2πk K ) ∑ K k=1 I(x, k) sin( 2πk K ) + π 2 , (5)
where I(x, k) is the intensity measured at pixel x in the phasemodulated image k. The large majority of contributions in FD FLIM exploits this estimator [START_REF] Spring | Image analysis for denoising full-field frequency-domain fluorescence lifetime images[END_REF][START_REF] Holub | Fluorescence lifetime imaging at video rate a new technique in photosynthesis research[END_REF][START_REF] Van Munster | Suppression of photobleachinginduced artifacts in frequency-domain flim by permutation of the recording order[END_REF][START_REF] Hutchinson | Fluorescence lifetime-based sensing in tissues: a computational study[END_REF] as much as commercial software such as LI-FLIM by Lambert Instrument. In order to handle the impact of photobleaching, [START_REF] Van Munster | Suppression of photobleachinginduced artifacts in frequency-domain flim by permutation of the recording order[END_REF] proposed to randomize the order of phase modulation. The noise induced by photobleaching is thus shown to be mixed into the uncorrelated measurement noise. Let us denote the true signal S θ : Ω × [1, K] → R + where S θ (x, k) is the signal value at a given pixel x ∈ Ω and frame k ∈ [1, K]. The equivalent statistical model put forward in those approaches is as follows:

I(x, k) = S θ (x, k) + ε(x, k), (6)
where ε(x, k) is assumed to be a white Gaussian noise such as E[ε(x, k)] = 0 and Var[ε(x, k)]=υ 2 . The associated optimal estimator is found as the minimization of the following least mean squares problem:

θ(x) = argmin θ(x) K ∑ k=1 (I(x, k) -S θ (x, k)) 2 . (7)
Very few studies have been focusing on improving this phase estimator. For instance, in [START_REF] Gadella | Fluorescence lifetime imaging microscopy: Pixel-by-pixel analysis of phase-modulation data[END_REF][START_REF] Clegg | Lifetime-resolved fluorescence imaging[END_REF], the authors proposed a linearization of the problem in the Fourier domain. Finally, given a two-dimensional map of phase estimation, the fluorescence lifetime map at pixel x ∈ Ω is given by (see Appendix A.2):

τ(x) = tan φ τ (x) -φ M 0 + φ E ω E . ( 8
)
Note that the phase delay (φ M 0 -φ E ) of the system must be calibrated by recording a reference sequence of a fluorescent sample with known lifetime τ re f while the frequency ω E is controlled during acquisition.

ROBUST ESTIMATION OF FD FLIM PARAMETERS OVER STATIC BACKGROUND

As described above, least mean squares minimization assumes that the errors are expected to be Gaussian distributed and homoscedastic. While the Gaussian approximation is correct in most cases due to the high counts of photons after intensification, the homoscedastic assumption does not hold in FD FLIM imaging because of photon counting noise, large intensity gradient and the specificity of the FD FLIM image formation process (see Figure 3). In this Section, we describe a novel estimator based on our previous work on ICCD noise analysis that takes this heteroscedasticity into account and provide a complementary scheme to detect outlier induced by small object motions.

A. Noise variance modeling and estimator

In [START_REF] Roudot | Noise modeling for intensified camera in fluorescence imaging: application to image denoising[END_REF], we described a study of instrument-induced noise based on sensor theory combined with experimental consideration. We showed that the quadratic relationship between the noise variance and intensity level demonstrated in [START_REF] Boddeke | Quantitative fluorescence microscopy: autofocusing, z-axis calibration, image sensors, fluorescence lifetime imaging[END_REF] is strongly perturbed by optical aberration induced by the intensifier setup. We thus proposed in our previous work a correction term to the parametric noise variance model υ 2 (x, k). The new model for the noise variance is thus defined as the product of two terms at pixel x and frame k:

υ 2 (x, k) = (aE[I(x, k)] 2 + bE[I(x, k)] + c) e - (x-x 0 ) 2 hx - (y-y 0 ) 2 hy + o . ( 9
)
where E[•] denotes the mathematical expectation and a, b, c are the coefficients of a quadratic polynome. The first term describes the overall contribution of the photon counting noise, read out noise and intensifier variability as described in [START_REF] Roudot | Noise modeling for intensified camera in fluorescence imaging: application to image denoising[END_REF][START_REF] Boddeke | Quantitative fluorescence microscopy: autofocusing, z-axis calibration, image sensors, fluorescence lifetime imaging[END_REF]. The second term is a spatial correction term determined empirically to cope with spherical aberration. In [START_REF] Roudot | Noise modeling for intensified camera in fluorescence imaging: application to image denoising[END_REF], we have proposed a Gaussian function parametrized by x 0 , y 0 , h x , h y and an offset o to take the blur observed on real images into account. The blurring effect increases with the distance from the center x 0 = (x 0 , y 0 ) T of the aberration (the center of the lens inside the intensifier). The x 0 , y 0 , h x and h y parameters are estimated from the same plain fluorescein sample used for phase calibration (see [START_REF] Roudot | Noise modeling for intensified camera in fluorescence imaging: application to image denoising[END_REF] for details).

Equipped with this non-stationary noise model, the FD FLIM parameters over the background are estimated by minimizing the following weighted least mean squares criterion:

θ(x) = argmin θ K ∑ k=1 I(x, k) -S θ (x, k) 2 υ 2 (x, k) ( 10
)
where S θ (x, k) is the signal defined as in [START_REF] Chessel | A maximum likelihood method for lifetime estimation in photon counting-based fluorescence lifetime imaging microscopy[END_REF]. The potential of this weighted least mean squares estimation method is discussed in Section 5.

B. Initialization and convergence

The initialization for the parameter θ(x) takes into account the set of measurements as well as the expected value taken by the fluorescent sample. Indeed, while the offset the amplitude of the sinusoidal signal are predicted using the mean and the range of the measure intensities, the phase is initialized using the lifetime of the fluorophore in the absence of FRET and the calibrated phase of the system. More formally, the initialization is carried out as follows:

⎧ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎨ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎩ C S (x) = 1 K K ∑ k=1 I(x, k) A S (x) = 1 2 max k∈[1,K] I(x, k) -min k∈[1,K] I(x, k) φ τ (x) = φ M 0 -φ E + atan(ω E τ re f )
The evaluation of the magnitude of our weighted residuals over a large parameter space constrained by φ τ (x) ∈ [0, 2π] highlights a single minimum (see Figure 4). Given our efficient initialization procedure, a conventional Gauss-Newton algorithm converges quickly towards the global minimum.

C. Motion-induced outlier management

Assuming a lifetime map τ : Ω → R + computed from ( 8)-( 10), we propose here a complementary approach to estimate the background lifetime in pixels "occluded" by moving objects (see Figs. 6 and12). To tackle this problem, we analyze the value τ(x) with respect to the lifetime values τ(y) estimated at pixels y in the neighborhood V(x) of x. We appeal to tools from robust statistics to determine how large the lifetime τ(x) can be before we consider this value to be an "outlier". We define the Median Absolute Deviation κ τ (x) to determine an adaptive threshold (see [START_REF] Rousseeuw | Robust Regression and Outlier Detection[END_REF]):

κ τ (x) = 1.4826 × median y∈V(x) | τ(y) -median y∈V(y) τ(y)|, (11) = 1.4826 × MAD ( τ(y), y ∈ V(x)),
where "MAD" denotes the median absolute deviation and V(x) ⊂ Ω denotes a spatial neighborhood (or window) cen-tered at location x. The constant in ( 14) is calculated from the fact that the MAD of a zero-mean normal distribution with unit variance is 0.6745 = 1/1.4826. Given an arbitrary real value λ > 0, we can detect outliers by using the usual Chebyshev's inequality which can be applied to any arbitrary distribution:

Prob (|T (x) -μ τ (x)| ≥ λκ τ (x)) ≤ 1 λ 2 (12)
where T (x) denotes a random variable (representing lifetime at pixel x) with finite (robust) expected value:

μ τ (x) = median y∈V(x) (τ(y)) (13)
and finite non-zero variance κ τ . Setting λ = 3 as chosen in all our experiments, means that the probability that the values τ(x) lying outside the interval

[μ τ (x) -3κ τ (x), μ τ (x) + 3κ τ (x)]
does not exceed 88.89% and the probability for a pixel to be classified as an outlier is 11.11%. The neighborhood size |V(x)| is typically set to 30×30 pixels to account for large objects.

Let O : Ω → {0, 1} be the binary "outlier" map such that:

O(x) = ⎧ ⎨ ⎩ 1 if | τ(x) -μ τ (x)| > λκ τ (x), 0 otherwise. ( 14
)
For each pixel x ∈ Ω such that O(x) = 1 ("outlier"), lets denote X NESW = {x N , x E , x S , x W } the set of 4 closest "inlier" lifetime locations along respectively the North, East, South and West directions such as

O(x N ) = O(x E ) = O(x S ) = O(x W ) = 0.
The fluorescence lifetime is then linearly interpolated as:

τ b (x) = ⎧ ⎪ ⎨ ⎪ ⎩ 1 C(x) ∑ y∈X NESW w(x, y) τ(y) if O(x) = 1 and x = y, τ(x) otherwise, with C(x) = ∑ y∈X NESW w(x, y), w(x, y) = x -y -2
2 and • 2 denotes the Euclidean norm. This interpolation technique can be very efficiently implemented as a two-pass filter and a demonstration of the method is illustrated in Figure 12 (simulation) and Figure 14 (experimental data). A more sophisticated but more computationally intensive interpolation method could be also considered as in [START_REF] Pécot | Background fluorescence estimation and vesicle segmentation in live cell imaging with conditional random fields[END_REF]. In this approach, the idea is to exploit the set of pixels located at the periphery of each connected component of the binary map O to compute a weighted average lifetime.

In the next section, we describe a fluorescence lifetime estimation method dedicated to moving objects and particles.

ESTIMATION OF FD FLIM PARAMETERS ON MOVING OBJECTS AND ENDOSOMES

One of the advantage of FD FLIM over time domain techniques is the speed of operation that puts time lapse measurements of fluorescence lifetime on dynamic sample within reach. However, in the conventional estimation methods presented in Section 2 B as well as the approach proposed above, the fluorescent objects are supposed to be immobile. Nevertheless, biological processes often involve motion at various temporal and spatial scales, and often faster than the acquisition time of a lifetime map (a few seconds typically). As a result, motions can occur on living cells during the phase modulation. Figure 5 illustrates the impact of motions on the phase estimation. The lifetime misestimation results in a so-called "doppler effect" on the lifetime map as illustrated in Figure 6 with an artefactual gradient in the lifetime map in the direction of motion. Indeed, when analyzing the intensity over the phase sequence to determine the lifetime, deviations from a sine function can be observed (see Figure 5 A) due to the object motion.

The issue has been previously diagnosed in [START_REF] Hanley | Fluorescence lifetime imaging: multi-point calibration, minimum resolvable differences, and artifact suppression[END_REF] and further analyzed in [START_REF] Lajevardipour | The effect of translational motion on flim measurements-single particle phasor-flim[END_REF] where the impact of moving object on lifetime variance is measured. However, to our knowledge, no method has been proposed to estimate the fluorescence lifetime on moving structures.

In this section, we propose a statistical framework for robust lifetime estimation on moving endosomes in order to provide a corrected reconstruction of the lifetime map τ. To alleviate complexity while taking into account the correlation of parameters, we propose an iterative estimation procedure where the estimation of FD FLIM and lifetime parameters and motion parameters are decoupled. The proposed method needs only the setting two non-critical input parameters: an upper-bound on the object number and the diffracted endosome size (e.g. a square window of 7 × 7 pixels is a typical choice and the scale may vary on microscope magnification). Fig. 6. "Doppler effect" on lifetime measurement (color bar in nanoseconds) due to motion (see Figure 5). Scale bar is 10 μm.

A. Modeling spatial and temporal variations

In order to model the lifetime of a moving object, we propose a parametrization of both spatial and temporal intensity variations of a sub-resolved object. In our approach the microscope point spread function (PSF), corresponding to a Bessel function, is approximated as a Gaussian function as it is a common practice in spot detection literature [START_REF] Smal | Quantitative comparison of spot detection methods in fluorescence microscopy[END_REF]. Secondly the signal is composed of a background and a pinpoint vesicle before microscope diffraction. The background follows the model ( 4) due to the membrane imaging and/or cytoplasmic fluorescence. Finally, the background is assumed to be smooth over the PSF support.

More formally, S(x, k) is the signal defined at location x ∈ Ω (pixel grid) as the sum of N object signals (e.g. vesicle, endosomes, ...) and a background signal:

S(x, k) = S θ b (x, k) + N ∑ i=1 S θ i (x i , k)e - x-x i,k 2 2 2σ 2 i,k , x ∈ Ω ( 15
)
where S θ b (x, k) ∈ R + defined in (4) is the spatially varying background signal parametrized by

θ b (x) = (C S b (x), A S b (x), φ τ b (x)) T , S θ i (k) is the i-th object intensity parametrized by θ i = (C S i , A S i , φ τ i ) T . In this modeling, σ i,k
is the standard deviation representing the isotropic Gaussian shape of object i at frame k and x i,k is the location of object i at frame k. Given an object trajectory {x i,1 , • • • , x i,K } over K frames, we consider the following approximation for x ∈ B(x i,k ) ⊂ Ω where B(x i,k ) is a local neighborhood (circle or square window) centered at x i,k whose size is constant for all objects and for all frames in our study:

S(x, k) ≈ S θ i (x i,k ) e - x-x i,k 2 2 2σ 2 i,k , x ∈ B i,k . (16)
In other words, we assume that the signal of an object i along its trajectory dominates locally other signal contributions.

B. Estimation procedure

Assuming a given set of images

I 1:K = {I(1, .) : Ω → R + , • • • , I(K, .) : Ω → R + }.
We must first detect the N objects (e.g. vesicles) at locations {x 1,1 , • • • , x N,1 } on the first frame (k = 1). Object detection can be performed by using any robust algorithms (e.g. see [START_REF] Smal | Quantitative comparison of spot detection methods in fluorescence microscopy[END_REF][START_REF] Olivo-Marin | Extraction of spots in biological images using multiscale products[END_REF][START_REF] Thomann | Automatic fluorescent tag detection in 3d with super-resolution: application to the analysis of chromosome movement[END_REF][START_REF] Basset | SLT-LoG: A vesicle segmentation method with automatic scale selection and local thresholding applied to TIRF microscopy[END_REF]). We have considered the method #10 in [START_REF] Chenouard | Objective comparison of particle tracking methods[END_REF] based on structure tensors [START_REF] Rao | Computing oriented texture fields[END_REF] and an optimal histogram based thresholding [START_REF] Matusita | Decision rule, based on the distance, for the classification problem[END_REF] since this combination of algorithms was able to reliably Given the initial set of objects, the trajectory x i,1:

K = {x i,1 , • • • , x i,K }, the shape parameter σ i,1:K = {σ i,1 , • • • , σ i,K }
and the parameters θ i of object i are individually estimated using an iterative procedure described in the next section.

B.1. Maximum likelihood estimation framework

The basic idea of estimation for a given object indexed by i is to consider the parameter x i,1:K , σ i,1:K and θ i on the same level. Let p θ i (x i,1:K , σ i,1:K , I 1:K ) the joint probability distribution of unknown variables and observations (i.e. measurements). The Generalized Maximum Likelihood (GML) estimation problem [START_REF] Mohammad-Djafari | On the estimation of hyperparameters in bayesian approach of solving inverse problems[END_REF] amounts to optimizing the following criteria:

( x i,1:K , σ i,1:K , θ i ) = arg max x i,1:K ,σ i,1:K ,θ i p θ i (x i,1:K , σ i,1:K , I 1:K ). (17)
Under the proper assumptions of independence, the probability distribution can be generally factorized as follows:

p θ i (x i,1:K , σ i,1:K , I 1:K ) ∝ K ∏ k=1 ∏ x∈B(x ik ) p θ i (x, x i,k , σ i,k , I(x, k)). ( 18
)
This optimization problem can be implemented by successive and alternative minimization of criteria with respect to θ i and (x i,1:K , σ i,1:K ) as:

⎧ ⎪ ⎪ ⎨ ⎪ ⎪ ⎩ θ (m) i = arg min θ i -log p θ i ( x (m) i,1:K , σ (m) i,1:K , I 1:K ), ( x (m+1) i,1:K , σ (m+1) i,1:K ) = arg min x i,1:K ,σ i,1:K -log p θ (m) i (x i,1:K , σ i,1:K , I 1:K ). ( 19
)
where m denotes the iteration number. Those minimizations, described in the following subsection, are carried out iteratively until convergence of the parameter of interest |( θ

(m) i - θ (m+1) i )/ θ (m) i | ≤
where is a constant (set to 10 -6 in our experiments). Given θ i , the decay time τ i for each detected vesicle is then estimated as (see Appendix A.2):

τ i = tan( φ τ i -φ M 0 + φ E ) ω E . ( 20
)
The final estimated value τ i of each object i is reported on the static lifetime map as a circular patch of value σ i,0 at location x i,0 for each vesicle.

B.2. Iterative and robust estimation of parameters

In the case of Gaussian model ( 10) with ( 18) and ( 16), we have

p θ i (x i,1:K , σ i,1:K , I 1:K ) ∝ exp ⎛ ⎝ - K ∑ k=1 ∑ x∈B(x i,k ) ρ r θ i (x, x i,k , σ i,k ) ⎞ ⎠ (21)
with

r θ i (x, x i,k , σ i,k ) = I(x, k) -S θ i (x i,k )e - x-x i,k 2 2 2σ 2 i,k υ(x, k) . ( 22
)
If the distribution of residuals {r θ i (x, x i,k , σ i,k )} in the neighborhood B(x i,k ) is a zero-mean Gaussian, then ρ(z) = z 2 , z ∈ R equivalent to a weighted least-mean square estimator. To increase robustness and reject possible outliers due to deadpixels or model approximation , we consider a ρ-function that increases less rapidly than z 2 and becomes constant beyond a fixed point determined by a "scale" parameter. The ρ-function minimizes the effect of large residuals induced by model approximations and ρ(z) ≈ z 2 for residuals smaller than the "scale" parameter. Our experiments showed that the Leclerc's influence function defined as ρ(z ; s) = 1e -z 2 /s 2 where s > 0 denotes the scale parameter, gives more stable results while compared with other M-estimators (e.g. Tukey Bi-weighted or German-McClure).

Object tracking. At iteration m of the algorithm, (x i,1:K , σ i,1:K ) are obtained from ( 19)-( 21) as follows:

( x (m) i,1:K , σ (m) i,1:K ) = arg min (x i,1:K ,σ i,1:K ) K ∑ k=1 ∑ x∈B(x i,k ) ρ r θ (m-1) i (x, x i,k , σ i,k ) ; s (m-1) i,k . (23)
Estimation can be performed frame-by-frame and the parameters at frame k are given by:

( x (m) i,k , σ (m) i,k ) = arg min (x i,k ,σ i,k ) ∑ x∈B(x i,k ) ρ r θ (m-1) i (x, x i,k , σ i,k ) ; s (m-1) i,k . (24)
The optimization ( 24) is thus performed using the iteratively reweighted least square method. The scale parameter is updated at each frame k as the median absolute deviation of the residuals as (see [START_REF] Gadella | Fluorescence lifetime imaging microscopy: Pixel-by-pixel analysis of phase-modulation data[END_REF] and [START_REF] Rousseeuw | Robust Regression and Outlier Detection[END_REF]):

s (m) i,k = λ × 1.4826 × MAD (r θ (m-1) i x, x (m) i,k , σ (m) i,k 2 , y ∈ B(x (m) i,k ) (25)
and λ is set to 3 as explained in Section 3 according to the Chebychev's inequality. For each frame k, the parameters are initialized with the estimators obtained on the previous frame:

σ (m) i,k+1 = σ (m) i,k , x (m) i,k+1 = x (m) i,k for k ∈ [1, K].
Let us note that in this framework, the motion is only a "nuisance parameter" and it is not required to track a object along the full FD FLIM measurement to estimate the fluorescence lifetime.

FD FLIM parameter estimation.

Given an estimation ( x

(m) i,1:K , σ (m)
i,1:K ), the FD FLIM parameters θ i are estimated at iteration m as:

θ (m) i = arg min θ i K ∑ k=1 ∑ x∈B(x (m) i,k ) ρ r θ i (x, x (m) i,k , σ (m) i,k ) ; s (m) i,k . (26)
In this case, no model approximation is carried out and no outliers has been detected experimentally. To save time and ease the estimation process, the Leclerc ρ-function is then approximated by a quadratic function, i.e. ρ(z) = z 2 ("small" residuals), and ( 26) is equivalent to a weighted least squares minimization problem. The solution can be then found by using a conventional Gauss-Newton algorithm as explained in Section 3.

B.3. Initialization of the estimation procedure

In our experiments, we considered small neighborhood |B(x i )| of about 7 × 7 pixels.

Initialization of object tracking parameters.

In order to initialize the iteration process, each object detected on the first frame is tracked over the K frames of the FD FLIM measurement using simple frame-by-frame ordinary least-square model fitting over a local neighborhood B(x i,k ) centered at x i,k . In this initialization step, we have considered the following intensity model:

S(x, k) ≈ a(x i,k )e - x-x ik 2 2 2σ 2 i,k + b(x i,k ), x ∈ B i,k , (27)
whose parameters a(x i,k ) and b(x i,k ) are constant values for x ∈ B(x i,k ), jointly estimated at each frame k with the parameters x i,k and σ i,k . We denote x 0)

i,k and σ (0) i,k those parameters which are used as initializations of [START_REF] Rao | Computing oriented texture fields[END_REF].

Initialization of FD FLIM parameters.

At each iteration m, θ i = (C S i , A S i , φ τ i ) T is initialized using the intensity followed by the tracked object and physical property of the fluorescent donor. More specifically, the mean and amplitude of the normalized intensity of each spot over the object support are used to initialize C S i and A S i respectively. The phase φ τ i is initialized by using the fluorophore lifetime τ re f in the absence of FRET. The fluorescence lifetime is determined from the literature or experimentally on fixed samples (e.g. 2.5 ns for the EGFP). For object i, (φ M 0 -φ E ) is obtained by calibration of the optical system. The computation is carried out as follows:

⎧ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎨ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎩ C S i = 1 |B(x i )|K K ∑ k=1 ∑ x∈B(x (m) i,k ) Ĩ(x, k) A S i = max x∈B(x (m) i,k ), k∈[1,K] Ĩ(x, k) -C S i φ τ i = φ M 0 -φ E + atan(ω E τ re f ) where Ĩ(x, k) = I(x, k)e - x-x (m) i,k 2 2 2(σ (m) i,k ) 2 .

C. Summary of the algorithm

The iterative framework is summarized Figure 8. In a nutshell, the full lifetime map reconstruction algorithm is composed of noise calibration, static lifetime estimation and dynamical lifetime estimation on intracellular structure as given in Algorithm 1. Estimations converge after three to five iterations.

EXPERIMENTAL RESULTS

We tested our estimation method on both simulated and experimental sequences. As there can be no ground truth for the fluorescence lifetime on vesicles in experimental sequences, we have been using diverse strategies depending on the scenario at hand. To test our static phase estimator, we exploit a plain reference stack acquisition of fluorescein where a single mode in the phase distribution is expected. To evaluate the performances of dynamic object lifetime estimation, we use simulated data and compare the motion estimation with competitive multiple particle tracking algorithm.

A. FD FLIM setup

FLIM measurements were done in frequency domain by phase modulation on a custom system exploiting a third generation intensifier from Lambert Instruments (model II18MD). The module was attached to a NikonTE2000 TIRF inverted microscope equipped with a Coolsnap HQ CCD camera (Photometrics), a 473 nm modulated laser diode (Omicron) used for excitation of the donor fluorophore and a 100x 1.49 NA TIRF objective. The laser light was coupled using an optical fiber to the TIRF illumination arm or to a Yokogawa CSU10 spinning-disk module. Fluorescence emission was selected by a band-pass filter (500-550 nm). FLIM images were acquired either in widefield illumination using the TIRF illumination set to send the laser perpendicularly to the objective, or in confocal mode using the spinningdisk module. Transfected cells were also imaged in widefield illumination by a mercury lamp with standard Nikon filter cubes and a Coolsnap Ez CCD camera (Photometrics) prior to FLIM measurements in order to estimate the expression level of both GFP-and mCherry-tagged proteins. All the experiments were performed at 37°C.

B. Lifetime map estimation on static background

B.1. Phase estimation on a reference FD FLIM measurement

Figure 10 shows the phase distribution estimated on a FD FLIM acquisition of a sample of fluorescein presenting a single fluorescence lifetime that is constant throughout the sample. As a result, a single mode is expected for the estimated phase distribution. The fluorescence sample has been acquired with a wide-field setup. On this graph, we compared the Fourier decomposition used in the literature (or equivalently, the least mean squares estimator) and our heteroscedastic modeling approach. Those results highlights the significant shift caused by the weighting of least squares residuals (means are 2.46 ns for Fourier vs 2.45 ns, modes are the same) and difference in standard deviation (0.021 ns for Fourier vs 0.015 ns) . Also, our correction results in a 10% gain in kurtosis which is expected for a plain fluorescein sample [START_REF] Spring | Image analysis for denoising full-field frequency-domain fluorescence lifetime images[END_REF].

B.2. Lifetime estimation on a control experimental sample

Figure 9 presents the lifetime distribution estimated on an epithelial cell membrane (RPE1) using the Fourier decomposition and our weighted least squares method on a control sample. Experimental data are a subset of the acquisition carried out in the context of the experiment described Section D. In this experiment, fluorescence from GFP-tagged receptors localized on plasma membrane and early endosomes was acquired in the absence of an acceptor. Experiments were conducted using a confocal microscope with spinning disk setup and the same intensifier as above.

Our method estimates a narrower distribution, highlighting a standard deviation 0.27 ns while the Fourier decomposition estimates a standard deviation of 0.48 ns (means are 2.11 ns for Fourier vs 2.15 ns, modes are 2.08 ns and 2.09 ns). It is hard to assess the quality of a lifetime distribution in experimental conditions. However, the sample under study cannot present any transfer of energy and should thus present a single lifetime distribution. It follows that a narrow lifetime distribution with a single mode is expected. One can note the significant difference in the impact of our noise model on the phase of reference sequence and the lifetime estimate on the living cells. Indeed, two phase estimates are required to obtain the lifetime map, one applied on the reference sequence and an other for the sample at hand. This accentuates the effect on distribution kurtosis.

B.3. Quantification of phase estimation accuracy with increasing photon count

To further quantify the performances of our weighted leastsquare approach on static lifetime estimation, we have simulated plain fluorescein reference FD FLIM measurement and measured the phase estimation accuracy as a function of the simulated photon count. As shown in Figure 11, our method shows a 40% decrease in root mean square error (0.062 ns vs 0.038 ns) when compared to the Fourier (or equivalently leastsquare technique) using a sinusoidal signal presenting an offset ranging from 32 to 1000 photon counts. In those simulations, the minimum sample in each sinusoidal signals is 10% of the offset, reaching a minimum photon count of three.

C. Lifetime map reconstruction on dynamical sequences

We first illustrate our lifetime reconstruction and motion estimation on both noise-free and noisy simulations, we then highlight our results on experimental biological data.

C.1. Lifetime map reconstruction on synthetic data

In order to focus the reconstruction of artifacts due to motion only, we applied our algorithm on simulations (see Figure 12) where the lifetime estimate is only affected by the vesicle motion without additional noise. The tracking is performed without issue, while parameter estimation takes 7 to 8 iterations to converge. The precision is about 10 -6 after the first iteration on this simple example and additional iterations are not necessary. We also tested the impact of our iterative approach on more challenging data. Those new simulations are based on our intensity model presented in Section A and the noise model presented in [START_REF] Roudot | Noise modeling for intensified camera in fluorescence imaging: application to image denoising[END_REF]. Background is simulated with an experimental vesicle-free sample and cytoplasmic fluorescence is simulated in the background with constant parameter θ b . Vesicle intensities follow [START_REF] Clegg | Lifetime-resolved fluorescence imaging[END_REF] and motions are limited to a two-pixel displacement. The non-stationary Gaussian noise parameters are estimated on an experimental image series and then imposed on synthetic vesicles. Table 1 shows the convergence of lifetime and motion parameter on simulated vesicles. Convergence of tracking and lifetime errors is reached after m = 3 iterations on 30 simulated vesicles and 95 % of vesicles tracked.Similarly to our experiment on static fluorescence lifetime estimation, we measured the evolution of the root mean square error of the phase estimate on moving object under the same conditions (see Figure ). In our iterative framework, we compare the use of a weighted least square estimator against a classic Fourier estimator for the iterative phase estimation. One can observe that the impact of the noise modeling is less strong on the dynamic sub-resolved objects (only 8%). This is expected given the impact of other variable such as motion estimation and model fitting error.

C.2. Endosome tracking accuracy in experimental conditions

We tested our reconstruction algorithm on living RPE1 cells expressing GFP-tagged Human Interferon Alpha Receptor Subunit-1 (IFN-R1), acquisitions that are further described in Section D. We highlight the result of our reconstruction method on experimental samples in Figure 14. The "doppler effect" can easily be seen on lifetime map computed using the usual Fourier method, and greatly reduced using our method. However, it is very challenging to assess the accuracy of the estimated endosome fluorescence lifetime on experimental samples. Nevertheless, the lifetime parameter is closely related to the quality of the displacement estimate, and we can much easily measure the quality of motion tracking. To test our algorithm in these conditions, 30 moving spots localizing the GFPtagged receptors have been imaged with our FD FLIM setup and manually tracked using the MtrackJ software [START_REF] Smal | Quantitative comparison of spot detection methods in fluorescence microscopy[END_REF]. In order to highlight the effect of the proposed intensity modeling, we first compared the performances of our method against a frameby-frame Gaussian fitting [START_REF] Anderson | Tracking of cell surface receptors by fluorescence digital imaging microscopy using a charge-coupled device camera[END_REF]. This simple method reuses the estimated position on the previous frame as an initialization for spot fitting. In order to compare our approach with a more competitive method, we also tracked those vesicles with the "probabilistic particle tracker" described in [START_REF] Genovesio | Multiple particle tracking in 3-d+ t microscopy: method and application to the tracking of endocytosed quantum dots[END_REF] and implemented in the ICY software [START_REF] De Chaumont | Icy: an open bioimage informatics platform for extended reproducible research[END_REF]. This multiple particle tracking algorithm detects spots using an undecimated wavelet transform [START_REF] Olivo-Marin | Extraction of spots in biological images using multiscale products[END_REF] on each frame before building tracks frame-by-frame. For each frame pair, the correspondences between the two sets of detections are optimized thanks to combinatorial optimization of the set of distances between detections. The distances are determined using the parameter of a Brownian model estimated through Kalman filtering . To address the problem of intensity variations that we modeled in our method, images are normalized using the Midway equalization method described in [START_REF] Angelini | Detection of glioma evolution on longitudinal mri studies[END_REF].

The detector proposed in [START_REF] Olivo-Marin | Extraction of spots in biological images using multiscale products[END_REF] cannot adapt to local changes of SNR. As a result, the lowest SNR images that are normalized by our histogram equalization present an higher number of detected objects due to noise artifacts. However, as objects are sparsely localized in our experimental dataset, the three compared algorithms managed to track the 30 objects on the entire stack of 12 frames. The mean localization errors are represented in Figure 15 and the tracking with our method is visually illustrated Figure 16. As one could expect, the localization errors occur as the signal-to-noise ratio drops at lower intensities. Thanks to our thorough intensity modeling, our method presents a more precise spot localization throughout the whole sequence, a difference that intensifies on more noisy images.

D. Application to protein interaction analysis at the plasma membrane on endosomes

An usual approach to determine molecular nm-scale proximity in living cells is the measuring the FRET (Forster Resonance Energy Transfer) efficiency between the fluorophores linked to the two proteins of interest. Thanks to our motion compensation method, we can measure lifetime loss on dynamical structures such as trafficking endosomes as well as at the membrane. To trigger the FRET phenomenon on protein proximity, we considered IFN-aR1 GFP-tagged receptor as the donor protein and its associated kinase Tyk2 tagged with mCherry at different tagging sites as its acceptor (see Appendix 3 for an illustration of Probabilistic particle tracking [START_REF] Genovesio | Multiple particle tracking in 3-d+ t microscopy: method and application to the tracking of endocytosed quantum dots[END_REF] and our method on fluorescently tagged protein in living epithelial cells. Endosomes were tracked manually using the MtrackJ software [START_REF] Smal | Quantitative comparison of spot detection methods in fluorescence microscopy[END_REF].

their emission and excitation spectra). The sub-resolved spot that are localized outside of the cluttered recycling compartment are the endosomes of interest and the plasma membrane is described by the cell background.

To differentiate the lifetime estimated on the early endosomes from the membrane, we exploit the lifetime estimated on tracked vesicles for the former and we automatically segment the cell background for the latter (Figure 17).

We summarize the results obtained on a set of acquisition using cells expressing IFN-aR1-GFP alone or with Tyk2 fused to mCherry in C-or N-terminus in Figure 18: • The first conclusion we can draw is that lifetime is lower in the condition where the acceptor is present, an indication that FRET might occur and that IFN-aR1 and its associated kinase Tyk2 may be in close proximity (at the nm scale) in RPE1 cells. This observation is consistent with the IFN-aR1/Tyk2 basal association described in the literature [START_REF] Kumar | Basal ubiquitin-independent internalization of interferon alpha receptor is prevented by Tyk2mediated masking of a linear endocytic motif[END_REF].

• No significant difference was found between endosomes and membranes, which validates also the lifetime estimation on moving endosomes. Without our algorithm, values would be probably different due to motion-induced artifact. We can also notice that the chosen site for tagging may have an impact on FRET efficiency. Indeed, both experiments with the two different tagging positions of mCherry on Tyk2 kinase show an higher lifetime loss.

CONCLUSION

In the first part of this paper, we have presented a dedicated signal processing method to reconstruct the fluorescein lifetime map corrupted by instrumental noise and intracellular motions in FD FLIM. In the second part, we proposed a fluorescence lifetime estimation method for dynamical analysis of endosomes inside the cell. Endosomes are detected on the first frame of the FD FLIM measurement. For each detected object, the lifetime and motion parameters are jointly estimated in an iterative and alternative fashion. On experimental sequences, the performance of the motion estimation is shown to outperform tracking algorithms that does not take the specific intensity model into account. The proposed methods are fully automatic and need only minimal and non-critical parameter setting.

The application of our method on the biological experiment presented in the last section demonstrated its efficiency in living cells. Thanks to those contributions we are hoping to pave the way for fluorescence lifetime imaging microscopy of dynamical intracellular structures. In particular, our method will allow the studies of protein-protein interactions in living cells during endocytosis process or signaling events on endosomes.

Two improvements could be envisaged for future work. To further improve the noise variance analysis, the next step is the study of the spatial correlation on the ICCD output. The estimation of multiple fluorescence lifetime on a pixel basis could be also interesting in the context of fluorophore with overlapping emission spectra. FD FLIM exploiting multiple frequencies have shown to be able to resolve multiple lifetimes [START_REF] Gratton | A continuously variable frequency cross-correlation phase fluorometer with picosecond resolution[END_REF][START_REF] Squire | Multiple frequency fluorescence lifetime imaging microscopy[END_REF]. To improve the tracking capacity, more diversified intracellular structures should be considered. A few leads are already in this direction work, such as pixel soft classification between object and background.

More generally, the long-term goal is to perform time lapse FD FLIM measurement on a longer acquisition time. The motivation is to quantify protein interaction throughout for several cell cycles. If our method can already handle time lapse FD FLIM measurements, the time between each measurement must be short enough for the local motion compensation to perform properly. By spacing out FD FLIM measurements, the reduced photo-toxicity and photo-bleaching would allow much longer time-lapse acquisitions. The motion estimation between each FD FLIM measurement should thus be carried out with multiple particle tracking method. Semi-global optimization strategies to optimize the set of trajectories needs then to be investigated.

MATERIAL AND METHODS

Human Retinal Pigment Epithelial (RPE1) cells were maintained in culture in DMEM/F-12 medium (Gibco) supplemented with 10% bovine fetal serum. RPE1 cells were transfected with plasmids with X-tremeGene HP DNA (Roche) according to manufacturer guidelines. DNA constructs were verified by sequencing. Human Interferon Alpha Receptor Subunit-1 (IFN-R1) cDNA was cloned into pCAGGS (kind gift of Pr. Fukata). Human Tyk2 cDNA was extracted from pRC-HuTyk2-VSV plasmid (kind gift of S. Pellegrini), and cloned into pmCherry-N1 or -C1 (Clontech).

The authors acknowledge PICT-BiSA UMR 144 CNRS for technical support. and a simple sinusoidal excitation of the following form:

E(t) = C E + A E sin(ω E t + φ E )], (28)
R(t) = R 0 e -t/τ , (29)
where R 0 denote the response at time t = 0. It follows that (1 + (ω E τ) 2 ) cos(atan(ω E τ)) .

F(t) = ∞ 0 E(t -t )R(t )dt = R 0 C E ∞ 0 e -t /
Finally, we have

F(t) = R 0 τ C E + A E sin(ω E t + φ E -atan(ω E τ)) cos(atan(ω E τ))(1 + (ω E τ) 2 ) = C F + A F sin(ω E t + φ E -atan(ω E τ)))
where C F = C E R 0 τ (see [START_REF] Holub | Fluorescence lifetime imaging at video rate a new technique in photosynthesis research[END_REF]) and

A F = A E 1 cos(atan(ω E τ))(1 + (ω E τ) 2 ) = A E 1 + (ω E τ) 2 (30)
This result shows that the fluorescence response of the sample induces a phase delay and an amplitude decrease in the excited signal which is coherent with physical intuition.

A.2. Lifetime estimation

As we control the frequency of the emission signal, [START_REF] Angelini | Detection of glioma evolution on longitudinal mri studies[END_REF] shows that the phase and amplitude of the signal yield the lifetime τ (using a calibration reference measurement to retrieve φ E ). Actually, the experimental frequency in such experiment is too high for the phase to be directly measurable on a CCD sensor. To overcome this problem, it is usual to modulate the detection sensitivity with an other sinusoidal signal M(t) (see (44)):

M(t) = C M + A M sin(ωt + φ M )]. (31
) In the homodyne case, we use the same frequency as the original signal: ω E = ω M . The result of this modulation then yields to:

M(t)F(t) = C M C F + A M A F sin(ω E t + φ M ) sin(ω E t + φ E -atan(ω E τ)) + A F C M sin(ω E t + φ E -atan(ω E τ)) + A M C F sin(ω E t + φ M )] = C M C F + A F A M cos(φ E -atan(ω E τ) -φ M ) + A M A F cos(2ω E t + φ E -atan(ω E τ) + φ M ) + A F C M sin(ω E t + φ E -atan(ω E τ)) + A M C F sin(ω E t + φ M ).
An intensified CCD is then required to acquire the lowest intensity images resulting from modulation. The higher frequencies of the phase modulated-signal F(t)M(t) are attenuated by the low pass effect of the CCD detector denoted G CCD resulting in the following approximation:

S(t) = G CCD (F(t)M(t)) ≈ C F C M + A F A M cos(φ E -atan(ω E τ) -φ M ) (32)
The measured signal is then time-independent. One solution to recover the value atan(ω E τ) is to adjust φ M as follows:

φ M = φ M 0 + 2πk K , k ∈ [1, K], ( 33
)
where φ M 0 denotes the phase for k = 0. We obtain K samples (K must be chosen > 3 ) defined as:

S(k) = C F C M + A F A M cos φ M 0 + 2πk K -φ E + atan(ω E τ) . (34) = C S + A S cos 2πk K + φ τ , k ∈ [1, K], ( 35
)
where

C S = C F C M , φ τ = (φ M 0 -φ E + atan(ω E τ)), ( 36
)
and

A S = A F A M = A M A E 1 + (ω E τ) 2 . ( 37
)
The observations S(k) can be then fitted with a cosine function (using various methods including Fourier decomposition and robust regression methods) to estimate the parameter controlling (34), most notably the phase φ τ and amplitude A S . We can thus estimate the fluorescence lifetime, provided that we can calibrate the optical setup to measure φ M 0 -φ E or alternatively A E A M . To do so, a FD FLIM measurement of reference sample is with known lifetime τ re f is acquired in the same condition. By fitting the S re f (k) with a cosine function, we obtain an estimate for the reference phase φ τ re f yielding :

(φ M 0 -φ E ) = φ τ re fatan(ω E τ re f ) (38)

A E A M = A F A Mre f 1 + (ω E τ re f ) 2 (39)
Finally, estimating the phase φ τ of the signal described by 34 yields an estimate of the fluorescence lifetime:

τ φ = tan( φ -φ τ re f + atan(ω E τ re f )) ω E , (40)
and denoted τ throughout the paper. Moreover, the estimation of the amplitude A S combined with equation 47 yield an other approach to estimate the fluorescence lifetime using:

τ M = 1 ω E ⎛ ⎝ A E A M A F A Mre f 1 + (ω E τ re f ) 2 ⎞ ⎠ 2 -1. (41)

Fig. 1 .

 1 Fig. 1. Instrumental setup showing the system components and the corresponding FD FLIM signals. We denote E(t) the excitation light, F(t) the fluorescence signal, G k (t) the detection modulation signal and S θ the FD FLIM measurement.

Fig. 2 .

 2 Fig. 2. Principle of FD FLIM measurement and post-processing A: Example of an image sequence showing the intensity variation (S θ (k)) at different phase shift values (K=12 values distributed on 2π radians). B: Example of S θ (k) observed at a given pixel (black dots) with a sine fitting curve (red). C: Fluorescence lifetime map reconstructed from the image sequence in B. Scale bar is 10 μm,

Fig. 3 .Fig. 4 .

 34 Fig. 3. Residuals of least-squares estimates of the parameters controlling the sinusoidal described by the phase modulation images (see (4)). The residuals associated to each pixel in the whole stack presented in Figure 2 are plotted. The variability of the residual amplitude highlights the presence of strong heteroscedasticity.

Fig. 5 .

 5 Fig. 5. Visualization of sub-cellular structure motion during image acquisition and effect on the FLIM images. A: Example of image sequence with moving sub-structures (i.e. series of zoomed thumbnails) and of modulation signal S θ with sine fitting in red. B. Global estimation of motion by maximum intensity projections (MIP) and ratio with initial intensity R 0 (right). C: Reconstructed lifetime image showing outliers in red. D. Distribution of lifetime values in the reconstructed lifetime map in C showing extreme, non physical (negative or too high) values. Scale bar is 10 μm.

Fig. 7 .

 7 Fig. 7. Moving vesicle model.

Fig. 8 .Algorithm 1 . 2 . 2 .i

 8122 Fig. 8. Iterative and alternative procedure for fluorescence lifetime estimation on moving endosomes.

Fig. 12 .

 12 Fig. 12. Noise-free simulation of motion impact on lifetime map and reconstruction. a) Simulated sequence with 5 vesicles exhibiting free Brownian motion modeled by a normal law for displacement with a standard deviation of 2 pixels (contrast has been enhanced for visualization). b) The estimated lifetime map on this FD FLIM measurement with the "doppler effect" footprint. c) Patches taking the value of the estimated fluorescence lifetime are added on the initial vesicle detection locus. The patch size is set to three times the size of the estimated spot scale d) A simple interpolation algorithm is applied to correct secondary artifacts for cosmetic purpose.e) True lifetime simulated on background and vesicles f) Motion estimation carried out during the iterative process. Scale bar is 10 μm.

Fig. 9 .

 9 Fig. 9. Phase histograms using the Fourier method and the weighted least squares method exploiting the estimated noise scale on a reference FD FLIM measurement of a plain fluorescein sample.

Fig. 10 .

 10 Fig. 10. Lifetime histograms on live cell RPE1 with GFPtagged receptor.

Fig. 11 .

 11 Fig. 11. Phase estimation error with respect to the minimum photon count on static FD FLIM measurement (4.10 4 FD FLIM measurements have been simulated for each photon count).

Fig. 13 .

 13 Fig. 13. Phase estimation error with respect to the minimum photon count on dynamic FD FLIM measurement (100 objects have been simulated for each photon count).

Fig. 14 .

 14 Fig. 14. Lifetime reconstruction on experimental data. Left: Lifetime map using Fourier decomposition. Right: lifetime map reconstruction using our method. Scale bar is 10 μm.

Fig. 15 .

 15 Fig.15. Mean localization error using Gaussian fitting[START_REF] Anderson | Tracking of cell surface receptors by fluorescence digital imaging microscopy using a charge-coupled device camera[END_REF], Probabilistic particle tracking[START_REF] Genovesio | Multiple particle tracking in 3-d+ t microscopy: method and application to the tracking of endocytosed quantum dots[END_REF] and our method on fluorescently tagged protein in living epithelial cells. Endosomes were tracked manually using the MtrackJ software[START_REF] Smal | Quantitative comparison of spot detection methods in fluorescence microscopy[END_REF].

Fig. 16 .

 16 Fig. 16. An example of IFN-aR1-GFP tracking on a 12 frames FD FLIM acquisition (images are normalized for visualization). Scale bar is 10 μm.

Fig. 17 .Fig. 18 .

 1718 Fig. 17. Tracked IFN-aR1-GFP on endosomes (left) and membrane lifetime map (right, from 0 to 5ns, 0 is assigned to background and endosomes) in the absence (up) or presence (down) of acceptor. Scale bar is 10 μm.

Table 1 .

 1 Tracking and lifetime errors computed on simulated noisy data.

	Iteration m	0	1	2	3
	Tracking mean error (pixel)	0.240 0.121 0.081 0.080
		(12%)	(6%)	(4%)	(4%)
	Lifetime mean error (ns)	0.387 0.248 0.247 0.247
		(15%)	(9%)	(9%)	(9%)

  τ dt + A E E (tt ) + φ E )e -t /τ dt = R 0 (C E τ + A E B(t)) E (tt ) + φ E )e -t /τ dt = -τ sin(ω E (tt ) + φ E )e -t /τ ∞ E (tt ) + φ E )e -t /τ dt = τ sin(ω E t + φ E ) --ω E τ 2 cos(ω E (tt ) + φ E )e -t /τ ) E (tt ) + φ E )e -t /τ dt = τ sin(ω E t + φ E ) -ω E τ 2 cos(ω E t + φ E ) -(ω E τ) 2 B(t). E τ) 2 (sin(ω E t + φ E ) -ω E τ cos(ω E τ + φ E )) = τ 1 + (ω E τ) 2 sin(ω E t + φ E ) -sin(atan(ω E τ)) cos(ω E τ + φ E )

					∞
	0 sin(ω where
	B(t) =	0	∞	sin(ω 0
					∞
			-ω E τ	0	cos(ω ∞
					0
	∞ sin(ω It follows that +(ω E τ) 2 0
	B(t) =	τ 1 + (ω

cos(atan(ω E τ)) = τ cos(atan(ω E τ)) sin(ω E t + φ E ) (1 + (ω E τ) 2 ) cos(atan(ω E τ)) -τ sin(atan(ω E τ)) cos(ω E τ + φ E ) (1 + (ω E τ) 2 ) cos(atan(ω E τ)) = τ sin(ω E t + φ Eatan(ω E τ))

ACKNOWLEDGMENTS

French National Research Agency (ANR-10-INBS-04-07, "Investments for the future").

APPENDIX A: DERIVATION OF THE PHASE MODU-LATED SIGNAL FOR FD FLIM

This appendix presents the derivation of the model used to recover the fluorescence lifetime from an excited sample. Mono-exponential decay and simple sinusoidal excitation is assumed. We start by deriving the equation of the excited signal, then we explain in details the experimental method to estimate the fluorescence lifetime of a given sample. To our knowledge, this derivation has not been described to such a full extent before in the literature.

A.1. The excited signal

Let E(t) and R(t) be an emission signal and the response of a fluorescent molecule respectively at time t and let τ be the unknown lifetime. The excited fluorescent signal is given by F(t) = E(t) * R(t) where * denotes the convolution operator. Assume a mono-exponential decay

APPENDIX B: NOISE VARIANCE MODELING

In this Appendix we summarize the derivation of a theoretical and experimental model for the variance of the noise induced by the ICCD setup previously published in [START_REF] Roudot | Noise modeling for intensified camera in fluorescence imaging: application to image denoising[END_REF]. Considering the gain and read out noise, Boddeke showed that the ICCD response I(x, k) at location x = (x, y) ∈ Ω is [START_REF] Boddeke | Quantitative fluorescence microscopy: autofocusing, z-axis calibration, image sensors, fluorescence lifetime imaging[END_REF]:

where ℵ(x) is the incident photon number on the ICCD which follows a Poisson law of parameter λ(x), ξ(x) is the CCD read out Gaussian noise such as ξ(x) ∼ (m ξ (x), σ 2 ξ (x)), g CCD is the gain of the CCD sensor and g INT is the gain of the intensifier assumed to be a realization of a random variable. Under the proper assumptions of independence of random variables it comes:

By taking the expectation of (42) and since Var[ℵ(x)] = λ(x) = E[ℵ(x)], it follows that:

The local variance thus follows a quadratic relationship with the local expectation that we summarize with:

with a,b and c ∈ R. We have shown experimentally in [START_REF] Roudot | Noise modeling for intensified camera in fluorescence imaging: application to image denoising[END_REF] that this model holds only locally. Indeed, when measuring the noise variance over the full spatial domain Ω, we notice a stronpg spatial dependency of the parameters (a, b, c). Further experiments in [START_REF] Roudot | Noise modeling for intensified camera in fluorescence imaging: application to image denoising[END_REF] have shown that the spatial dependency follows a bell-shape curve due to optical aberration induced by the lens in-between the phosphore screen and the CCD sensor in the intensifier. To correct for this aberration we have proposed the following model: where x 0 , y 0 , σ x , σ y and o are the center, scale and offset of the correction term. In [START_REF] Roudot | Noise modeling for intensified camera in fluorescence imaging: application to image denoising[END_REF], we propose a two-step approach for the estimation of the eight parameter controlling (47). Despite the simplicity of our correction model, the magnitude of the residuals of our estimation process indicate that our model is suitable to take into account the lens-induced aberration and fully explain the measured noise footprint.

APPENDIX C: EXCITATION AND EMISSION SPECTRA OF DONOR AND ACCEPTOR


Nicolas Bertheuil

Sophie Thienot

Vincent Huguier

Cédric Ménard

Eric Watier

Medial Thighplasty After Massive Weight Loss: Are There Any Risk Factors for Postoperative Complications?

Background

With the increasing incidence of obesity in western societies, more and more patients undergo bariatric surgery that leads to functional and aesthetic sequelae related to a rapid and massive weight loss, namely skin and fat excess. The goal of surgical management after massive weight loss is to optimize the functional results obtained from bariatric surgery or diet with removal of the redundant skin folds. We report our experience with medial thighplasty following massive weight loss and identify predictive factors of postoperative complications.

Methods

This is a retrospective study of 53 patients operated for medial thighplasty after massive weight loss. Data were collected through chart review and Submitted to statistical analysis.

Results

The average weight loss before thighplasty was 43.5(±14.9) kg, 39.6% of patients experienced at least 1 complication. Complication rates were as follows: seroma 9.4%, wound dehiscence 20.8%, scar migration 17%, wound infection 5.7% and partial skin necrosis 1.9%.

Pre-massive weight loss body mass index (BMI) and pre-medial thighplasty BMI were found to be risk factors for developing a postoperative complication.

Conclusion

Few guidelines exist on the optimization of care and safety in this complex patient population.

The informations provided herein by this retrospective study complement the current data of the literature and can help surgeons to select patients eligible for medial thighplasty. We show that pre-medial thighplasty BMI should be taken into consideration for patients willing to go on body contouring procedure.

Evidence Based Medicine Level III: Evidence obtained from well-designed cohort or casecontrol analytic studies, preferably from more than one center or research group.

Introduction

Obesity is a major health risk factor because of associated medical complications such as cardiovascular disease, hyperlipidemia, diabetes, high blood pressure, cancer, degenerative arthritis and sleep apnea [START_REF] Gallagher | Insulin resistance in obesity as the underlying cause for the metabolic syndrome[END_REF][START_REF] Lavie | Obesity and cardiovascular disease: risk factor, paradox, and impact of weight loss[END_REF][START_REF] Van Gaal | Mechanisms linking obesity with cardiovascular disease[END_REF][START_REF] Montague | The perils of portliness: causes and consequences of visceral adiposity[END_REF]. Obesity is a multifactorial disease whose one of the components seems to be from genetic origin. Several studies have identified genetic polymorphisms associated with obesity [START_REF] Bochukova | Large, rare chromosomal deletions associated with severe early-onset obesity[END_REF][START_REF] Clement | A mutation in the human leptin receptor gene causes obesity and pituitary dysfunction[END_REF][START_REF] Emilsson | Genetics of gene expression and its effect on disease[END_REF][START_REF] Menzaghi | A haplotype at the adiponectin locus is associated with obesity and other features of the insulin resistance syndrome[END_REF]. Others have been involved in the progression of obesity such as environmental and behavioural factors and socio-economic status [START_REF] Strawford | Adipose tissue triglyceride turnover, de novo lipogenesis, and cell proliferation in humans measured with 2H2O[END_REF].

Obesity may lead to difficulties in daily life (clothing, walking, sexual activity, hygiene problems), altering the quality of life and generating psychological disorders such as devalued self-image and depression. Bariatric surgery was developed to combat against obesity and is indicated for patients with a body mass index (BMI) ≥ 40kg/m 2 or BMI ≥ 35kg/m 2 with an additional related disease and an acceptable surgical risk [START_REF]Gastrointestinal surgery for severe obesity: National Institutes of Health Consensus Development Conference Statement[END_REF]. With the increasing incidence of obesity in western societies, more and more patients undergo this type of surgery that leads to functional and aesthetic sequelae related to a rapid and massive weight loss, namely skin excess and poor skin tone. The goal of surgical management after massive weight loss is to optimize the functional outcomes obtained from bariatric surgery or diet with surgical removal of the redundant skin folds.

Within this field, the medial thighplasty is now well acknowledged. This kind of surgery was first described by Lewis et al in 1957 [START_REF] Lewis | The thigh lift The[END_REF][START_REF] Lewis | Correction of ptosis of the thighs: the thigh lift[END_REF] and allows to reduce redundant skin and fat inside the thighs. Minor complications often occur such as scar enlargement, scar migration, recurrence of ptosis and vulvar problems [START_REF] Louarn | Internal faces lifting legs complications[END_REF]. However, improvements were made to the initial technique [START_REF] Schultz | Medial thigh lift[END_REF][START_REF] Vilain | Aesthetic surgery of the medial thigh[END_REF][START_REF] Louarn | The concentric medial thigh lift[END_REF] by anchoring the superficial fascial layer to the Colles' fascia as described by Lockwood et al which has allowed to significantly reduce the number of cases of ptosis [START_REF] Lockwood | Fascial anchoring technique in medial thigh lifts[END_REF].

Two types of incision may be done depending on the location of the skin fold [START_REF] Mathes | Current concepts in medial thighplasty[END_REF]: a horizontal approach for the excess on the upper third of the thigh or a vertical approach in case of important excess particularly after massive weight loss [START_REF] Labardi | Medial thighplasty: horizontal and vertical procedures after massive weight loss[END_REF][START_REF] Kenkel | Medial thigh lift[END_REF]. Over the past decade, liposuction has been more and more often combined to medial thighplasty in order to eliminate the fat excess to such and extent that it is now often the first step of the thigh lift.

Concomitant liposuction has the advantage of improving the aesthetic results of the intervention, and reduces the healthcare cost by combining two operations into one. But this technique is regarded by some as risky as it may lead to ischemia and necrosis of the flaps.

There are presently few reports on the postoperative complications of medial thighplasty after massive weight loss. It seems necessary to improve the knowledge about this technique in order to have a better management of patients.

The aim of our retrospective study, including 53 patients, was to report our experience with medial thighplasty following massive weight loss and identify factors predictive of postoperative complications.

Materials and methods

We report the experience of the Department of Plastic Surgery of Rennes University Hospital through a retrospective clinical study of 53 patients operated for medial thighplasty consecutively between September 2004 and December 2011. All of the patients who underwent a medial thighplasty were included. Data collection was done through analysis of medical records for each patient, including age, sex, weight and height, pre-bariatric surgery body mass index (BMI), pre-medial thighplasty BMI and weight loss before thighplasty. History of bariatric surgery was sought as well as diabetes, high blood pressure and osteoarthritis. Data regarding surgery included the type of medial thighplasty (horizontal or vertical incision), the use of liposuction, the pannus weight, operating duration, postoperative drainage duration, postoperative complications and duration of hospital stay. All operated patients had their weight stabilized at a plateau over a period of at least 6 months. Active smokers were asked to stop smoking for 1 month before the reconstructive procedure and for 1 month postoperatively.

This study has been performed in accordance with the principles of the Declaration of Helsinki (1964) and was in agreement with French bioethics laws of 7 july 2011.

Operative procedure

Drawings were made preoperatively in the standing and prone position with lower extremities in abduction. The superior incision is marked in the perineal-thigh fold. The anterior extent of the incision was stopped just before to the femoral triangle. The posterior extent continues into the medial buttocks fold. The location of excessive skin and fat is determined that allows to choose between horizontal scar (n=50) or a combinaison of horizontal and vertical scar (n=3). Then, we identify the point of maximum resection by pinching the skin. This point was marked. Then we drew an ellipse between this point and the anterior and posterior end of the drawing. Finally, suction areas are identified and marked.

Surgery was performed under general anesthesia and the patient was set up in the supine position combined with gynecological position. We performed an infiltration by lidocaineadrenaline to reduce bleeding during surgery and to facilitate liposuction. Liposuction, was always performed as the first step of surgical procedure, with the use of 3 or 4 mm-diameter cannulas. Liposuction was performed all along the thigh under the superficial fascial layer, and in addition under the skin pattern where all the fat under the dermis were eliminated.

Then, the incision was lead until the Colles' fascia. We realized the removal of excess skin "on demand" as previously described [START_REF] Louarn | The concentric medial thigh lift[END_REF]. The anchoring was performed with non-absorbable stitches beetween superficial fascial layer and Colles's fascia, then subcutaneous/sub-dermal planes were sutured with absorbable stitches (Monocryl® 3.0 and 4.0). In all cases, suction drains were placed and skin closure was subsequently achieved without tension. An adhesive occlusive dressing (OpSite®) was placed to prevent any contamination. Finally, compressive clothing (Lipo-Panty® standard S001 Medical Z) had to be worn for one month postoperatively.

Postoperative care

In all patients, phlebitis was prevented by the subcutaneous administration of low molecular weight heparin (Lovenox®) for a period of 10 days postoperatively and wearing compressive stockings for a period of 15 days. Antibiotics were given only in case of signs of local infection.

Statistical analysis

Statistical analysis was performed using Prism 5 GraphPad Software. A descriptive analysis of all data was first carried out. Then a comparative study was conducted between the patients who presented with a complication and those who did not. Categorical variables were compared with chi-square analysis or Fisher exact test. Quantitative variables were compared by Mann Whitney U non parametric test. Differences between groups were considered significant if p<0.05.

Results

Our series includes 53 patients who underwent surgery in order to improve aesthetic appearance and functional outcome of their thighs and who received a surgical treatment by medial thighplasty. Patient characteristics are summarized in Table 1.

The mean age was 47.2 (±11.4) years. The mean body mass index was 29 (±4.8) with an average weight loss before operation of 43.5 kg (±14.9, ranging from 20 to 84 kg). Fortyseven point two percent of our patients had a history of bariatric surgery either through laparoscopic adjustable gastric banding or by Roux-en-Y gastric bypass. Liposuction was performed in 71.7% of patients and the horizontal technique was used in 94.3% of patients.

The mean weight of bilateral pannus was 760 (±455) grams. The mean drainage duration was 5.6 (±2) days. The mean follow up was 51.1 months (±21.6 months), ranging from 12 to 112 months. The mean length of stay was 6 +/-1.8 days.

In our series, 39.6% of patients experienced at least 1 complication. Early complications were seroma (5 cases, 9.4%), posterior wound dehiscence (11 cases, 20.8%) treated through healing by secondary intention, wound infection (3 cases, 5.7%) and 1 case of partial skin necrosis (1.9%). Late complication was the scar migration (9 cases, 17%). There were no cases of large necrosis, lymphoedema, deep vein thrombosis, pulmonary embolism, nor any vulvar problems.

Comparing patients who developed a postoperative complication with those who did not, there were few statistically significant variables. Only a high massive weight loss BMI (Figure 1) and a high pre-medial thighplasty BMI (Figure 2) were found to be statistically associated with an increased postoperative complication rate (p = 0.042 and 0.044 respectively). Likewise, a low complication rate was associated with a low BMI (p = 0.048, Chi-square analysis, Figure 3). No difference was found regarding at in the age, the weight, the weight loss before thighplasty, BMI decrease, the history of bariatric surgery, diabetes, blood pressure, the use of liposuction, the operating duration and the vertical or horizontal approach.

Discussion

Obesity is a major public health concern in western societies. Beside medical treatments, bariatric surgery was developed for patients with class III obesity (Body Mass Index > 40 kg/m 2 ). Presently, it is the best treatment for these patients to reduce mortality and morbidity [START_REF] Adams | Long-term mortality after gastric bypass surgery[END_REF], particularly type 2 diabetes [START_REF] Buchwald | Weight and type 2 diabetes after bariatric surgery: systematic review and meta-analysis[END_REF], which is corrected even before weight loss [START_REF] Rubino | The early effect of the Roux-en-Y gastric bypass on hormones involved in body weight regulation and glucose metabolism[END_REF]. With the spread of obesity worldwide, the number of patients operated by this kind of surgery is steadily increasing.

The effects of massive weight loss after this kind of surgery often require multiple operations to correct the redundant skin and fat excess, such as: breast reduction mammaplasty, abdominoplasty, medial thighplasty and brachioplasty. Overall management of these patients is difficult because their obesity and the consequences of bariatric surgery may generate significant psychological disorders. These disorders may explain an increased suicide rate in this population versus obese patient without bariatric surgery (2.6 vs 0.9 suicide/10 000 patients/years) [START_REF] Adams | Long-term mortality after gastric bypass surgery[END_REF], approximately 70% of these suicide attempts occuring in the first three years after operation [START_REF] Tindle | Risk of suicide after long-term follow-up from bariatric surgery[END_REF]. It is likely that patients are not adequately preoperatively informed of the sequelae of massive weight loss and of existing solutions to treat these sequelae.

Clearly, skin laxity of the medial thigh is both an aesthetic and a functional problem possibly causing skin maceration and interfering with sexual activity, clothing and walking. Medial thighplasty has been established for several decades since it was first described by Lewis [START_REF] Lewis | The thigh lift The[END_REF].

Anchoring techniques are now well described [START_REF] Lockwood | Fascial anchoring technique in medial thigh lifts[END_REF] which help to stabilize the results of surgery over time. In addition, it is now accepted for liposuction to be integrated to the thigh lift technique. However, no study compared in an objective manner the benefits that could provide liposuction within the same operative procedure.

We herein report one of the largest series in the literature of medial thighplasty. In our study, postoperative complication rate is comparable to those found in the literature. In our series, only pre-bariatric BMI and pre-medial thighplasty BMI were statistically found to be associated with an increased postoperative complication rate. Our analysis does not take into account the influence of one variable on another. This is to take into account in the interpretation of results. Pre-medial thighplasty BMI is for us the best predictive factor for postoperative complication. Thus, we propose this parameter to be taken into account in the decision to operate the patient, contrariwise to the pre-bariatric BMI that is not the current ratio on the day of consultation. Currently, this study is the only one that investigated the risk factors for complications after medial thighplasty. Our results are consistent with Arthurs et al [START_REF] Arthurs | Post-bariatric panniculectomy: prepanniculectomy body mass index impacts the complication profile[END_REF] who found the same predictive risk factor on a retrospective series of 126 abdominal dermolipectomy.

In our study, patients with a BMI of less than 25 kg/m 2 had a null complication rate while it increased up to 57% for BMI ≥ 35 (Figure 3). This confirms the importance of BMI in the selection of patients eligible for thighplasty, but there is currently no limiting value of BMI above which it would seem unreasonable to go on surgery. Colwell et al [START_REF] Colwell | Optimization of patient safety in postbariatric body contouring: a current review[END_REF] in their review of the literature suggests that patients with a BMI greater than 32 should not be operated, but their assumption does not rely on a high level of evidence. Patients should be informed that the higher the BMI, the greater the risk of complications is important. In this case they have to accept a higher risk of complications during the surgical procedure. In our opinion, based on the literature, the most important point is to get a stable weight on a period of at least 6 months before thighplasty and a stop of the tobacco consumption [START_REF] Colwell | Optimization of patient safety in postbariatric body contouring: a current review[END_REF][START_REF] Van Der Beek | Complications after body contouring surgery in post-bariatric patients: the importance of a stable weight close to normal[END_REF].

It is interesting to note that the choice of bariatric surgery does not affect the rate of complications. Indeed, there is not more complications in patients who lost weight by dieting than in those who underwent bariatric surgery. On the other hand, neither the surgical approach nor the use of liposuction does affect the results. Unlike Duff et al [START_REF] Duff | Fleur-de-Lys abdominoplasty--a consecutive case series[END_REF] who found the weight of resected tissue to impact complications rate in their series of 68 patients who underwent an abdominal dermolipectomy, this is not the case in our series of 53 patients for other pannus location.

There are mainly two schools on why and when liposuction should be used in massive weight loss population [START_REF] Shrivastava | Body contouring surgery in a massive weight loss patient: an overview[END_REF]. Like us, some teams perform concomitant liposuction. This step should be done with careful and selective undermining of the flap. Indeed, it allows for a better mobility of the flap and achievement of skin closure without tension. Liposuction can be safely used as examplified in our series and others [START_REF] Gusenoff | Brachioplasty and concomitant procedures after massive weight loss: a statistical analysis from a prospective registry[END_REF][START_REF] Kolker | The circumferential thigh lift and vertical extension circumferential thigh lift: maximizing aesthetics and safety in lower extremity contouring[END_REF][START_REF] Kolker | Maximizing aesthetics and safety in circumferentialincision lower body lift with selective undermining and liposuction[END_REF] if patient's weight has reached a plateau for at least of 6 months before surgery and if they dropped smocking in perioperative procedure. Another strategy is to perform liposuction six months prior to pannus removal. This can decrease morbidity of the flap associated with oedema caused by liposuction but this strategy needs 2 surgical procedures and this increases the cost of the procedure.

Limitations of this study need to be kept in mind; its retrospective design subjects it to the associated errors and biases inherent in this type of study. These data should be confirmed in a prospective study on a larger patient cohort.

Conclusion

The aim of this large retrospective study of 53 patients was to report our experience and evaluate our surgery outcomes concerning medial thighplasty in post massive weight loss.

Few guidelines exist on the optimization of patient care and safety in this complex patient population. The information provided by this retrospective study complements the current data in the literature and can help surgeons to select patients eligible for surgery by medial thighplasty. We show that pre-medial thighplasty BMI should be taken into consideration for developing scoring the risk of complication in the postoperative period to select the best time for patients to go on body contouring procedure.

Figure 1

 1 Figure 1 Pre-massive weight loss BMI risk factor. The values for Body Mass Index in the Y axis are in kg/m 2 and were analysed with Mann Whitney U test. Mean is depicted by a scale bar.

Figure 2

 2 Figure 2 Pre-medial thighplasty BMI risk factor. The values for Body Mass Index in the Y axis are in kg/m 2 and were analysed with Mann Whitney U test. Mean is depicted by a scale bar.

Table 1

 1 The values for categorical variables are numbers and percentage and the values for quantitative variables are mean +/-SD.

	Tables

Table 2

 2 The values for categorical variables are numbers, statistical analysis was performed by the Fisher exact test. The values for quantitative variables are mean +/-SD and were analysed by the Mann Whitney test.

	-Statistical analysis
		complications	no complications	p
	Age	47,7 ±10,0	47,2 ±12,4	0,961
	Weight	79,1 ±9,8	73,4 ±12,6	0,066
	Pre-massive weight loss BMI	48,3 ±7,9	43,4 ±6,8	0,042
	Weight loss (Kg)	44,9 ±16,1	40,4 ±17,1	0,651
	Pre-medial thighplasty surgery BMI	30,8 ±4,5	27,9 ±4,8	0,044
	Delta BMI	17,5 ±6,2	14,8 ±6,9	0,275
	history of bariatric surgery	12	13	0,168
	history of diabetes	5	5	0,475
	history of high pressure blood	4	8	1,000
	horizontal approach	20	30	0,282
	liposuction	15	24	1,000
	weight of skin resection (grams)	844,0 ±491	657,0 ±400	0,393
	operative time (minutes)	111,5 ±22	107,3 ±27	0,492

Conflict of interest disclosure statement: The authors declare that they have no conflict of interests. Overall postoperative complication rate is expressed as a number of patients in the Y. The Body Mass Index in the X axis is in kg/m 2 and differences between groups were analysed by the Chi-square test.