Optimizing Laguerre expansion based deconvolution methods for analysing biexponential fluorescence

: Fast deconvolution is an essential step to calibrate instrument responses in big fluorescence lifetime imaging microscopy (FLIM) image analysis. This paper examined a computationally effective least squares deconvolution method based on Laguerre expansion (LSD-LE), recently developed for clinical diagnosis applications, and proposed new criteria for selecting Laguerre basis functions (LBFs) without considering the mutual orthonormalities between LBFs. Compared with the previously reported LSD-LE, the improved LSD-LE allows to use a higher laser repetition rate, reducing the acquisition time per measurement. Moreover, we extended it, for the first time, to analyze bi-exponential fluorescence decays for more general FLIM-FRET applications. The proposed method was tested on both synthesized bi-exponential and realistic FLIM data for studying the endocytosis of gold nanorods in Hek293 cells. Compared with the previously reported constrained LSD-LE, it shows promising results.


Introduction
Fluorescence lifetime image microscopy (FLIM) has been showing great potential for visualizing biomolecules or physiological parameters (such as Ca 2+ , O 2 , pH, temperature) in fixed and living cells in biology and medicine [1][2][3][4][5][6][7][8]. It can localize specific fluorophores as usual fluorescence intensity imaging, but more importantly it can probe local environments around the fluorophores. It can be implemented in scanning confocal or multi-photon microscopes, or in wide-field microscopes and endoscopes [4][5][6][7]. In this paper we will focus on time-domain time-correlated single-photon counting (TCSPC) or time-gated FLIM techniques, where the measured fluorescence response from biological tissues to a laser excited pulse is a convolution of the intrinsic fluorescence impulse response function (fIRF) (emitted from tissues) with the instrument impulse response function (iIRF) (contributed by the distorted light pulse, detection mechanisms, instrument electronics and other delay components [9]).
In FLIM experiments using Förster resonance energy transfer (FRET) techniques to study cellular protein interaction networks [10-12], fIRF can be modelled as a sum of two exponentials (or more than two if fluorophores have complicated fluorescence decay profiles). Accurately measuring lifetime components is very crucial for detecting protein-protein interactions and protein conformational changes inside living cells [13][14][15][16][17]. In this paper we will examine the applicability of the proposed approach for bi-exponential fluorescence decays, especially for those having a fast (sub-nanosecond) and a slow (nanosecond) components, common situations encountered when we use FLIM-FRET techniques to study the endocytosis of gold nanorods (GNR) in living cells. The study results will be valuable for exploring research in drug delivery or disease treatment.
Deconvolution techniques used to recover the fIRF from the fluorescence histograms measured by TCSPC systems are very important and typical in FLIM analysis. Although tailfitting is widely used for fast analysis, scientists are still keen to know exact estimations. Numerous deconvolution techniques have been proposed [18][19][20][21], and among them least squares deconvolution based on Laguerre expansion (LSD-LE) has been proven effective showing superior sensitivity in disease detection [22][23][24][25][26]. There are, however, two mysterious factors to be chosen in order to use LSD-LE properly in diagnosis or parameter identification applications [27-31]. These LSD-LE methods using ordinary least squares analysis (OLSD-LE) can easily cause 'overfitting' (i.e. fitting the noise instead of the true signals), a serious problem especially when the acquisition has to be fast for diagnosis applications where the detected photon signal is contaminated by noise. Liu et al. presented a constrained LSD-LE (CLSD-LE) to avoid this problem, and they concluded that the chosen Laguerre basis functions (LBFs) should be mutually orthonormal [9] within the observation window (T) to perform CLSD-LE. To meet this condition, however, T needs to be much larger than the slowest decay to ensure that the LBFs, fIRF, and the derivatives of LBFs are 'sufficiently close to zero' at t ~T (t being the time delay with respect to the laser pulse). This would require shinning pulsed lasers with a low duty-cycle, reducing the efficiency of photon collection.
There are two major contributions in this paper. First, we introduce new criteria for selecting LBFs according to the residual level of Laguerre expansion instead of the mutual orthonormalities between LBFs. This will allows LSD-LE to be applicable even when T is comparable to the largest lifetime component. Second, the selection criteria are expanded to study bi-exponential decays for more general FLIM-FRET applications instead of only for diagnosis applications where single-exponential approximations might not produce enough contrast. To demonstrate the performances, we will test the proposed approaches on both synthesized data and realistic FLIM data.

Theory and method
In a time-domained FLIM experiment the measured fluorescence decay y(t) is the convolution of the fIRF h(t) and iIRF I(t): where (t) where b l (k; ) is determined by the Laguerre dimension L and the scale , and c l is the l th expansion coefficient. It is well known that LBFs form an orthonormal basis set only when N→∞. Previous studies concluded that the parameters L and should be chosen such that the LBFs and their corresponding derivatives should be 'sufficiently close to zero' at t ~T [9]. This condition would need a much larger T (compared to the largest lifetime component; 2 in this case) by using a pulsed laser with a lower duty cycle. In fact the expansion of fIRF with LBFs is simply a fitting problem where the optimal criteria should be the extent to which the sum of squared errors (SSE) can be minimized regardless of the orthonormalities between the LBFs used within the observation window 0 ≤ t ≤ T. Here, we define the normalized SSE (NSSE) for the fitting as: Minimizing NSSE h would be a straightforward way to assess the performances for fitting an fIRF with different lifetimes by Laguerre expansion. The Laguerre scale determines the rate of the exponential asymptotic decline of LBFs. The fIRF with a small lifetime prefers a small , whereas the one with a larger lifetime prefers a larger . When the field of view contains fluorescence decays with a wide range of lifetimes, a strategy to find the optimized should be in place. To facilitate the discussions, we rewrite the LSD-LE equations here. Substituting Eq. (4) into Eq. (2), we can obtain The deconvolution is to estimate c l , and Eq. (6) can be rewritten as . =+ yV c (7) With linear optimization, we have Equation (8) is called OLSD-LE. The main problem of OLSD-LE is that it causes overfitting easily when a higher L is applied [9]. To avoid this a smaller L was often suggested, however a smaller L is not able to resolve small lifetimes and therefore lowers the lifetime dynamic range. Liu et al. introduced a constrained LSD-LE called CLSD-LE to overcome this overfitting problem, and a higher L can then be applied [9]. Here only the conclusion is given, where ˆ is the solution to non-negative least square problem The deconvolutions and the computations of f D and 1, 2 τ are all curve-fitting problems. As Eq.
(5), to assess their performances, the NSSE of (k), D (k) and E (k) are defined as  Figure 1 shows the flow diagram summarizing how the performances of OLSD-LE and CLSD-LE on bi-exponential decays (0 < f D < 1, 0.1ns ≤ 1 ≤ 0.9ns and 2ns ≤ 2 ≤ 3ns) were assessed in four steps, highlighted in different colours. The flow diagram can also be applicable to realistic FLIM data by replacing the synthesized data with the measured data.

Symthesized FLIM data analysis
Firstly, LBFs were chosen to ensure that a given NSSE h , Eq. To reduce instability (ill-conditioned V T V deteriorates stability), computation and overfitting (for OLSD-LE), the smallest L within the shadowed area is suggested. Therefore, Fig. 2(b) shows that the optimal LBFs for OLSD-LE has L = 12 and = 0.924, whereas Fig. 3(b) suggests L = 16 and = 0.912 for CLSD-LE. Obviously CLSD-LE needs a larger L and is slightly more complicated computationally.  In general, LSD-LE methods need to specify L and properly for robust analysis. Usually the lifetime range in the field of view can be known before experiments; this information can be used to obtain new residual plots similar to Figs. 2 and 3. For a given residual requirement, L is suggested to be as small as possible to ensure a faster analysis speed. On the other hand, to improve the resolvability for the smaller lifetime 1 , it is required that L cannot be too small. The selection of L is a trade-off between the speed and the lifetime resolvability, whereas determines the accuracy of the fitting. For these reasons L = 16 and = 0.912 and L = 12 and = 0.924 are chosen for CLSD-LE and OLSD-LE, respectively.
Secondly, the synthesized decays, y(k), were generated according to the paramters listed in Table 1. There are nine different h(k), with (f D , 1 ) = (0.8, 0.2ns), (0.8, 0.5ns), (0.8, 0.8ns), …, and (0.2, 0.8ns) respectively, generated from Eq. (3), and nine corresponding y(k) (k = 1,2,…,9) were generated from Eq. (2). Thirdly, CLSD-LE (L = 16, = 0.912), OLSD-LE (L = 12, = 0.924 and L = 16, = 0.912) were applied to y(k) to obtain the recovered fIRF, D (k), and NSSE y and NSSE hD were used to assess the performances as shown in Fig. 4 where axis x corresponds to k in y(k). Being different from the previous analysis in Fig. 3 where Poisson noise was not included. This analysis shows how a larger L is more likely to cause overfitting after Poisson noise sources are included. In this analysis, 500 Monte Carlo simulations were performed for each y(k). Because of overfitting, µ NSSE,y for OLSD-LE (L = 16) is larger than µ NSSE,y for OLSD-LE (L = 12). Although µ NSSE,y for OLSD-LE (L = 12) is almost equal to that for CLSD-LE, µ NSSE,hD of OLSD-LE is in general larger than that of CLSD-LE, showing that CLSD-LE performs better and produces a D (k) closer to h(k). Finally, LSE is applied to D (k) to obtain f D and 1, 2 τ . Figure 5 shows NSSE hE for CLSD-LE and OLSD-LE. Again, µ NSSE,hE for CLSD-LE is smaller than that for OLSD-LE. It shows that all E (k) obtained by CLSD-LE are closer to h(k), giving much better estimations.  Figure 6 shows the performances of the estimated f D and 1, 2 τ . All µ x (x = f D or 1, 2 τ ) obtained by OLSD-LE and CLSD-LE are close to x r and all Variance are bigger than the corresponding Bias 2 , suggesting that both methods are effective. And x and Variance for CLSD-LE are smaller than those for OLSD-LE, indicating that CLSD-LE is more robust against the noise. Figure 6 shows the performances of f D and 1, 2 τ : (a) and (b) for f D , (c) and (d) for 1 , and (e) and (f) for 2 . Figures 6(a), 6(c), and 6(e) show that with a fixed f D a larger 1 gives less precise estimations, and with a fixed 1 a reduced f D gives less precise 1 but more precise 2 . These trends are reasonable. Figures 6(b), 6(d), and 6(f) show that although OLSD-LE produces less biased results and each case is Variance-limited, CLSD-LE produces smaller variances and therefore performs better in all cases.  x (x = f D , 1 , or 2 ), N C is the photon count; it is used to characterize the photon efficiency of an algorithm [37]) and the bias ( x/x) using Liu's CLSD-LE and our CLSD-LE. Figures 7(b), 7(d) and 7(f) shows that our CLD-LE has comparable or better F-value performances than Liu's CLSD-LE for 2 = 2.5 (T/ 2 = 4). However, Figs. 7(a), 7(c) and 7(e) shows our CLSD-LE has superior bias performances. Liu's CLSD-LE needs to meet the requirement that the LBFs should be orthonormal and therefore the largest they can use is 0.877 when L = 16 (in order to compare with our method). A lower usually contributes a larger bias. To demonstrate how the ratio T/ 2 affects Liu's CLSD-LE, we reduced T/ 2 to 3.3 by setting 2 = 3ns. Figure 7 shows that Liu's CLSD-LE has worse bias performances in all parameters, whereas the proposed CLSD-LE has similar bias performances as the previous example (T/ 2 = 4). The F-value of Liu's method seems smaller for 1 and 2 , but this should not be misled to conclude that its photon efficiency is better [38]. Instead, it is due to the seriously biased estimations [38]. Compared with Liu's approach, the proposed CLSD-LE performs more consistently.

Real FLIM data analysis
The proposed method was also tested on two-photon FLIM images of Cy5-ssDNA-GNRs labelled Hek293 cells. The images are for evaluating the endocytosis of gold nanorods (GNR) in living cells. The detailed synthesis of GNR-based RNA nanoprobes can be found elsewhere [39]. In brief, GNRs were functionalized with thiolated oligonucleotides (ssDNA) labeled with Cy5 through ligand exchange and salting aging process. After the incubation with Cy5-ssDNA-GNRs, Hek293 cells were washed and fixed with paraformaldehyde. Two-photon FLIM experiments were performed on an LSM 510 confocal microscope (Carl Zeiss) using the SPC-830 TCSPC acquisition system (Becker & Hickl GmbH). A Ti:sapphire laser (Chameleon, Coherent) was used (at 800 nm) to generate laser pulses with a duration less than 200 fs. The timing resolution of the TCSPC is 0.039ns, and measured histograms with 256 time bins (T = 256 × 0.039 = 10ns) were recorded.   in the figure shows 1 histograms within (0, 0.5ns), and it explains why a larger L is required for resolving the lifetimes of GNRs ( 1 < 100ps). The discrepancy in 2 histograms between the proposed and the Liu's CLSD-LE is due to the fact that a large number of pixels show no energy transfer and contain a larger 2 around 3ns. For Liu's CLSD-LE, the lower T/ 2 (~3) limits its resolvability for 2 (unable to resolve 2 > 3ns) causing misinterpretation that there is energy transfer at these pixels. This observation is in good agreement with Fig. 7(e), T/ 2 = 3.3. In Fig. 8(d), there is a population of pixels showing 1 < 100ps, indicating that there is energy transfer between GNRs and Cy5. For Liu's CLSD-LE, however, a smaller L (L = 8) results in biased estimations of 1 , not able to allocate GNRs (green curve). The maximum can be applied is 0.877 (for L = 16) for Liu's CLSD-LE to meet the orthonormality requirement (note that it is only quasi-orthonormal). This lower leads to a bigger bias in 2 . For our methods, although there is a slight discrepancy between CLSD-LE and OLSD-LE, they are still be able to provide similar contrast. Compared with our previous report [37], the results show that considering the iIRF in the analysis would improve locating GNRs. The results also show that the proposed CLSD-LE and OLSD-LE produce similar results, and both work robustly even when the ratio T/ 2 is less than 4. Unlike previously reported LSD-LE [22-29] and BCMM [37] requiring a much larger T/ 2 or extra bias correction procedures (for BCMM), the proposed method can reduce the acquisition time per measurement. Figure 8(f) also shows that our CLSD-LE and OLSD-LE produce similar f D histograms, whereas for Liu's CLSD-LE a smaller T/ 2 causes biased f D estimations, see Fig. 7(a). The analysis results show that the proposed OLSD-LE and CLSD-LE are effective and have potential to be used to analyze FLIM-FRET data, with the latter showing better performances.

Conclusion
We presented new criteria to choose LBFs for LSD-LE based only on how close the Laguerre expansion can approximate the fIRF. Different from the conclusion suggested by previous studies, the proposed criteria do not need to consider the mutual orthonormalities between LBFs. The new criteria do not require that the LBFs and the corresponding derivatives to be close to zero at the end of the measurement window, and they allow using a smaller T/ 2 ratio and therefore reducing the acquisition time per measurement. We applied this upgraded method to analyzing bi-exponential decays and its performances (on both CLSD-LE and OLSD-LE) were accessed and compared against the original CLSD-LE. The results show that both the upgraded CLSD-LE and OLSD-LE can be applicable to our studies, but the former performs slightly better. Both synthesized and realistic experimental FLIM data show that the proposed CLSD-LE has better performance than the original CLSD-LE when T/ 2 is small and suggest that the proposed CLSD-LE can be an effective tool to analyze bi-exponential FLIM-FRET data. It can be further extended to study multi-exponential decays in the future. The proposed methods should be able to encourage wider applications of fast FLIM technologies and gold nanoparticles for cancer therapy [37,[39][40][41].