Open Access
Add to BibSonomy
Abstract
We propose a histogram clustering (HC) method to accelerate fluorescence lifetime imaging (FLIM) analysis in pixel-wise and global fitting modes. The proposed method’s principle was demonstrated, and the combinations of HC with traditional FLIM analysis were explained. We assessed HC methods with both simulated and experimental datasets. The results reveal that HC not only increases analysis speed (up to 106 times) but also enhances lifetime estimation accuracy. Fast lifetime analysis strategies were suggested with execution times around or below 30 μs per histograms on MATLAB R2016a, 64-bit with the Intel Celeron CPU (2950M @ 2GHz).
Published by The Optical Society under the terms of the Creative Commons Attribution 4.0 License. Further distribution of this work must maintain attribution to the author(s) and the published article's title, journal citation, and DOI.
1. Introduction
Fluorescence lifetime imaging (FLIM) [1] is a crucial technique for assessing microenvironments of fluorophores, such as pH, Ca$^{2+}$, O$_2$, viscosity, or temperature [2–5]. Combining with Förster Resonance Energy Transfer (FRET) techniques [6], FLIM can be a powerful "quantum ruler" to measure protein conformations and interactions [7]. In contrast to fluorescence intensity imaging, FLIM is independent of fluorescence intensities and fluorophore concentrations, making FLIM a robust quantitative imaging technique for life sciences applications [8,9], medical diagnosis [10], drug developments [11,12], and flow diagnosis [13–15].
A fluorescence decay is usually modeled as a sum of exponential decay functions:
where $A$ is the amplitude, $q_p$ and $\tau _p$ are the fraction and lifetime of the $p^{th}$ component, $p = 1, \dots , P$. In vector forms, $\boldsymbol{q} = [q_1, \dots , q_P]^T$ and $\boldsymbol{\tau} = [\tau _1, \dots , \tau _P]^T$. In reality, the measured signal is a convolution of $f(t)$ and the instrument response function (IRF) $ir\!f(t)$, where $\epsilon (t)$ is noise.FLIM analysis is equivalent to solving the inverse problem from Eq. (2) with the measured $h(t)$ to obtain $\boldsymbol{q}$ and $\boldsymbol{\tau}$. FLIM experiments can be conducted either in time- or frequency-domain manners [8]. In time-domain approaches, samples are illuminated with ultrashort laser pulses. $h(t)$ is measured using a time-correlated single-photon counting (TCSPC) system [16,17] with photomultiplier tubes, delay line anode detectors [18] or single-photon avalanche diodes [19] in scanning or widefield modes. $h(t)$ can also be measured with time-gated cameras [20,21] and streak cameras [22,23]. There are also frequency-domain approaches [24,25], but we will focus on time-domain approaches in this report.
A fluorescence decay histogram measured by a TCSPC system can be:
whereWith $\boldsymbol{h}$ and $\boldsymbol{ir\!f}$ already measured, $A$, $\boldsymbol{q}$ and $\boldsymbol{\tau}$ can be extracted with a lifetime determination algorithm by solving a nonlinear minimization problem $\arg \min \|\boldsymbol{h} - \hat {\boldsymbol{h}}\|^2$, where $\hat {\boldsymbol{h}}$ is the estimated histogram. The iterative convolution (IC) is commonly used with the least-squared method (LSM) [26,27] for solving the inverse problem, denoted as IC-LSM. Still, IC-LSM suffers from low photon efficiency and slow analysis. Several deconvolution approaches have been developed to enhance the analysis, such as the Laguerre expansion [28–30], the non-fitting and the global fitting [31,32] methods. The Laguerre expansion methods speed up deconvolution procedures by converting the nonlinear-fitting problem to a linear-fitting problem estimating a Laguerre basis set’s expansion coefficients. The non-fitting methods, including the centre-of-mass method (CMM) [33–35], the integral extraction method (IEM) [36,37], the phasor method [38,39], or the rapid lifetime determination method [40,41], can provide rapid average lifetime analysis [42]. The global fitting methods can accelerate analysis by changing the estimation mode from the pixel-wise mode to a global fitting mode and using spatial lifetime invariances of fluorescent species in an image to reduce the degree of freedom significantly. There are two strategies, IC [31] and the variable projection (VP) method [32], for implementing global fitting.
However, the Laguerre expansion, the non-fitting and the global fitting methods are not fast enough for growing demands for real-time FLIM. This work presents a histogram clustering (HC) method for improving FLIM analysis in analysis speed and accuracy. Section 2 (Methods) summarizes the workflows for decay parameter image reconstructions with and without HC. We will then introduce and demonstrate the HC method. Besides the algorithms used in this work, HC can also accelerate other algorithms, such as the maximum likelihood method [43], Bayesian methods [44], and deep-learning methods [45]. In Section 3 (Results and Analysis), synthetic and experimental TCSPC datasets will be used to evaluate the HC method’s performances. Suggestions of the fastest algorithms for different outcomes will be given.
2. Methods
2.1 Modes for decay parameter image reconstructions
Figure 1(a) shows the Pixel-Wise (PW) mode’s workflow. $N_{v\!p}$ is the number of valid pixels in a TCSPC dataset whose intensities are beyond a threshold. Histogram $s$, denoted as $\boldsymbol{h}^{(s)}$, is sent into an algorithm for PW along with $\boldsymbol{ir\!f}$, $s = 1, \dots , N_{v\!p}$. After $N_{v\!p}$ histograms are analyzed pixel-by-pixel, decay parameter images are produced. The total execution time $t_{exe}^{P\!W} = N_{v\!p} \times t_{A}^{P\!W}$, $t_{A}^{P\!W}$ is the adopted algorithm’s execution time for PW.
Fig. 1. Flow diagrams of (a) the pixel-wise (PW) mode and (b) the global fitting mode for all pixels (GF-P).
Figure 1(b) shows the workflow of the Global Fitting mode for all Pixels, denoted as GF-P. Instead of estimating decay parameters individually for each pixel, GF-P assumes lifetimes $\boldsymbol{\tau}$ are constants, and $A$ and $\boldsymbol{q}$ vary across the image. $t_{exe}^{G\!F\!-P} = t_ {A}^{G\!F}$. $t_{A}^{G\!F}$ is the adopted algorithm’s execution time for GF.
Figure 2 shows the workflows where the HC method is embedded. Figure 2(a) shows the Cluster-Wise (CW) mode, which combines PW and HC; likewise, Fig. 2(b) shows the Global Fitting mode for all Clusters (GF-C), which combines GF-P and HC.
Fig. 2. Flow diagrams of (a) the cluster-wise mode (CW) and (b) the global fitting mode for all clusters (GF-C).
In CW, $N_{v\!p}$ histograms are first sorted by HC, whose execution time is $t_{HC}$, into $N_c$ classes with $N_c$ cluster-histograms $\bar{\boldsymbol{h}}^{(s)}$, $s = 1,\dots , N_c$. $\bar{\boldsymbol{h}}^{(s)}$ is used to estimate decay parameters for Cluster $s$. Then, the decay parameters are assigned to the corresponding cluster’s pixels with a parameter assignment function, whose execution time is $t_{P\!A}$. Therefore, $t_{exe}^{C\!W} = t_{H\!C} + N_c \times t_{A}^{P\!W} + t_{P\!A}$.
In GF-C, $N_{v\!p}$ histograms are processed with HC first, and the output $N_c$ histograms are sent into an algorithm for GF. Decay parameters for all clusters are obtained and assigned to the pixels in corresponding clusters with the parameter assignment function. Therefore, $t_{exe}^{G\!F\!-C} = t_{H\!C} + t_{A}^{G\!F} + t_{P\!A}$.
The algorithms used in this work are reviewed in Supplement 1.
2.2 Histogram clustering
In reality, there are always many pixels within the field of view showing similar histogram profiles, and it’s unnecessary to analyze them individually because it would be time-consuming. The idea of HC is to sort histograms with similar profiles and to divide them into $N_c$ clusters. If histograms have similar decay profiles, they are supposed to show similar decay parameters. Therefore, we can average similar decay profiles in a cluster into one profile to estimate decay parameters and then assign them to all pixels. With this arrangement, we only need to process $N_c$ instead of $N_{v\!p}$ histograms. HC significantly speedups FLIM analysis.
For simplicity, we only discuss bi-exponential decays widely used in practice. Figure 3(a) shows an IRF and normalized signal profiles $h(t)$ following a bi-exponential decay model, and Fig. 3(b) shows corresponding cumulative signals, $H(t)=\int _0^t h(t)dt$, which is not sensitive to Poisson noise due to the integration. Signal decay parameters are also labelled in Fig. 3(a). If we choose an intensity bound, $I_{bound}$, then each signal has a corresponding time delay, $t_{I\!b}$, to reach $I_{bound}$, as shown in Fig. 3(b).
Fig. 3. Illustrations of (a) an IRF and normalized signal profiles $h(t)$ following a bi-exponential decay model and (b) cumulative signal profiles $H(t)$.
It is straightforward for mono-exponential decays, $f(t) = A\exp {(-t/\tau )}$, that $t_{I\!b}$ has an approximately linear relationship with $\tau$. Figure 4 shows $t_{I\!b}$ curves with different IRFs which introduce time-shifts (assuming that IRFs for all histograms are the same in a scanning system). If a multichannel sensor is used, IRF alignments are required before using HC.
Fig. 4. $t_{I\!b}$ of signals following mono-exponential models depending on $\tau$ under different $ir\!f(t)$.
However, it is less straightforward for bi-exponential decays. Thus, we used numerical methods to conduct analysis, as shown in Fig. 5, in which three cases were simulated to explain how the proposed concept works. Case A has $q_1$ = 0.5, $0.1 \leq \tau _1 \leq 3$ ns, and $\tau _2$ = 3 ns as shown in Fig. 5(a); Case B has $0 \leq q_1 \leq 1$, $\tau _1$ = 0.5 ns, and $\tau _2$ = 3 ns as shown in Fig. 5(b); Case C has $q_1$ = 0.5, $\tau _1$ = 0.5 ns, and $0.5 \leq \tau _2 \leq 10$ ns as shown in Fig. 5(c). The IRF follows a Gaussian distribution with an FWHM of 0.5 ns.
Fig. 5. $t_{I\!b}$ for (a) Case A: $q_1$ = 0.5, $\tau _1$ = 0.1 $\sim$ 3 ns, and $\tau _2$ = 3 ns, (b) Case B: $q_1$ = 0 $\sim$ 1, $\tau _1$ = 0.5 ns, and $\tau _2$ = 3 ns, (c) Case C: $q_1$ = 0.5, $\tau _1$ = 0.5 ns, and $0.5 \leq \tau _2 \leq 10$ ns with different $I_{bound}$.
For Case A, $t_{I\!b}$ is not monotonic with $\tau _1$, and the monotonic range and the slope are functions of $I_{bound}$. As $I_{bound}$ increases, the slope increases with a smaller monotonic range. The profiles with $\tau _1$ outside the range are wrongly sorted into a cluster with a larger $\tau _1$. The monotonic ranges for $I_{bound}$ = 0.2 and 0.6 are 0.4 $\sim$ 3 ns and 1 $\sim$ 3 ns, respectively. For Cases B and C, $t_{I\!b}$ is monotonic (decreasing and increasing) with $q_1$ and $\tau _2$ for all $I_{bound}$, respectively. For the signals like Cases A $\sim$ C (which only have one variable), we can cluster the signals by $t_{I\!b}$ with a proper $I_{bound}$ considering the monotonic range. For example, for Case A, if the shortest lifetime is around 0.5 ns, $I_{bound}$ = 0.2 is a proper choice; for Cases B and C, $I_{bound}$ can be set arbitrarily in 0.1 $\sim$ 1. We use $I_{bound}$ = 0.2 hereafter.
However, it is not realistic that the signals in a dataset have one variable and two constant decay parameters. For instance, in FRET-FLIM applications, donors without FRET have a constant lifetime and donors interacting with acceptors have shorter lifetimes due to FRET. Therefore, the short and long lifetimes, $\tau _1$ and $\tau _2$, are donors’ lifetimes with and without FRET, and $q_1$ is the portion of the donors undergoing FRET among all donors. $q_1$ and $\tau _1$ are variables depending on FRET efficiency.
For FRET-FLIM datasets, such as Case D: $q_1$ = 0 $\sim$ 1, $\tau _1$ = 0.5 $\sim$ 3 ns, and $\tau _2$ = 3 ns, with two variables, it is not enough to divide the histograms only depending on $t_{I\!b}$, as histograms with different profiles would have the same $t_{I\!b}$ and be wrongly divided into one cluster. Figure 6(a) shows the resulting clusters ( $N_1$ = 15) in different colors for Case D with $M$ = 256 depending on $t_{I\!b}$. Figures 6(b) and (c) show the cumulative signals in Clusters 14 and 7 in red, respectively, and the averaged cumulative histograms $\bar {\boldsymbol{H}}$ for the clusters (green dash lines). When $q_1$ $\sim$ 0 (or $q_1$ $\sim$ 1) or $\tau _1$ $\sim$ $\tau _2$ (such as Clusters 14 and 15), the signals are nearly mono-exponential and have similar profiles, as shown in Fig. 6(b). However, the signals for other clusters (for example, Cluster 7) have the same $t_{I\!b}$, but the profiles after $t_{I\!b}$ diverge. At $t = t_{bound}$ in Fig. 6(c), the cumulative intensity $I_{t\!b}$ in Cluster 7 is within [0.7, 0.9]. Therefore, we can further divide each cluster into $N_2$ sub-clusters depending on $I_{t\!b}$.
Fig. 6. For Case D with $M$ = 256, (a) clusters, cumulative signals (red solid lines) and averaged cumulative signals (green dash lines) for (b) Cluster 14 and (c) Cluster 7.
Setting a larger $N_2$ can result in a higher clustering precision, which means histogram profiles in one cluster are more similar. Another way to increase clustering precision is interpolating $\boldsymbol{h}$ with $M$ time-bins to $\boldsymbol{h}^{inter\!p}$ with $M \cdot N_{inter\!p}$ time-bins. $N_{inter\!p}$ ($\geq$ 1) is an interpolation factor, and $\boldsymbol{h}^{inter\!p}$ can be expressed as
Figure 7(a) shows the clusters for Case A with $M$ = 256, $N_{inter\!p}$ = 1, and $N_1$ = 6. Figure 7(b) shows the clusters for $M$ = 256, $N_{inter\!p}$ = 2, and $N_1$ = 11. The histograms in each cluster have a smaller range of $\tau _1$, leading to a higher clustering precision.
Fig. 7. For Case A, clusters with (a) $M$ = 256 and $N_{inter\!p}$ = 1 and (b) $M$ = 256 and $N_{inter\!p}$ = 2.
The HC workflow is summarized in Fig. 8(a). There are three steps: 1) depending on $t_{I\!b}$, $N_{v\!p}$ histograms are divided into $N_1$ clusters, which can be adjusted by setting $N_{inter\!p}$; 2) depending on $I_{tb}$, histograms in each of the $N_1$ clusters are further divided into $N_2$ sub-clusters, and $N_c = N_1 \times N_2$ clusters are finally produced; 3) $N_c$ histograms are generated by obtaining the averaged histogram in each cluster.
Fig. 8. (a) Workflow of the HC method. (b) Boxplot of $\chi ^2$ for Case D with different $N_{inter\!p}$ and $N_2$.
To assess HC in terms of $N_{inter\!p}$ and $N_2$, we can define:
3. Results and analysis
Synthetic and experimental TCSPC datasets were used to assess the performances of HC. Table 1 summarizes different output types of algorithms: (1) the fitting method LE-LSM in PW and CW modes can produce $\boldsymbol{q}$, $\boldsymbol{\tau}$, $\tau _A$, and $\tau _I$ images; (2) the non-fitting methods LE-IEM and CMM in PW and CW modes can produce $\tau _A$ and $\tau _I$ images; (3) the global fitting methods IC and VP in GF-P and GF-C modes can produce $\boldsymbol{q}$, $\tau _A$, and $\tau _I$ images and constant $\boldsymbol{\tau}$. $\tau _A$ and $\tau _I$ are two types of average lifetimes, amplitude- and intensity-weighted average lifetimes,
Using IEM to calculate $\tau _A$ requires conducting the deconvolution first. The Laguerre expansion method is employed for deconvolution when IEM is used; therefore, we denote the whole process as LE-IEM.
Table 1. Outputs of algorithms for different modes.
3.1 Simulated data
The synthetic TCSPC dataset has an image size of 150 $\times$ 150 pixels and $M$ = 256. The simulated signals are bi-exponential ($P$ = 2). Figure 9(a) shows the $\log _{10}(I_i)$ image consisting of three regions with integrated intensities of $I_1$ = 500, $I_2$ =1000, and $I_3$ = 10000, respectively. Possion noise is included in the dataset. Figures 9(b), 9(c), and 9(d) shows the $q_1$, $\tau _1$, and $\tau _2$ images, respectively. The $q_1$ image has three regions with mean values of [0.2, 0.5, 0.8] and relative standard deviations of 10%; the $\tau _1$ image has two regions with mean values of [0.5, 1] ns and relative standard deviations of 10%; the $\tau _2$ image has a mean value of 3 ns and a relative standard deviation of 10%.
Fig. 9. (a) $\log _{10}(I_i)$, (b) $q_1$, (c) $\tau _1$, and (d) $\tau _2$ images of the synthetic TCSPC dataset.
The execution times ($t_{exe}$) and the mean squared errors (MSE) of the results evaluated by different algorithms without and with HC are summarized in Table 2. $MSE$ is defined as:
where $\hat {a}$ represents the estimated $a$ ($a = q_1, \tau _1, \tau _2, \tau _A, \tau _I$). The results for the fitting and non-fitting methods in PW and CW modes and the global fitting methods in GF-P and GF-C modes are illustrated and analyzed in the following sections.Table 2. $t_{exe}$ and $MSE$ evaluated by algorithms without and with HC.
3.1.1 Fitting method LE-LSM in PW and CW modes
Figure 10 shows $\hat {q}_1$, $\hat {\tau}_1$, $\hat {\tau}_2$, $\hat {\tau}_A$, and $\hat {\tau}_I$ images estimated with LE-LSM in (a1) – (a5) PW (without HC) and (b1) – (b5) CW (with HC) modes, respectively. The pixel brightness of each image represents the intensity. The parameter histograms ($\hat {q}_1$, $\hat {\tau}_1$, $\hat {\tau}_2$, $\hat {\tau}_A$, and $\hat {\tau}_I$) of different intensity regions (blue, red, and magenta lines for $I_1$, $I_2$, and $I_3$ respectively) and the true histogram (black dash line) are attached to each image. HC improves the estimated images, especially for Regions $I_1$ and $I_2$, as the histograms are closer to the truth than those estimated without HC. $MSE$s are reduced from 0.019 to 0.011 for $q_1$, from 0.173 ns$^2$ to 0.102 ns$^2$ for $\tau _1$, from 0.198 ns$^2$ to 0.104 ns$^2$ for $\tau _2$, from 0.110 ns$^2$ to 0.037 ns$^2$ for $\tau _A$, and from 0.027 ns$^2$ to 0.025 ns$^2$ for $\tau _I$. $t_{exe}$ is significantly reduced from 389.45 s to 3.36 s.
Fig. 10. $\hat {q}_1$, $\hat {\tau}_1$, $\hat {\tau}_2$, $\hat {\tau}_A$, and $\hat {\tau}_I$ images with LE-LSM in (a1) – (a5) PW (without HC) and (b1) – (b5) CW (with HC) modes. Histograms of different intensity regions (blue, red, and magenta for $I_1$, $I_2$, and $I_3$, respectively) and the true histogram (black dash line).
3.1.2 Non-fitting methods LE-IEM and CMM in PW and CW modes
Figure 11 shows $\hat {\tau}_A$ and $\hat {\tau}_I$ images and histograms produced by LE-IEM and CMM (a) – (b) without and (c) – (d) with HC. LE-IEM is for estimating $\hat {\tau}_A$. $MSE(\tau _A)$ is improved from 0.102 ns$^2$ to 0.038 ns$^2$, and $t_{exe}$ is reduced from 62.30 s to 0.63 s with HC.
Fig. 11. $\hat {\tau}_A$ and $\hat {\tau}_I$ images and histograms produced by LE-IEM and CMM (a) – (b) without and (c) – (d) with HC.
CMM is for estimating $\hat {\tau}_I$ with the shortest $t_{exe}$ either in PW or CW, around 0.20 s, but it has a bias, as shown in Figs. 11(b) and 11(d). $MSE(\tau _I)$ is around 0.180 ns$^2$. There is a way to correct the bias, as described in [46].
3.1.3 Global fitting methods IC and VP in GF-P and GF-C modes
Figure 12 shows $\hat {q}_1$, $\hat {\tau}_A$, and $\hat {\tau}_I$ images estimated with (a1) – (a3) IC and (a4) – (a6) VP in GF-P (without HC) mode and with (b1) – (b3) IC and (b4) – (b6) VP in GF-C (with HC) mode. The estimated constants ($\hat {\tau}_1$, $\hat {\tau}_2$) are labelled in corresponding $\hat {q}_1$ images. The estimations of $q_1$ in Region $I_1$ with IC in GF-P are mostly inaccurate, as shown in Fig. 12(a1) and ($\hat {\tau}_1$, $\hat {\tau}_2$) = (0.59, 2.74) ns. As a result, $\hat {\tau}_A$ and $\hat {\tau}_I$ are also not correct in Region $I_1$, as shown in Figs. 12(a2) and (a3). In GF-C, IC performs better with a successfully estimated Region $I_1$, a significantly reduced $t_{exe}$ from 724.82 s to 11.85 s, and a reduced $MSE(q_1)$ from 0.100 to 0.017, a reduced $MSE(\tau _A)$ from 0.678 to 0.102, and a reduced $MSE(\tau _I)$ from 1.098 to 0.050, as shown in Figs. 12(b1) - (b3).
Fig. 12. $\hat {q}_1$, $\hat {\tau}_A$, and $\hat {\tau}_I$ images and histograms with (a1) – (a3) IC and (a4) – (a6) VP in GF-P mode and with (b1) – (b3) IC and (b4) – (b6) VP in GF-C mode. Constants ($\hat {\tau}_1$, $\hat {\tau}_2$) are labelled in corresponding $\hat {q}_1$ images.
HC also accelerates VP from 3.34 s to 0.31 s with a reduced $MSE(q_1)$ from 0.033 to 0.014, a reduced $MSE(\tau _A)$ from 0.122 to 0.048, and a reduced $MSE(\tau _I)$ from 0.178 to 0.093, as shown in Figs. 12(b4) – (b6).
Although VP has some invalid estimations with $\hat {q}_1$ < 0 (pixels in white) when $q_1$ is small, as shown in Figs. 12(a4) and (b4), its $\hat {\tau}_A$ and $\hat {\tau}_I$ images are accurately evaluated without invalid pixels, as shown in Figs. 12(a5), (a6), (b5), and (b6). Thus, VP in GF-C is a promising choice for fast average lifetime estimations for its short execution time ($t_{exe}$ = 0.31 s).
In conclusion, HC not only accelerates analysis but also enhances accuracy (MSE). HC sorts histograms with similar profiles into a cluster and takes the average of histograms for lifetime determination, equivalent to increasing the number of photon counts and reducing noise. Therefore, the decay parameters estimated with the average cluster histogram by lifetime determination algorithms have higher accuracy than those of individual histograms. Although the decay parameters of the histograms in one cluster have a deviation from those estimated with the average cluster histogram, the results indicate that the error introduced by HC is smaller than that introduced by processing original histograms with a relatively lower photon count.
3.2 Experimental data
Mouse raw macrophage cells were routinely cultured in DMEM (Dulbecco’s Modified Eagle Medium) supplemented with 10$\%$ FCS (Fetal Calf Serum) under 5$\%$ CO$_2$ at 37$^oC$. Cells were seeded on glass cover slips in 24-well plates and cultured overnight for bacterial infection. Bacteria engineered to express GFP (Green Fluorescent Protein) were harvested from an early exponential phase and added to the cells with an MOI (Multiplicity of Infection) = 100. Cells were washed with PBS (Phosphate-Buffered Saline) and stained for actin with phalloidin Alexa Flour 546 (Thermo Fisher Scientific). The scanning FLIM used in this work is LSM510 (Carl Zeiss), equipped with a TCSPC module (SPC-830, Becker $\&$ Hickl GmbH). The sample was excited by a tunable femtosecond Ti: Sapphire laser (Chameleon, Coherent) at 850 nm as a two-photon excitation source. The repetition rate is 80 MHz, and the pulse width is less than 200 fs. The emitted photons were collected through a 63$\times$ water-immersion objective lens (N.A. = 1.0) and a 500 $\sim$ 550 nm bandpass filter and transferred into a photomultiplier tube.
Figure 13 shows the intensity image and $t_{exe}$ for all algorithms without or with HC for a TCSPC dataset. The estimations of three output types are shown as follows.
3.2.1 Type 1: $\hat {q}_1$, $\hat {\tau}_1$, and $\hat {\tau}_2$ images with LE-LSM in PW and CW modes
Figure 14 shows (a) $\hat {q}_1$, (b) $\hat {\tau}_1$, and (c) $\hat {\tau}_2$ images with LE-LSM in PW, (d) - (f) the results with LE-LSM in CW, and the histograms of (g) $\hat {q}_1$, (h) $\hat {\tau}_1$, and (i) $\hat {\tau}_2$ in PW (blue) and CW (red) modes. LE-LSM shows similar lifetime estimation performances in PW and CW, whereas LE-LSM in CW ($t_{exe}$ = 5.87 s) is faster than LE-LSM in PW ($t_{exe}$ = 632.32 s). Therefore, LE-LSM used in CW is a better choice for Type 1.
Fig. 14. $\hat {q}_1$, $\hat {\tau}_1$, and $\hat {\tau}_2$ images from LE-LSM (a) - (c) without and (d) – (f) with HC. Histograms of (g) $\hat {q}_1$, (h) $\hat {\tau}_1$, and (i) $\hat {\tau}_2$ in PW (blue) and CW (red) modes.
3.2.2 Type 2: $\hat {q}_1$ image and constants ($\hat {\tau}_1$, $\hat {\tau}_2$) with IC and VP in GF-P and GF-C modes
Figure 15 shows $\hat {q}_1$ images with (a) IC and (b) VP in GF-P, (c) IC and (d) VP in GF-C. Figure 15(e) shows histograms of $\hat {q}_1$ with IC (dash blue) and VP (dash red) in GF-P and IC (solid blue) and VP (solid red) in GF-C. The constants ($\hat {\tau}_1$, $\hat {\tau}_2$) of each approach are attached in $\hat {q}_1$ images.
Fig. 15. $\hat {q}_1$ images from IC and VP (a) – (b) without and (c) – (d) with HC. (e) Histograms of $\hat {q}_1$ with IC (dash blue) and VP (dash red) in GF-P and IC (solid blue) and VP (solid red) in GF-C.
3.2.3 Type 3: $\hat {\tau}_A$ and $\hat {\tau}_I$ images
Figure 16 shows $\hat {\tau}_A$ images (a) – (d) without and (e) – (h) with HC. Figure 16(i) shows the histograms of $\hat {\tau}_A$ with LE-LSM and LE-IEM in PW and CW. Figure 16(j) shows the histograms of $\hat {\tau}_A$ with IC and VP in GF-P and GF-C. Figure 17 shows $\hat {\tau}_I$ images (a) – (d) without and (e) – (h) with HC. Figure 17(i) shows the histograms of $\hat {\tau}_I$ with LE-LSM and CMM in PW and CW. Figure 17(j) shows the histograms of $\hat {\tau}_I$ with IC and VP in GF-P and GF-C.
Fig. 16. $\hat {\tau}_A$ images from the algorithms (a) – (d) without and (e) – (h) with HC. (i) Histograms of $\hat {\tau}_A$ with LE-LSM and LE-IEM in PW and CW. (j) histograms of $\hat {\tau}_A$ with IC and VP in GF-P and GF-C.
Fig. 17. $\hat {\tau}_I$ images from the algorithms (a) – (d) without and (e) – (h) with HC. (i) Histograms of $\hat {\tau}_I$ with LE-LSM and CMM in PW and CW. (j) histograms of $\hat {\tau}_I$ with IC and VP in GF-P and GF-C.
Like the conclusions drawn from simulations, LE-LSM in CW is the fastest for Type 1 with $t_{exe}$ = 5.87 s; VP in GF-C is the fastest for Type 2 with $t_{exe}$ = 0.41 s. For average lifetime images, VP in GF-C is the fastest for both $\tau _A$ and $\tau _I$ with $t_{exe}$ = 0.41 s, LE-IEM in CW is the second one for $\tau _A$ with $t_{exe}$ = 0.94 s; meanwhile, CMM in CW is the fastest for $\tau _I$ with $t_{exe}$ = 0.20 s.
4. Conclusion
We developed a histogram clustering (HC) method to accelerate FLIM analysis. HC can improve both the speed and the accuracy for FLIM analysis by sorting histograms with similar profiles in a dataset into several clusters and significantly reducing the number of histograms to be analyzed. The HC method implements clustering with two features of a histogram. Several commonly used lifetime determination algorithms’ performances for producing decay parameter images without and with HC were compared using synthetic and experimental datasets. For different output types, the fastest FLIM analysis methods are suggested: 1) LE-LSM with HC for all lifetime component images with an execution time ($t_{exe}$) of 5.87 s, 106-fold shorter than $t_{exe}$ without HC; 2) VP with HC for constant lifetimes, $q_1$, $\tau _A$, and $\tau _I$ images with $t_{exe}$ = 0.41 s, 32-fold shorter than $t_{exe}$ without HC; 3) LE-IEM with HC as the second choice for $\tau _A$ with $t_{exe}$ = 0.94 s, 78-fold shorter than $t_{exe}$ without HC, and CMM as the second choice for $\tau _I$ with $t_{exe}$ = 0.2 s without or with HC (biased if the largest lifetime > $T$/4). The analysis was conducted in Matlab, and it can be translated to C or other environments to speed up the analysis. We believe the proposed HC method can benefit applications demanding real-time FLIM such as clinical diagnosis and fast screening.
Funding
Science Foundation of the Chinese Academy of Sciences (CXJJ-21S006); Strategic Priority Research Program of Chinese Academy of Sciences (XDA25030900); National Natural Science Foundation of China (12075312); Medical Research Scotland (1179-2017); Engineering and Physical Sciences Research Council (EP/L01596X/1).
Acknowledgments
N. Sapermsap acknowledges the Development and Promotion of Science and Technology Talents (DPST) Project under the Institute for the Promotion of Teaching Science and Technology (IPST), Thailand, for a PhD scholarship.
Disclosures
The authors declare no conflicts of interest.
Data availability
Data underlying the results presented in this paper are not publicly available and can be made available by the authors without undue reservation.
Supplemental document
See Supplement 1 for supporting content.
References
1. J. R. Lakowicz, Principles of Fluorescence Spectroscopy (Springer, 2006).
2. E. F. Fornasiero, S. Mandad, H. Wildhagen, M. Alevra, B. Rammner, S. Keihani, F. Opazo, I. Urban, T. Ischebeck, M. S. Sakib, M. K. Fard, K. Kirli, T. P. Centeno, R. O. Vidal, R. U. Rahman, E. Benito, A. Fischer, S. Dennerlein, P. Rehling, I. Feussner, S. Bonn, M. Simons, H. Urlaub, and S. O. Rizzoli, “Precisely measured protein lifetimes in the mouse brain reveal differences across tissues and subcellular fractions,” Nat. Commun. 9(1), 4230 (2018). [CrossRef]
3. K. Okabe, N. Inada, C. Gota, Y. Harada, T. Funatsu, and S. Uchiyama, “Intracellular temperature mapping with a fluorescent polymeric thermometer and fluorescence lifetime imaging microscopy,” Nat. Commun. 3(1), 705 (2012). [CrossRef]
4. C. Ma, W. Sun, L. Xu, Y. Qian, J. Dai, G. Zhong, Y. Hou, J. Liu, and B. Shen, “A minireview of viscosity-sensitive fluorescent probes: design and biological applications,” J. Mater. Chem. B 8(42), 9642–9651 (2020). [CrossRef]
5. D. Aigner, R. I. Dmitriev, S. M. Borisov, D. B. Papkovsky, and I. Klimant, “pH-sensitive perylene bisimide probes for live cell fluorescence lifetime imaging,” J. Mater. Chem. B 2(39), 6792–6801 (2014). [CrossRef]
6. T. Gadella, FRET and FLIM Techniques (Elsevier B.V., 2009).
7. J. A. Levitt, S. P. Poland, N. Krstajic, K. Pfisterer, A. Erdogan, P. R. Barber, M. Parsons, R. K. Henderson, and S. M. Ameer-Beg, “Quantitative real-time imaging of intracellular FRET biosensor dynamics using rapid multi-beam confocal FLIM,” Sci. Rep. 10(1), 5146 (2020). [CrossRef]
8. W. Becker, “Fluorescence lifetime imaging - techniques and applications,” J. Microsc. 247(2), 119–136 (2012). [CrossRef]
9. K. Suhling, L. M. Hirvonen, J. A. Levitt, P. H. Chung, C. Tregidgo, A. Le Marois, D. A. Rusakov, K. Zheng, S. Ameer-Beg, S. Poland, S. Coelho, R. Henderson, and N. Krstajic, “Fluorescence lifetime imaging (FLIM): Basic concepts and some recent developments,” Medical Photonics 27, 3–40 (2015). [CrossRef]
10. F. Poulon, J. Pallud, P. Varlet, M. Zanello, F. Chretien, E. Dezamis, G. Abi-Lahoud, F. Nataf, B. Turak, B. Devaux, and D. Abi Haidar, “Real-time Brain Tumor imaging with endogenous fluorophores: a diagnosis proof-of-concept study on fresh human samples,” Sci. Rep. 8(1), 14888 (2018). [CrossRef]
11. N. Hirmiz, A. Tsikouras, E. J. Osterlund, M. Richards, D. W. Andrews, and Q. Fang, “Highly multiplexed confocal fluorescence lifetime microscope designed for screening applications,” IEEE J. Sel. Top. Quantum Electron. 27(5), 1–9 (2021). [CrossRef]
12. C. Garciá, A. Losada, M. A. Sacristán, J. F. Martínez-Leal, C. M. Galmarini, and M. P. Lillo, “Dynamic cellular maps of molecular species: Application to drug-target interactions,” Sci. Rep. 8(1), 1140 (2018). [CrossRef]
13. A. Ehn, J. Zhu, X. Li, and J. Kiefer, “Advanced laser-based techniques for gas-phase diagnostics in combustion and aerospace engineering,” Appl. Spectrosc. 71(3), 341–366 (2017). [CrossRef]
14. M. Jonsson, A. Ehn, M. Christensen, M. Aldén, and J. Bood, “Simultaneous one-dimensional fluorescence lifetime measurements of OH and CO in premixed flames,” Appl. Phys. B 115(1), 35–43 (2014). [CrossRef]
15. A. Ehn, O. Johansson, J. Bood, A. Arvidsson, B. Li, and M. Aldén, “Fluorescence lifetime imaging in a flame,” Proc. Combust. Inst. 33(1), 807–813 (2011). [CrossRef]
16. W. Becker, Advanced Time-Correlated Single Photon Counting Techniques (Springer, 2005).
17. W. Becker, “Fluorescence lifetime imaging by multi-dimensional time correlated single photon counting,” Medical Photonics 27, 41–61 (2015). [CrossRef]
18. L. M. Hirvonen, W. Becker, J. Milnes, T. Conneely, S. Smietana, A. Le Marois, O. Jagutzki, and K. Suhling, “Picosecond widefield time-correlated single photon counting fluorescence microscopy with a delay line anode detector,” Appl. Phys. Lett. 109(7), 071101 (2016). [CrossRef]
19. D. D.-U. Li, J. Arlt, J. Richardson, R. Walker, A. Buts, D. Stoppa, E. Charbon, and R. Henderson, “Real-time fluorescence lifetime imaging system with a 32 × 32 0.13μm CMOS low dark-count single-photon avalanche diode array,” Opt. Express 18(10), 10257 (2010). [CrossRef]
20. K. Dowling, S. Hyde, J. Dainty, P. French, and J. Hares, “2-D fluorescence lifetime imaging using a time-gated image intensifier,” Opt. Commun. 135(1-3), 27–31 (1997). [CrossRef]
21. S. Yu, T. Yao, and B. Yuan, “An ICCD camera-based time-domain ultrasound-switchable fluorescence imaging system,” Sci. Rep. 9(1), 10552 (2019). [CrossRef]
22. R. V. Krishnan, E. Biener, J. H. Zhang, R. Heckel, and B. Herman, “Probing subtle fluorescence dynamics in cellular proteins by streak camera based fluorescence lifetime imaging microscopy,” Appl. Phys. Lett. 83(22), 4658–4660 (2003). [CrossRef]
23. L. Camborde, A. Jauneau, C. Brière, L. Deslandes, B. Dumas, and E. Gaulin, “Detection of nucleic acid–protein interactions in plant leaves using fluorescence lifetime imaging microscopy,” Nat. Protoc. 12(9), 1933–1950 (2017). [CrossRef]
24. H. Chen and E. Gratton, “A practical implementation of multifrequency widefield frequency-domain fluorescence lifetime imaging microscopy,” Microsc. Res. Tech. 76(3), 282–289 (2013). [CrossRef]
25. P. J. Verveer and Q. S. Hanley, “Chapter 2 Frequency domain FLIM theory, instrumentation, and data analysis,” in Laboratory Techniques in Biochemistry and Molecular Biology, T. W. J. Gadella, ed. (Elesvier B.V., 2009), pp. 59–94.
26. W. R. Ware, L. J. Doemeny, and T. L. Nemzek, “Deconvolution of fluorescence and phosphorescence decay curves. A least-squares method,” J. Phys. Chem. 77(17), 2038–2048 (1973). [CrossRef]
27. V. V. Apanasovich and E. G. Novikov, “Deconvolution method for fluorescence decays,” Opt. Commun. 78(3-4), 279–282 (1990). [CrossRef]
28. J. Liu, Y. Sun, J. Qi, and L. Marcu, “A novel method for fast and robust estimation of fluorescence decay dynamics using constrained least-squares deconvolution with Laguerre expansion,” Physics in Medicine and Biology 57(4), 843–865 (2012). [CrossRef]
29. Y. Zhang, Y. Chen, and D. D.-U. Li, “Optimizing Laguerre expansion based deconvolution methods for analyzing bi-exponential fluorescence lifetime images,” Opt. Express 24(13), 13894 (2016). [CrossRef]
30. J. Jo, Q. Fang, and L. Marcu, “Ultrafast method for the analysis of fluorescence lifetime imaging microscopy data based on the Laguerre expansion technique,” IEEE J. Sel. Top. Quantum Electron. 11(4), 835–845 (2005). [CrossRef]
31. S. Pelet, M. Previte, L. Laiho, and P. C. So, “A fast global fitting algorithm for fluorescence lifetime imaging microscopy based on image segmentation,” Biophys. J. 87(4), 2807–2817 (2004). [CrossRef]
32. S. C. Warren, A. Margineanu, D. Alibhai, D. J. Kelly, C. Talbot, Y. Alexandrov, I. Munro, M. Katan, C. Dunsby, and P. M. W. French, “Rapid Global Fitting of Large Fluorescence Lifetime Imaging Microscopy Datasets,” PLoS One 8(8), e70687 (2013). [CrossRef]
33. J. Xu, A. Zhang, Z. Gao, K. Nie, and J. Qiao, “A low data-rate fluorescence lifetime imaging system with CMM pre-processed in-pixel,” Microelectron. J. 69, 28–34 (2017). [CrossRef]
34. D. D.-U. Li, B. Rae, R. Andrews, J. Arlt, and R. Henderson, “Hardware implementation algorithm and error analysis of high-speed fluorescence lifetime sensing systems using center-of-mass method,” J. Biomed. Opt. 15(1), 017006 (2010). [CrossRef]
35. S. P. Poland, A. T. Erdogan, N. Krstajić, J. Levitt, V. Devauges, R. J. Walker, D. D.-U. Li, S. M. Ameer-Beg, and R. K. Henderson, “New high-speed centre of mass method incorporating background subtraction for accurate determination of fluorescence lifetime,” Opt. Express 24(7), 6899–6915 (2016). [CrossRef]
36. D. D.-U. Li, E. Bonnist, D. Renshaw, and R. Henderson, “On-chip, time-correlated, fluorescence lifetime extraction algorithms and error analysis,” Journal of the Optical Society of America A 25(5), 1190 (2008). [CrossRef]
37. D. D.-U. Li, R. Walker, J. Richardson, B. Rae, A. Buts, D. Renshaw, and R. Henderson, “Hardware implementation and calibration of background noise for an integration-based fluorescence lifetime sensing algorithm,” Journal of the Optical Society of America A 26(4), 804–814 (2009). [CrossRef]
38. M. A. Digman, V. R. Caiolfa, M. Zamai, and E. Gratton, “The phasor approach to fluorescence lifetime imaging analysis,” Biophys. J. 94(2), L14–L16 (2008). [CrossRef]
39. S. Ranjit, L. Malacrida, D. M. Jameson, and E. Gratton, “Fit-free analysis of fluorescence lifetime imaging data using the phasor approach,” Nat. Protoc. 13(9), 1979–2004 (2018). [CrossRef]
40. B. B. Collier and M. J. McShane, “Dynamic windowing algorithm for the fast and accurate determination of luminescence lifetimes,” Anal. Chem. 84(11), 4725–4731 (2012). [CrossRef]
41. S. P. Chan, Z. J. Fuller, J. N. Demas, and B. A. DeGraff, “Optimized gating scheme for rapid lifetime determinations of single-exponential luminescence lifetimes,” Anal. Chem. 73(18), 4486–4490 (2001). [CrossRef]
42. Y. Li, S. Natakorn, Y. Chen, M. Safar, M. Cunningham, J. Tian, and D. D.-U. Li, “Investigations on average fluorescence lifetimes for visualizing multi-exponential decays,” Frontiers in Physics 8, 576862 (2020). [CrossRef]
43. T. A. Laurence and B. A. Chromy, “Efficient maximum likelihood estimator fitting of histograms,” Nat. Methods 7(5), 338–339 (2010). [CrossRef]
44. M. I. Rowley, A. C. C. Coolen, B. Vojnovic, and P. R. Barber, “Robust Bayesian fluorescence lifetime estimation, decay model selection and instrument response determination for low-intensity FLIM imaging,” PLoS One 11(6), e0158404 (2016). [CrossRef]
45. J. T. Smith, R. Yao, N. Sinsuebphon, A. Rudkouskaya, N. Un, J. Mazurkiewicz, M. Barroso, P. Yan, and X. Intes, “Fast fit-free analysis of fluorescence lifetime imaging via deep learning,” Proc. Natl. Acad. Sci. 116(48), 24019–24030 (2019). [CrossRef]
46. D. D.-U. Li, H. Yu, and Y. Chen, “Fast bi-exponential fluorescence lifetime imaging analysis methods,” Opt. Lett. 40(3), 336–339 (2015). [CrossRef]
