Lifetime Determination Algorithms for Time-Domain Fluorescence Lifetime Imaging: A Review

Abstract

Fluorescence lifetime imaging (FLIM) is powerful for monitoring cellular microenvironments, protein conformational changes, and protein interactions. It can facilitate metabolism research, drug screening, DNA sequencing, and cancer diagnosis. Lifetime determination algorithms (LDAs) adopted in FLIM analysis can influence biological interpretations and clinical diagnoses. Herein, we discuss the commonly used and advanced time-domain LDAs classified in fitting and non-fitting categories. The concept and explicit mathematical forms of LDAs are reviewed. The output lifetime parameter types are discussed, including lifetime components, average lifetimes, and graphic representation. We compare their performances, identify trends, and provide suggestions for end users in terms of multi-exponential decay unmixing ability, lifetime estimation precision, and processing speed.

1. Introduction

Fluorescence lifetime imaging (FLIM) is a vital and versatile technique for assessing molecular microenvironments of fluorophores in living cells, such as pH, O₂, viscosity, temperature, or ion concentrations [1, 2]. FLIM can be a powerful “quantum ruler” to measure subnanometer protein conformational changes and interactions by quantifying the occurrence of Förster Resonance Energy Transfer (FRET) [3, 4, 5, 6, 7]. FLIM has been used in diverse disciplines, including biology, chemistry and biophysics [8, 9, 10, 11, 12]; however, it is an indirect imaging technique that needs sophisticated data analysis to deliver meaningful information. FLIM analysis can profoundly impact the interpretations of biochemical and physical phenomena.

In time-domain approaches, FLIM usually measures a three-dimensional data cube (x-y-t), obtained with a time-correlated single-photon counting (TCSPC) system [13, 14, 15, 16], a time-gated camera [17, 18, 19, 20, 21], or a streak camera [22, 23, 24, 25, 26, 27, 28]. A time-resolved histogram h_m at (x, y) in a measured data cube can be expressed as:

$$h_m=\sum _{k=0}^m\mathit{ir}{f}_{k-m}·f_m+{ϵ}_m,$$E1

$$i{\mathit{rf}}_m={\int }_{m\mathrm{\Delta }t}^{\left(m+1\right)\mathrm{\Delta }t}\mathit{irf}\left(t\right)\mathit{dt},f_m={\int }_{m\mathrm{\Delta }t}^{\left(m+1\right)\mathrm{\Delta }t}f\left(t\right)\mathit{dt},m=0,\dots ,M-1,$$E2

where M is the number of time-bins and ∆t is the bin width.

f(t) in Eq. (1) is the underlying fluorescence decay, usually following a multi-exponential decay model,

$$f\left(t\right)=A\sum _{p=0}^{P-1}{q}_{p}exp\left(-t/{\tau }_p\right),\sum _{p=0}^{P-1}{q}_p=1,$$E3

where A is the amplitude, q_p and τ_p are the fraction and lifetime of the p^th component (p = 0, …, P-1).

irf(t) is the instrument response function (IRF), often measured using a sample with a much shorter lifetime than the width of the excitation pulse or a scattering solution. The full width at half maximum (FWHM) of IRF, ∆t_IRF, is given by [16, 29].

$$\mathrm{\Delta }{t}_{\mathrm{IRF}}^2=\mathrm{\Delta }{t}_{\text{optical}}^2+\mathrm{\Delta }{t}_{\mathrm{tts}}^2+\mathrm{\Delta }{t}_{\text{jitter}}^2,$$E4

where ∆t_optical is the optical pulse width, ∆t_tts is the detector transit time spread, and ∆t_jitter is the detection and timing electronics jitter. With h = (h₀, …, h_M-1) and irf = (irf₀, …, irf_M-1) already measured, A, q = (q₀, …, q_P-1) and τ = (τ₀, …, τ_P-1) can be extracted with a lifetime determination algorithm (LDA).

LDAs can be divided into two categories: fitting and non-fitting methods.

Fitting methods solve a nonlinear minimization problem argmin χ², where χ² is a merit function revealing the goodness of fit. χ² is determined by statistics models, such as least square estimation (LSE) or maximum likelihood estimation (MLE), using pixelwise or global fitting modes [30]. Fitting methods suffer from slow analysis due to extensive intrinsic convolutions. The Laguerre expansion method converts the nonlinear-fitting problem to a linear-fitting problem [31, 32], and speeds up deconvolution procedures.

Non-fitting methods can provide lifetime parameters much faster than fitting methods, but some can only provide average lifetimes or graphic representation rather than specific lifetime components. Non-fitting methods should be used carefully according to applications. As discussed in [33, 34, 35], the intensity-weighted average lifetime, ${\tau }_I={\sum }_{p=0}^{P-1}{q}_p{\tau }_p^2/{\sum }_{p=0}^{P-1}{q}_p{\tau }_p$, can estimate the average collisional constant k_q from the Stern–Volmer constant K_D, whereas the amplitude-weighted average lifetime, ${\tau }_A={\sum }_{p=0}^{P-1}{q}_p{\tau }_p$, can estimate FRET efficiency and assess dynamic quenching behaviors, as described by the Stern–Volmer equation. ${\tau }_I$ was previously misused in FRET efficiency estimation, introducing significant bias when the decay follows a multi-exponential decay model.

The practical use of diverse LDAs depends on applications. Due to the time-consuming estimation procedure, fitting methods are generally used for offline analysis, providing more decay information, including lifetime components and average lifetimes, than most non-fitting methods. Meanwhile, non-fitting methods are suitable for real-time FLIM applications, as they are much faster and more hardware-friendly.

This work attempts to review widely used and newly developed cutting-edge LDAs for time-domain FLIM analysis. The subsequent sections are arranged as follows. Section 2 reviews fitting methods, such as the least squares estimation, maximum likelihood estimation, global fitting using iterative convolution and variable projection approaches, and Laguerre expansion deconvolution methods. Section 3 reviews non-fitting methods, including rapid lifetime determination, center-of-mass, integral extraction, phasor, τ_A/τ_I, deep learning, and histogram clustering methods. Finally, we conclude this review in Section 4 and speculate on future research directions.

2. Fitting methods

Fitting methods use iterative constrained optimization algorithms to estimate fitting parameters based on a specified decay kinetics model. Feedback is provided at each iteration for updating or terminating the process with a criterion. The weighted least squares estimation (LSE) and the maximum likelihood estimation (MLE) are commonly used for FLIM analysis, and these methods have been compared using synthetic and experimental data [36, 37, 38, 39]. Although it is well known that MLE is more efficient and accurate than LSE, as FLIM data is Poisson distributed rather than Gaussian distributed [40, 41, 42, 43, 44], LSE is still more popular than MLE in FLIM analysis. Attempts have been made to provide robust MLE fitting algorithms, which are counterparts of the Levenberg–Marquardt procedure for LSE [45, 46, 47].

The global fitting method can accelerate analysis by assuming spatial lifetime invariances to reduce the degree of freedom significantly. The Laguerre expansion method can accelerate analysis by converting the nonlinear-fitting problem to a linear-fitting problem to speed up deconvolution procedures.

2.1 Weighted least squares estimation (LSE)

The weighted LSE defines a merit function as [48, 49],

$${\chi }^2=\sum _m\frac{{\left(h_m-{\widehat{h}}_m\right)}^2}{{\sigma }_m^2},$$E5

where h_m is the measured counts, ${\widehat{h}}_m$ is the predicted counts, and σ_m is the count deviation in Bin m. By minimizing ${\chi }^2$, lifetime parameters can be obtained. For Poisson distributed data, the variance is equal to the actual value h_m^true of Bin m, that is, ${\sigma }_m=\sqrt{h_m^{\mathit{\text{true}}}}$. As this actual value is not available, ${\chi }^2$ can be approximated as:

$${\chi }_P^2=\sum _m\frac{{\left(h_m-{\widehat{h}}_m\right)}^2}{h_m},{\chi }_N^2=\sum _m\frac{{\left(h_m-{\widehat{h}}_m\right)}^2}{\mathit{max}\left(h_m,1\right)},$$E6

where ${\chi }_P^2$ and ${\chi }_N^2$ are Pearson’s and modified Neyman’s ${\chi }^2$. Studies show that Neyman weighting underestimates the target answer, whereas Pearson weighting affords an acceptable answer when the total count is more than 1000 [39].

2.2 Maximum likelihood estimation (MLE)

MLE maximizes the probability that the data can occur given a model and a set of parameters. The likelihood L can be expressed as [50, 51].

$$L=\prod _{m}P\left(h_m;{\widehat{h}}_m\right),$$E7

where $P\left(h_m;{\widehat{h}}_m\right)$ is the probability that Bin m has $h_m$ counts if the actual value is ${\widehat{h}}_m$. Because the measurements are Poisson distributed, the Poisson likelihood becomes.

$$L_P=\prod _m\frac{{{\widehat{h}}_m}^{h_m}}{h_m!}{\mathrm{e}}^{-{\widehat{h}}_m}.$$E8

If Eq. (8) is divided by the maximum possible likelihood$L_P\left(h_m;h_m\right)$, then

$$\lambda =L_P\left(h_m;{\widehat{h}}_m\right)/L_P\left(h_m;h_m\right).$$E9

A merit function can be defined as:

$${\chi }_{\mathrm{mle}}^2=-2ln\lambda =2\sum _m\left({\widehat{h}}_m-h_m\right)-2\sum _{m,h_m\ne 0}{h}_{m}ln\left({\widehat{h}}_m/h_m\right).$$E10

Based on the Poisson likelihood function, MLE can ensure unbiased estimations.

2.3 Global fitting (GF)

The global fitting method uses a least squares estimate ${\chi }_{\mathit{GF}}^2$ for all histograms as the merit function,

$${\chi }_{\mathit{GF}}^2=\sum _s\sum _m\frac{{\left(h_m^{\left(s\right)}-{\widehat{h}}_m^{\left(s\right)}\right)}^2}{{\sigma }_m^{\left(s\right)2}},$$E11

where the superscript s represents Histogram s. GF treats the lifetime components τ as constants, but the amplitude A and fraction q as variables for all histograms. There are two strategies for implementing GF, the iterative convolution and the variable projection approaches. The variable projection approach appears to be faster than the iterative convolution method, as investigated in [52].

2.3.1 Iterative convolution

The underlying decay is estimated with.

$${\widehat{f}}_m^{\left(s\right)}={\widehat{A}}^{\left(s\right)}\sum _{p=0}^{P-1}{{\widehat{q}}_p}^{\left(s\right)}{e}^{\frac{-t_m}{{\widehat{\tau }}_p}},$$E12

where ${\widehat{\tau }}_p$ are estimated constant lifetimes for all histograms with ${\widehat{A}}^{\left(s\right)}$ and ${{\widehat{q}}_p}^{\left(s\right)}$ being the parameters for Histogram s.

Then the estimated signal can be expressed as:

$${\widehat{h}}_m^{\left(s\right)}=\sum _{k=0}^m\mathit{ir}{f}_{k-m}·{\widehat{f}}_m^{\left(s\right)},\mathit{m}=0,\dots ,\mathit{M}-1.$$E13

Using Eq. (13), we can minimize Eq. (11) with constrained LSE. The analysis speed is significantly affected by chosen initial conditions. S. Pelet et al. compared different strategies for initial guesses [30] and proposed an efficient image segmentation method.

2.3.2 Variable projection

The idea of global fitting with variable projection is to minimize a projection function that depends only on nonlinear parameters τ and obtain linear parameters A^(s) and q^(s). A matrix whose columns only depend on τ is constructed,

$$\mathit{\Phi }\left(\widehat{\mathit{\tau }}\right)=\left[{\mathit{\phi }}_0\left({\widehat{\tau }}_0\right),\dots ,{\mathit{\phi }}_{P-1}\left({\widehat{\tau }}_{P-1}\right)\right],$$E14

where${\mathit{\phi }}_p\left({\widehat{\tau }}_p\right)={\left[{\phi }_p\left({\widehat{\tau }}_p;t_0\right),\dots ,{\phi }_p\left({\widehat{\tau }}_p;t_{M-1}\right)\right]}^T$, ${\phi }_p\left({\widehat{\tau }}_p;t_m\right)={\sum }_{k=0}^m\mathit{ir}{f}_{k-m}·exp\left(-\frac{t_m}{{\widehat{\tau }}_p}\right)$, and p = 0, …, P-1.

Then the estimated signal can be written as:

$${\widehat{h}}_m^{\left(s\right)}=\sum _{p=0}^{P-1}{\widehat{a}}_p^{\left(s\right)}{\phi }_p\left({\widehat{\tau }}_p;t_m\right),$$E15

where ${{\widehat{a}}_p}^{\left(s\right)}={\widehat{A}}^{\left(s\right)}{{\widehat{q}}_p}^{\left(s\right)}$. Eq. (11) becomes.

$${\chi }_{\mathit{GF}}^2=\sum _s\sum _m\frac{{\left(h_m^{\left(s\right)}-{\sum }_{p=0}^{P-1}{\widehat{a}}_p^{\left(s\right)}{\phi }_p\left({\widehat{\tau }}_p;t_m\right)\right)}^2}{{\sigma }_m^{\left(s\right)2}},$$E16

$${\chi }_{\mathit{GF}}^2=\sum _s{∥{\overline{\mathit{h}}}^{\left(s\right)}-\overline{\mathit{\Phi }}\left(\widehat{\mathit{\tau }}\right){\widehat{\mathit{a}}}^{\left(s\right)}∥}^2,$$E17

where ${\overline{\mathit{h}}}^{\left(s\right)}$ is ${\mathit{h}}^{\left(s\right)}$ weighted by ${\sigma }_m^{\left(s\right)}$, $\overline{\mathit{\Phi }}\left(\widehat{\mathit{\tau }}\right)$ is $\mathit{\Phi }\left(\widehat{\mathit{\tau }}\right)$ weighted by ${\sigma }_m^{\left(s\right)}$, and ${\widehat{\mathit{a}}}^{\left(s\right)}={\left[{{\widehat{a}}_0}^{\left(s\right)},\dots ,{{\widehat{a}}_{P-1}}^{\left(s\right)}\right]}^T$.

For a given $\widehat{\mathit{\tau }}$, Eq. (17) is minimized when ${\widehat{\mathit{a}}}^{\left(s\right)}=\overline{\mathit{\Phi }}{\left(\widehat{\mathit{\tau }}\right)}^{-}{\overline{\mathit{h}}}^{\left(s\right)}$, where $\overline{\mathit{\Phi }}{\left(\widehat{\mathit{\tau }}\right)}^{-}$ is the symmetric generalized inverse of $\overline{\mathit{\Phi }}\left(\widehat{\mathit{\tau }}\right)$. Then Eq. (17) can be rewritten as:

$${\chi }_{\mathit{GF}}^2=\sum _s{∥{\mathit{P}}_{\mathit{\Phi }\left(\widehat{\mathit{\tau }}\right)}^{\perp }{\overline{\mathit{h}}}^{\left(s\right)}∥}^2,$$E18

where ${\mathit{P}}_{\mathit{\Phi }\left(\widehat{\mathit{\tau }}\right)}^{\perp }=\mathit{I}-\overline{\mathit{\Phi }}\left(\widehat{\mathit{\tau }}\right)\overline{\mathit{\Phi }}{\left(\widehat{\mathit{\tau }}\right)}^{-}$ [52], which can be calculated by matrix decomposition of $\overline{\mathit{\Phi }}\left(\widehat{\mathit{\tau }}\right)$ using the QR method. Eq. (18) only has the nonlinear parameter, reducing the minimization parameter space considerably. Once $\widehat{\mathit{\tau }}$ is obtained by minimizing Eq. (18), the linear parameter ${\widehat{\mathit{a}}}^{\left(s\right)}$ can be obtained as a solution of ${\mathit{h}}^{\left(s\right)}=\mathit{\Phi }\left(\widehat{\mathit{\tau }}\right){\widehat{\mathit{a}}}^{\left(s\right)}$. The implementation can be achieved with the VARP2 code [53].

2.4 Laguerre expansion (LE)

Due to the IRF introduced in FLIM data, iterative fitting methods include an enormous amount of convolutions, which is time-consuming [54]. Furthermore, iterative deconvolution methods require acquiring a considerable number of counts, which would increase the acquisition time. Numerous mathematical tools have been devised for deconvolution [55, 56, 57, 58, 59], and LE is faster and more robust than others. LE estimates the underlying fluorescence decay f with an ordered set of discrete-time Laguerre basis functions (LBFs) [60, 61],

$${\widehat{f}}_m=\sum _{l=0}^{L-1}{\widehat{c}}_l{b}_l\left(m;\alpha \right),$$E19

where L and α are the basic parameters and ${\widehat{c}}_l$ is the estimated l^th expansion coefficient. The l^th discrete-time LBF is defined as:

$$b_l\left(m;\alpha \right)={\alpha }^{\frac{\left(m-l\right)}{2}}{\left(1-\alpha \right)}^{\frac{1}{2}}\sum _{i=0}^l{\left(-1\right)}^i\left(\begin{array}{c}m\\ i\end{array}\right)\left(\begin{array}{c}l\\ i\end{array}\right){\alpha }^{l-i}{\left(1-\alpha \right)}^i,$$E20

where 0 < α < 1 and l = 0, …, L-1.

Substituting Eq. (19) for f_m in Eq. (1), the estimated signal becomes:

$${\widehat{h}}_m=\sum _{i=0}^m\sum _{l=0}^{L-1}{\widehat{c}}_l·\mathit{ir}{f}_{m-i}·b_l\left(i;\alpha \right)=\sum _{l=0}^{L-1}{\widehat{c}}_l·{\upsilon }_l\left(m;\alpha \right),$$E21

where ${\upsilon }_l\left(m;\alpha \right)={\sum }_{i=0}^m\mathit{ir}{f}_{m-i}·b_l\left(i;\alpha \right)$. Then, Eq. (5) becomes.

$${\chi }_{\mathit{LE}}^2=\sum _m\frac{{\left(h_m-{\sum }_{l=0}^{L-1}{\widehat{c}}_l·{\upsilon }_l\left(m;\alpha \right)\right)}^2}{{\sigma }_m^2},$$E22

where V = [v₀, …, v_L-1] and v_l = [v_l(0, α), …, v_l(M-1; α)]^T. Eq. (22) can be minimized with the ordinary and constrained LSE, as demonstrated in [61]. Once $\widehat{\mathit{c}}={\left[{\widehat{c}}_0,\dots ,{\widehat{c}}_{L-1}\right]}^T$ is determined, $\widehat{\mathit{f}}$ can be recovered with Eq. (19). Then decay parameters can be extracted from $\widehat{\mathit{f}}$ using the fitting methods mentioned above or the following non-fitting methods.

Setting proper L and α depends on the lifetime dynamic range and the measurement window T = MΔt [32]. Automated Laguerre deconvolution methods have been reported to optimize L and α during the deconvolution routine [62, 63].

3. Non-fitting methods

Fitting methods usually require experienced users to set proper initial conditions and require time-consuming iterative computations, making unsupervised non-fitting methods attractive for robust, real-time FLIM analysis. We review the rapid lifetime determination, center-of-mass, integral extraction, phasor, τ_A/τ_I, deep learning, and histogram clustering methods.

3.1 Rapid lifetime determination (RLD)

RLD needs two time-gated signals for the lifetime estimation of mono-exponential decays. The initial RLD method utilizes two consecutive time-gates of equal widths [64], as shown in Figure 1(a). By integrating the signal within the time-gated windows, the lifetime can be determined using,

Figure 1.

Rapid lifetime determination (RLD) schemes with (a) equal and (b) unequal time-gate widths.

$${\tau }_{\mathit{RLD}}=\frac{-\Delta t_{\mathit{TG}}}{\mathit{ln}\left(D_1/D_0\right)},$$E23

where D₀ and D₁ denote the integrated signal in Gates 1 and 2, and Δt_TG denotes the time-gate width. This calculation is simple and fast; however, it has limited estimation accuracy and precision in a narrow lifetime range determined by the gate width. To optimize RLD, several strategies with overlapping and unequal time gates, as shown in Figure 1(b), have been proposed [65, 66, 67], suggesting optimized gating schemes. In optimized RLD schemes, the ratio of D₁ and D₀ can be expressed as:

$$\mathit{DR}=\frac{D_1}{D_0}=\frac{\mathit{exp}\left(-Y\frac{{\mathit{\Delta t}}_{\mathit{TG}}}{\tau }\right)-\mathit{exp}\left[-\left(Y+W\right)\frac{{\mathit{\Delta t}}_{\mathit{TG}}}{\tau }\right]}{1-\mathit{exp}\left(-\frac{{\mathit{\Delta t}}_{\mathit{TG}}}{\tau }\right)},$$E24

where Δt_TG is the width of Gate 1 and Y and W are the rising edge and width coefficient of Gate 2. Lifetimes τ_RLD can be solved by applying Newton’s method.

RLD can deal with multiple types of FLIM data, including TCSPC, streak camera, time-gated counting, and time-gated imaging data. However, RLD performs differently for uncorrelated and correlated datasets, as analyzed in ref. [66]. For multi-exponential decays, RLD can provide average lifetimes which are neither intensify- nor amplitude-weighted. Additionally, with three-time gates, the lifetime estimation dynamic range of RLD can be expanded for mono-exponential decays [21]. With four time-gates, we can extract lifetime components of bi-exponential decays [68, 69]. IRF is usually neglected in RLD, introducing estimation bias, especially for short lifetimes comparable with IRF’s width.

3.2 Center-of-mass method (CMM)

CMM provides intensity-weighted average lifetimes for multi-exponential decays [70, 71, 72],

$${\tau }_{\mathit{CMM}}=\frac{{\int }_0^{\infty }t·h\left(t\right)\mathit{dt}}{{\int }_0^{\infty }h\left(t\right)\mathit{dt}}-\frac{{\int }_0^{\infty }t·\mathit{irf}\left(t\right)\mathit{dt}}{{\int }_0^{\infty }\mathit{irf}\left(t\right)\mathit{dt}}=\frac{{\sum }_{p=0}^{P-1}{q}_p{\tau }_p^2}{{\sum }_{p=0}^{P-1}{q}_p{\tau }_p}\approx \frac{{\sum }_{m=0}^{M-1}{t}_m·h_m}{{\sum }_{m=0}^{M-1}{h}_m}-\frac{{\sum }_{m=0}^{M-1}{t}_m·\mathit{ir}{f}_m}{{\sum }_{m=0}^{M-1}\mathit{ir}{f}_m}.$$E25

Li et al. proposed two versions of bi-exponential CMM (BCMM) that provide lifetime components information [73].

3.3 Integral extraction method (IEM)

For IEM, deconvolution is required to obtain $\widehat{\mathit{f}}$ with which the average lifetime can be determined as [74, 75].

$${\tau }_{\mathit{IEM}}=-\frac{{\int }_0^{\infty }g\left(t\right)\mathit{dt}}{{\int }_0^{\infty }{g}^{’}\left(t\right)\mathit{dt}}\approx \sum _{p=0}^{P-1}{q}_p{\tau }_p\approx -\frac{{\sum }_{m=0}^{M-1}{S}_m·{\widehat{f}}_m}{{\sum }_{m=0}^{M-1}\frac{{\widehat{f}}_m-{\widehat{f}}_{m-1}}{\mathit{\Delta t}}}=-\frac{\mathit{\Delta t}{\sum }_{m=0}^{M-1}{S}_m·{\widehat{f}}_m}{{\widehat{f}}_{M-1}-{\widehat{f}}_0},$$E26

where S_m = [1/3, 4/3, 2/3, …, 4/3, 1/3] are the coefficients for numerical integration based on Simpson’s rule,

$$g\left(t\right)=A\sum _{p=0}^{P-1}{q}_p{\tau }_p{e}^{-\frac{t}{{\tau }_p}}\left[1-e^{-\frac{\mathit{\Delta t}}{{\tau }_p}}\right].$$E27

τ_IEM is an estimator of amplitude-weighted lifetimes.

3.4 Phasor method

The phasor method transforms each histogram into a phasor, like a vector. The sine-cosine transforms of decays are represented in a phasor plot as a two-dimensional histogram [76, 77, 78, 79, 80]. Phasor components g and s for time-domain FLIM can be expressed as:

$$g=\frac{{\int }_0^{\infty }f\left(t\right)·cos\left(\mathit{\omega t}\right)\mathit{dt}}{{\int }_0^{\infty }f\left(t\right)\mathit{dt}}=\frac{{\sum }_{p=0}^{P-1}\frac{q_p{\tau }_p}{1+{\omega }^2{{\tau }_p}^2}}{{\sum }_{p=0}^{P-1}{q}_p{\tau }_p}=\frac{\mathit{Re}+s\bullet B}{A},$$E28

$$s=\frac{{\int }_0^{\infty }f\left(t\right)·sin\left(\mathit{\omega t}\right)\mathit{dt}}{{\int }_0^{\infty }f\left(t\right)\mathit{dt}}=\frac{{\sum }_{p=0}^{P-1}\frac{\omega q_p{{\tau }_p}^2}{1+{\omega }^2{{\tau }_p}^2}}{{\sum }_{p=0}^{P-1}{q}_p{\tau }_p}=\frac{A\bullet \mathit{Im}-B\bullet \mathit{Re}}{A^2+B^2},$$E29

where

$$\mathit{Re}=\frac{{\int }_0^{\infty }h\left(t\right)·cos\left(\mathit{\omega t}\right)\mathit{dt}}{{\int }_0^{\infty }h\left(t\right)\mathit{dt}}\approx \frac{{\sum }_{m=0}^{M-1}{h}_m·cos\left(\omega t_m\right)}{{\sum }_{m=0}^{M-1}{h}_m},$$E30

$$\mathit{Im}=\frac{{\int }_0^{\infty }h\left(t\right)·sin\left(\mathit{\omega t}\right)\mathit{dt}}{{\int }_0^{\infty }h\left(t\right)\mathit{dt}}\approx \frac{{\sum }_{m=0}^{M-1}{h}_m·sin\left(\omega t_m\right)}{{\sum }_{m=0}^{M-1}{h}_m},$$E31

$$A=\frac{{\int }_0^{\infty }\mathit{irf}\left(t\right)·cos\left(\mathit{\omega t}\right)\mathit{dt}}{{\int }_0^{\infty }\mathit{irf}\left(t\right)\mathit{dt}}\approx \frac{{\sum }_{m=0}^{M-1}\mathit{ir}{f}_m·cos\left(\omega t_m\right)}{{\sum }_{m=0}^{M-1}\mathit{ir}{f}_m},$$E32

$$B=\frac{{\int }_0^{\infty }\mathit{irf}\left(t\right)·sin\left(\mathit{\omega t}\right)\mathit{dt}}{{\int }_0^{\infty }\mathit{irf}\left(t\right)\mathit{dt}}\approx \frac{{\sum }_{m=0}^{M-1}\mathit{ir}{f}_m·sin\left(\omega t_m\right)}{{\sum }_{n=0}^{M-1}\mathit{ir}{f}_m}.$$E33

Phasor plots’ interpretation is usually user-dependent by manually selecting regions of interest in a phasor plot to find corresponding regions in fluorescence images. Based on the feature that pixels with similar fluorescence decays tend to congregate and form a cluster in a phasor plot, machine learning techniques have been developed with clustering methods [81] to automatically organize phasors into sensible groupings. Up to four-lifetime components can be resolved from a phasor plot using the rule of a linear combination of phasors [82], a graphical approach [83], and a computational nonlinear minimization algorithm [83].

A weighted average lifetime, τ_Phasor, can also be derived using phasors,

$${\tau }_{\text{Phasor}}=\frac{s}{\mathit{g\omega }}=\frac{{\sum }_{p=0}^{P-1}\frac{q_p{{\tau }_p}^2}{1+{\omega }^2{{\tau }_p}^2}}{{\sum }_{p=0}^{P-1}\frac{q_p{\tau }_p}{1+{\omega }^2{{\tau }_p}^2}},$$E34

where ω = 2π/T, T = MΔt is the measurement window. The weights of τ_Phasor are $\frac{q_p{\tau }_p}{1+{\omega }^2{{\tau }_p}^2}$. If τ_p < < T, then the weights are approximately equal to $q_p{\tau }_p$, i.e. τ_Phasor is close to intensity-weighted lifetimes.

3.5 τ_A/τ_I method

We proposed the τ_A/τ_I method, a multi-exponential decay visualization method using two types of average lifetimes, τ_I and τ_A. As mentioned previously (discussed in ref. [35]), a fluorescence decay can be approximated by a bi-exponential decay model, so that the ratio of τ_A and τ_I can be expressed as:

$$\frac{{\tau }_A}{{\tau }_I}=\frac{{\left(1+q_1\left(R-1\right)\right)}^2}{1+q_1\left(R^2-1\right)},$$E35

where R = τ₁/τ₂. The distribution of τ_A/τ_I (Figure 2) shows that when R ≈ 1 or q₁ ≈ 0 or 1, τ_A/τ_I ≈ 1. With a decrease in R or an increase of q₁, τ_A/τ_I decreases. Therefore, the ranges of q₁ and R of a pixel can be determined by τ_A/τ_I.

Figure 2.

Distribution of τ_A/τ_I with q₁ = 0 ∼ 1 and R = 0 ∼ 1 [35].

With a TCSPC dataset of tSA201 cells, Figures 3 and 4 show the results of the selected pixels within different τ_A/τ_I ranges in (a) τ_I, (b) τ_A, (c) τ_A/τ_I images, (d) histograms, (e) phasor plots, and (f) τ_A/τ_I plots. For the pixels within τ_A/τ_I = 0.2 ∼ 0.5 (Figure 3), the histograms clearly show that τ_A is much smaller than τ_I, which means the difference between τ₁ and τ₂ is significant. Figure 3 (f) shows that the ranges of q₁ and R are approximately 0.5 ∼ 1 and 0 ∼ 0.2, respectively. For the pixels with τ_A/τ_I = 0.5 ∼ 1 (Figure 4), τ_A is closer to τ_I, meaning the pixels have decays close to mono-exponential.

Figure 3.

(a) τ_I-intensity image, (b) τ_A-intensity image, (c) τ_A/τ_I ratio image, (d) histograms of τ_I (yellow) and τ_A (blue), (e) phasor plot, and (f) distribution of τ_A/τ_I of the selected pixels in (c) with τ_A/τ_I = 0.2 ∼ 0.5 [35].

Figure 4.

(a) τ_I-intensity image, (b) τ_A-intensity image, (c) τ_A/τ_I ratio image, (d) histograms of τ_I (yellow) and τ_A (blue), (e) phasor plot, and (f) distribution of τ_A/τ_I of the selected pixels in (c) with τ_A/τ_I = 0.5 ∼ 1 [35].

τ_A/τ_I is an intuitive tool for visualizing multi-exponential decays in a lifetime image. Separating the average lifetime images with τ_A/τ_I is easier than phasor plots because τ_A/τ_I is one-dimensional, and phasors are two-dimensional. Furthermore, τ_A/τ_I can intuitively show the q₁ and R ranges.

3.6 Deep learning (DL)

Recently, DL-based FLIM analysis methods have been reported. DL features hierarchical representation learning by extracting high-level features through multiple nonlinear transformations of low-level features. DL shows a powerful ability to learn complex data and functions. The advent of DL breaks the conventional “model-driven” paradigm and offers us a new “data-driven” approach to solving general optimization problems. Given sufficient labeled training data, DL algorithms can directly map the raw input data to their corresponding results, thus avoiding time-consuming iterative optimization processes. Wu et al., first employed DL with a multi-layer perception (MLP) network using one input layer, one output layer, and two hidden layers [84]. The input layer has 57 entries depending on the number of time-bins in a histogram. The output layer has four neurons mapping four-lifetime parameters for bi-exponential decays. The investigation reveals that MLP can provide comparable or better performances and generate lifetime images at least 180-fold faster than conventional LSE. However, IRF was not considered in their work. Smith et al. proposed a 3D convolutional neural network (CNN) architecture named fluorescence lifetime imaging network (FLI-Net) to quantify fluorescence decays at fast speeds [85]. FLI-Net is designed and trained for TCSPC FLIM, and gated intensity charged-coupled device (ICCD) based FLIM. FLI-Net’s input is a 3D data cube (x, y, t), and output is bi-exponential decay parameters. Synthetic data is used for training, avoiding acquiring massive training datasets experimentally. FLI-Net is about 30-fold faster than SPCImage, the widely used FLIM processing software. FLI-Net outperforms LSE and MLP, especially with low photon counts range from 25 to 100. Xiao et al. proposed an easier and faster trained 1D CNN architecture named 1D-ConvResNet [86]. Compared with 2D or 3D CNNs, 1D-ConvResNet is more hardware-friendly and can be implemented on field-programmable gate array (FPGA) devices. Synthetic data is used for training and the training time with a CPU (Intel i7–4790) is about 0.5 hours, 8-fold faster than FLI-Net with a GPU (NVIDIA TITAN Xp GPU). Experimental FLIM datasets with an intensity threshold of 100 counts per pixel were used for validation, and analyzing a 256 × 256 image takes several seconds on a laptop. To generate high-quality FLIM images under photon-starved conditions (50 counts per pixel), Chen et al. introduced a method called flimGANE [87]. It was derived from the Wasserstein GAN algorithm, where an “artificial” high-photon-count fluorescence decay histogram can be produced with a generator with a low-photon-count input. Using a well-trained generator and an estimator, a low-quality decay histogram can be mapped to a high-quality counterpart and provides bi-exponential lifetime parameters within 0.32 ms/pixel using a CPU, which is 258-fold faster than LSE and 2800-fold faster than MLE in generating a 512 × 512 FLIM image. However, it takes up to 500 h to fully train the network. Additionally, using time-resolved single-pixel datasets, a deep CNN named Net-FLICS (fluorescence lifetime imaging with compressed sensing) was reported [88].

DL methods for spatial resolution enhancement of FLIM images are also developed, including SRI-FLIMnet for reconstructing high-resolution images from low-resolution 3D FLIM data [89] and CNN-based denoising method removing noise in phasor plots after the K-means clustering segmentation [90].

Although DL is promising in real-time FLIM analysis even under starve photon-starved conditions and needs no user-dependent initial conditions, it suffers from long training times (hours to days) when retaining is required due to the change of IRF. Whenever the laser source for excitation and the detector are changed, the IRF can also vary, and the network should be retrained, reducing the universality. Recently, many attempts have been made to alleviate this disadvantage. For example, Zang proposed an online neural network training method using the extreme learning machine (ELM) [91]. ELM does not require a back-propagation process in the training phase. Therefore, it provides a much faster training speed, enabling online network training for any system configuration. The phasor coordinates of a fluorescence decay can also serve as the network’s inputs, which can be viewed as a simple feature engineering to reduce the training time significantly [92]. Despite some trade-offs in current DL algorithms, it is no doubt that DL has excellent potential in a wide range of FLIM-related applications. DL algorithms can also be implemented on edge-computing platforms like FPGAs and smartphones to develop intelligent and portable FLIM devices [93].

3.7 Histogram clustering (HC)

The HC method devised by Li et al. can improve FLIM analysis speed and accuracy by sorting histograms with similar profiles in a dataset into several clusters and significantly reducing the number of histograms to be analyzed [94]. HC implements clustering with two features of a histogram. It is a preprocessing method that can work with the LDAs mentioned above. The performances for producing decay parameter images without and with HC using synthetic and experimental datasets were investigated [94]. The execution time t_exe and the mean squared error (MSE) of a FLIM dataset following a bi-exponential decay model with 150 × 150 pixels and 256 time-bins are shown in Table 1. LE-LSE and LE-IEM stand for LSE and IEM with the Laguerre expansion method; IC and VP represent iterative convolution and variable expansion global fitting approaches.

LDA	texe (s)	MSE
		q1	τ1 (ns^2)	τ2 (ns^2)	τA (ns^2)	τI (ns^2)
Without HC
LE-LSE	389.45	0.019	0.173	0.198	0.110	0.027
LE-IEM	62.30	X	X	X	0.102	X
CMM	0.20	X	X	X	X	0.185
IC	724.82	0.100	X	X	0.678	1.098
VP	3.34	0.033	X	X	0.122	0.178
With HC
LE-LSE	3.36	0.011	0.102	0.104	0.037	0.025
LE-IEM	0.63	X	X	X	0.038	X
CMM	0.20	X	X	X	X	0.180
IC	11.85	0.017	X	X	0.102	0.050
VP	0.31	0.014	X	X	0.048	0.093

Table 1.

t_exe and MSE evaluated by lifetime determination algorithms without and with histogram clustering [94].

For different output types, we suggested the fastest FLIM analysis methods: 1) LE-LSE with HC for all lifetime component images with t_exe = 3.36 s, 116-fold shorter than t_exe without HC; 2) VP with HC for constant lifetimes, q₁, τ_A, and τ_I images with t_exe = 0.31 s, 10-fold shorter than t_exe without HC; and 3) LE-IEM with HC as the second choice for τ_A with t_exe = 0.63 s, 98-fold shorter than t_exe without HC, and CMM as the second choice for τ_I with t_exe = 0.2 s without or with HC (biased if the most significant lifetime > T/4). The analysis was conducted in MATLAB, and it can be translated to C or other environments to speed up the analysis.

HC can benefit applications demanding real-time FLIM, such as clinical diagnosis and fast screening. In the future, deep learning methods can be employed for unsupervised histogram clustering.

4. Conclusions

Lifetime determination algorithms (LDAs) play a vital role in FLIM analysis. The results provided by diverse LDAs profoundly impact the interpretation of observed phenomena. Numerous approaches have been developed with notable features in terms of photon efficiency, estimation speed, accuracy, and precision. Fast FLIM analysis under low-photon conditions is universally desirable, especially for real-time and single-molecular FLIM. This mini review focuses on the popular and cutting-edge LDAs used in time-domain FLIM. The progress in LDAs has shown possibilities for real-time FLIM under low photon conditions.

Fitting methods estimate lifetime parameters by solving an iterative minimization problem. The least squares estimation (LSE) introduces bias even with high photon counts (>1000), while the maximum likelihood estimation (MLE) is unbiased and photon efficient, as FLIM data is Poisson distributed. With robust algorithms, MLE instead of LSE can be used as a standard. The iterative convolution of IRF with exponential decays in LSE and MLE is time-consuming. The global fitting and Laguerre expansion methods have been proposed to accelerate analysis.

Non-fitting methods are attractive, as they can be much faster and more hardware-friendly than fitting methods because of simple calculations. However, non-fitting methods for a specific application should be careful. RLD using two time-gates is suitable for mono-exponential decays but may introduce a bias for lifetimes comparable with IRF’s width. CMM and IEM deliver intensity- and amplitude-weighted average lifetimes, respectively. The phasor method can graphically show exponential decays in a phasor plot. The τ_A/τ_I method is a multi-exponential decay visualization method that uses two types of average lifetimes, which can intuitively show lifetime components’ fraction and lifetime ratio. The deep learning method employs fully connected networks, 3D CNN, 1D CNN, and GAN algorithms in rapid FLIM analysis under photon-starved conditions, down to 50 counts per pixel. However, its long retraining time (hours to days) limits its universality when IRF is changed. The histogram clustering (HC) method suggests that histograms having similar profiles can be clustered for one calculation. By embedding HC, the algorithms mentioned above can be significantly accelerated with an execution time down to 30 μs/pixel.

In summary, the existing LDAs reveal different features in FLIM analysis, and it is hard to judge which one is replaceable. We believe that software with multiple embedded LDAs is desirable for end-users, providing comprehensive and complementary information for robust analysis. In addition, the newly developed “data-driven” DL-based LDAs and LDA accelerating method HC are promising, as they can not only enhance speed but also accuracy, especially for photon-starved applications. More efforts are expected to shrink retraining time for DL and increase clustering performance for HC via deep learning.

Acknowledgments

Supported by the Scientific Instrument Developing Project of the Chinese Academy of Sciences (GJJSTD20220006), Science Foundation of the Chinese Academy of Sciences (CXJJ-21S006), China Scholarship Council, Biotechnology and Biological Sciences Research Council (BB/V019643/1) and Medical Research Scotland.

Conflict of interest

The authors declare no conflict of interest.

Abbreviations

FLIM	Fluorescence Lifetime Imaging
FRET	Förster Resonance Energy Transfer
IRF	Instrument Response Function
TCSPC	Time-Correlated Single-Photon Counting
LDA	Lifetime Determination Algorithm
LSE	Least Squares Estimation
MLE	Maximum Likelihood Estimation
GF	Global Fitting
LE	Laguerre Expansion
LBF	Laguerre Basis Function
RLD	Rapid Lifetime Determination
CMM	Center-of-Mass Method
IEM	Integral Extraction Method
DL	Deep Learning
HC	Histogram Clustering
MLP	Multi-Layer Perception
CNN	Convolutional Neural Network
GAN	Generative Adversarial Network
ELM	Extreme Learning Machine
FPGA	Field-Programmable Gate Array