Investigations on Average Fluorescence Lifetimes for Visualizing Multi-Exponential Decays

Intensity- and amplitude-weighted average lifetimes, denoted as τI and τA hereafter, are useful indicators for revealing Förster resonance energy transfer (FRET) or fluorescence quenching behaviors. In this work, we discussed the differences between τI and τA and presented several model-free lifetime determination algorithms (LDA), including the center-of-mass, phasor, and integral equation methods for fast τI and τA estimations. For model-based LDAs, we discussed the model-mismatch problems, and the results suggest that a bi-exponential model can well approximate a signal following a multi-exponential model. Depending on the application requirements, suggestions about the LDAs to be used are given. The instrument responses of the imaging systems were included in the analysis. We explained why only using the τI model for FRET analysis can be misleading; both τI and τA models should be considered. We also proposed using τA/τI as a new indicator on two-photon fluorescence lifetime images, and the results show that τA/τI is an intuitive tool for visualizing multi-exponential decays.


INTRODUCTION
Fluorescence lifetime imaging (FLIM) is a crucial technique for assessing microenvironments of fluorescent molecules [1,2], such as pH [3], Ca 2+ [4,5], O 2 [6], viscosity [7], or temperature [8]. Combining with Förster Resonance Energy Transfer (FRET) techniques, FLIM can be a powerful "quantum ruler" to measure protein conformations and interactions [9][10][11][12]. Compared with fluorescence intensity imaging, FLIM is independent of the signal intensity and fluorophore concentrations, making FLIM a powerful quantitative imaging technique for applications in life sciences [13], medical diagnosis [14][15][16], drug developments [17][18][19], and flow diagnosis [20][21][22]. FLIM techniques can build on time-correlated single-photon counting (TCSPC) [23][24][25], time-gating [26][27][28], or streak cameras [29]; they record time-resolved fluorescence intensity profiles to extract lifetimes with a lifetime determination algorithm (LDA) [1]. There is a rapid growth of real-time applications that fast analysis is sought after [12,30]. Traditional LDAs usually use the least square method (LSM) or maximum likelihood estimation (MLE) [31] to analyze decay models chosen by users, and model-fitting analysis follows a reduced chi-squared criterion [1]. In reality, however, it is difficult to know the exact decay model as fluorescent molecules in biological systems can demonstrate complex multi-exponential decay profiles. For instance, a mixture of fluorophores, a multi-tryptophan protein, single fluorophores in varied environments, and single-tryptophan proteins in multiple conformational states [1] can show multiexponential decays as where A represents the amplitude, q i and τ i (i = 1, . . . , p) denote the amplitude fractions and lifetimes, respectively, and p is the number of lifetime components. There are time-domain or frequency-domain [32][33][34][35] FLIM systems to measure a fluorescence decay. In this work, we focus on timedomain approaches. Suppose the instrument response function (IRF) of the measurement system is irf (t), the task performed by FLIM analysis tools is to extract f (t) from the measured decay h(t), as The problems with traditional LSM or MLE are two-fold. (1) It is challenging to categorize a fluorescence emission into a specific exponential model described by Equation (1) in complex biological processes. An arbitrary choice of p in Equation (1) simply based on reduced chi-squared tests [36] would lead to totally different interpretations. As the fitting routine is not mathematically unique; a measured decay could be fitted equally well with a bi-exponential or a tri-exponential model.
(2) To ensure the accuracy, it usually needs a high photon count (long acquisition time) when p ≥ 2 [37]. Instead of completely extracting q i and τ i (i = 1, . . . , p), which is doubtful as mentioned above and time-consuming, in many applications, it is often useful to determine only the average lifetime which can be expressed in two forms [1]: the intensity-weighted average lifetime and the amplitude-weighted average lifetime The question about which average lifetime we should use according to the applications has been investigated in [38]. For instance, they suggested: (a) τ A can estimate the energy transfer efficiency in FRET [39], where E is the energy transfer efficiency, I DA and I D are the fluorescence intensities of the donor in the presence and absence of energy transfer, respectively, and τ DA,A and τ D,A are τ A of the donor in the presence and absence of energy transfer, respectively. E can further estimate the donor-acceptor distance.
(b) τ A can also assess dynamic quenching behaviors, described by the Stern-Volmer equation [40], where I 0 and I 1 are fluorescence intensities, τ 0,A and τ 1,A are τ A of the fluorophore in the absence and presence of the quencher, respectively, K D is the Stern-Volmer quenching constant, and [Q] is the concentration of the quencher. Additionally, the average radiative rate constant can be expressed as, k r = QE/τ 1,A , where QE is the quantum yield.
(c) τ I can be used to estimate the average collisional constant k q from the Stern-Volmer constant K D .
In this work, we theoretically investigated two types of average lifetimes evaluated by model-free LDAs, examined the performances of τ A and τ I estimations using different LDAs, and suggested the choices of LDAs in terms of accuracy, precision, and estimation speeds according to the applications. We also described a multi-exponential decay visualization tool using the ratio τ A /τ I . Experimental results demonstrate the performance of τ A /τ I in comparison with Phasor.

THEORY
In this section, we derived the average lifetimes determined by the model-free methods, CMM, Phasor, and IEM and described the general work flow of average lifetime estimations with the model-free and model-based LDAs.
As Equation (2), the measured signal h(t) is the convolution of f (t) with irf (t). Here we focus on the signal h m and irf m obtained from a TCSPC system, as shown in Figure 1, where h m is the photon count collected in Bin m at t m = (m + 1/2) t, M is the number of bins, and t is the time resolution.
(a) CMM Frontiers in Physics | www.frontiersin.org The average lifetime evaluated with CMM is which is equal to τ I . The derivation of Equation (8) is shown in the Appendix.

(b) Phasor
The average lifetime evaluated with Phasor is where ω = 2π/T, T = M t is the measurement window, and g and s are the phasor components expressed as τ P is a weighted average lifetime whose weights are q i τ i /(1 + ω 2 τ 2 i ). If τ i ≪ T, then the weights are approximately equal to q i τ i , i.e., τ P is close to τ I .

(c) IEM
For IEM, the underlying exponential decay should be extracted by a model-free deconvolution method. With the estimated exponential decayf m , the average lifetime with IEM is where S m = [1/3, 4/3, 2/3, 4/3, 1/3] are the coefficients for numerical integration based on Simpson's rule. τ IEM is actually an estimator for τ A . Figure 2 summarizes the flow diagram for τ I and τ A estimations with different algorithms used in this study. The simulated signals h m and irf m are directly sent into CMM and Phasor blocks to estimate τ I . The estimated f m (from h m and irf m with the Laguerre expansion deconvolution method with L = 16 and α = 0.912 [54,55]) is sent to IEM to estimate τ A and sent to the bi-decay center-of-mass method (BCMM; j = 2) [56], the variable projection method (VPM; j = 2) [57], or LSM with a j-exponential model (denoted as LSM-j), to estimate τ I and τ A . CMM and Phasor are fast as no deconvolution routine is needed, whereas IEM, BCMM, VPM, and LSM are direct or iterative estimation approaches once f m is extracted. Artificial neural network assisted analysis tools [58,59] can be included in this diagram, but they are out of the scope of this work.

Simulations
In reality, it is difficult to characterize a real fluorescence profile with a suitable exponential model described in Equation (1). To demonstrate how model-free analysis can be beneficial, we examined two scenarios. Case A: we used exponential decay signals with p = 1 ∼ 4 to assess the influence of the model mismatch between the signal and the algorithm on τ I and τ A estimations. This study is to investigate the scenario when users select a j-exponential model to analyze a p-exponential decay (p can be different from j). Case B: we generated synthetic biexponential (p = 2) decay signals to assess the performances of τ I and τ A estimations with the model-free and model-based LDAs.
The performances of lifetime estimations can be assessed in two aspects: (1) the accuracy B n = |τ n − τ n |/τ n and the precision F n = √ N tot στ n /τ n [60], where n = I or A for the intensity-or the amplitude-weighted lifetimes, τ n andτ n are actual and estimated values, στ n is the standard deviation ofτ n , and N tot is the total photon count. The lower the F, the higher the precision (F = 1 for the ideal case).

Case A: Model Mismatch
Ideally, a bi-exponential signal should be analyzed by a biexponential model. For instance, BCMM, VPM, and LSM-2 are used for bi-exponential decay models, and LSM-j for jexponential models, j > 2. However, in realistic biological processes, it is difficult to know precisely how many lifetime components a decay profile contains. In traditional FLIM analysis tools, users usually need to select an exponential model to fit measured decays and use the reduced chi-squared to evaluate the goodness-of-fit. If the reduced chi-squared is not satisfactory, then a different exponential model is chosen. This process continues until the reduced chi-squared is acceptable. Often different exponential models can produce similar reduced chisquared values, and the question is which fitting we should use? It is quite common that a j-exponential model might analyze a signal containing p lifetime components and j = p. We would like to know if p is unknown to the user, whether using a different analysis model (j = p) would lead to a different biological story. We generated exponential decay signals h m (m = 0,. . ., M-1) to test the LDAs for τ I and τ A estimations. h m can be artificially generated with 3,4], q i = 1/p, and the IRF is approximated with a Poisson distribution irf m = exp(λ)λ m /m! with λ = 500 ps, FWHM ≃ 300 ps, and M = 256. The measurement window T = 10 ns, and the total photon count N tot = 10 3 . τ i , τ I , and τ A for each p are summarized in Table 1.
The performances of τ I and τ A estimations with the simulated exponential decays are shown in Figures 3A-D Figures 3E-G, respectively. The χ 2 values are insensitive to p for the three LDAs so that we conclude that a bi-exponential decay is suitable to approximate an arbitrary p-exponential decay (p ≤ 4).
Therefore, if the decay model of the signal is inaccessible, model-free and model-based LDAs, BCMM, VPM, and LSM-2 are enough for τ I and τ A estimations.
In practice, users can choose an optimization algorithm and set initial conditions to analyze FLIM images when LSM-2 is used. We would like to know how they can affect τ I and τ A estimations. Four bi-exponential decays, Decays 1 ∼ 4, with different parameters (q 1 , τ 1 , τ 2 ) were analyzed using LSM-2 with different initial conditions (q 10 , τ 10 , τ 20 ), denoted as Init. 1 ∼ 4 listed in Table 2 with N tot = 10 3 . When either of the estimated τ 1 and τ 2 is larger than T (10 ns), we say that the estimation fails. The probabilities of producing a failed trial, P(τ 1 or τ 2 > 10) and producing biased τ I and τ A with B I and B A > 0.3, i.e., P(B n > 0.3), n = I or A, are shown in Figure 4. Figures 4A-F are the LSM-2 results with the unconstrained and constrained trustregion-reflective (TRR) algorithms, respectively. The constraints are 0 < q 1 < 1 and 0 < τ 1 , τ 2 < 10 ns. Figures 4G-I are the LSM-2 results using the Levenberg-Marquardt (LM) algorithm. For the unconstrained TRR, the performances are relatively sensitive to initial conditions. P(τ 1 or τ 2 > 10) for Init. 4 is quite significant which results in large P(B n > 0.3), n = I  or A, for all four decays. Although Init. 3 leads to a low P for Decays 2 ∼ 3, P(B A > 0.3) for Decay 1 rises to 0.7. Thus, if the initial conditions are not chosen properly, the quality of τ I and τ A images cannot be guaranteed. The constrained TRR and LM are insensitive to initial conditions. Although the LM has failed trials, they barely affect P(B n > 0.3), n = I or A. Therefore, to ensure accurate τ I and τ A estimations, the constrained TRR and LM are recommended for LSM-2.

Case B: Performances of Average Lifetime Estimations
As mentioned above, it might be challenging to use a proper exponential model to describe realistic biological processes; a biexponential model might well approximate them. Here we will use a bi-exponential model to explain why model-free LDAs have the benefits of higher photon efficiency and faster analysis than model-based LDAs for τ I and τ A estimations. h m can be artificially generated with the same IRF used in Case A and f (t) = A[q 1 exp(−t/τ 1 ) + (1 − q 1 )exp(−t/τ 2 )], where τ 1 < τ 2 and q 1 is the amplitude fraction of τ 1 . Figure 5A shows the signal and IRF. In FRET and dynamic quenching applications, the fluorescence lifetime of the donor fluorophore is in general decreasing, and we assume τ 2 = 2.5 ns and τ 1 varying from 0.1 to 2.5 ns to emulate FRET or quenching. The theoretical τ I , τ A , and τ P with q 1 = 0.5 are shown in Figure 5B. τ P has a negative bias from τ I . With T/τ 2 increasing, τ P approaches τ I . Figure 5B that two different (τ 1 , τ 2 ) sets can deliver the same τ I , for instance, (0.32, 2.5) ns and (2.1, 2.5) ns have the same τ I of 2.3 ns.
Therefore, only estimating τ I can be misleading. Figure 5B also shows that the dynamic range of τ I is only 2.5-2.23 = 0.27 ns and within which the above problem persists. Whereas τ A does not have this problem for this case. We conducted Monte Carlo simulations to estimate τ I and τ A with the simulated signals, including Poisson noise under different conditions q 1 = 0.2, 0.5, and 0.8.
The performances of τ I and τ A estimations with biexponential decay signals are shown in Figures 6A-D for B I , F I , B A , and F A , respectively. Forτ I , B I,CMM , and B I,BCMM are roughly 10 and 8%, respectively determined by T/τ 2 . The larger T/τ 2 is, the smaller B I becomes (with F I,CMM and F I,BCMM being closer to 1). Phasor has a lower accuracy when q 1 becomes larger and τ 1 smaller, and it is less precise than CMM. VPM and LSM-2 both have a smaller B I = 3% but higher F I (1.5 ∼ 5) than CMM and BCMM. Forτ A , B A is 7% except for τ 1 = 0.1 ns, and F A is around 5 for the four LDAs. Figures 6C,D show that if only τ A is needed, there is no need to resort to slower model-based LDAs.
For τ I estimations, LSM-2 and VPM are preferred when high accuracy is required. Still, they are slower and have lower photon efficiency than CMM and BCMM which means the photon count should be higher to have similar precision, for instance, a relative standard deviation of 5% can be reached with N tot = 3,600 for LSM-2 and N tot = 500 for CMM and BCMM. When the  accuracy of CMM or BCMM (10% @ T/τ 2 = 4) is acceptable, CMM or BCMM should be employed for their high photon efficiency and estimation speeds. CMM is faster than BCMM as it can work without deconvolution. For τ A estimations, since the performances of IEM, BCMM, VPM, LSM-2 are similar, IEM can be the right candidate for fast analysis. Notice that the τ A method is less photon efficient than the τ I method as F A is higher than F I .

Experimental Results
tSA201 cells, which are a transformed human kidney cell line, were co-transfected with hP2Y 12 -eCFP and hP2Y 1 -eYFP receptors. After 48 h of transfection, the cells on the coverslips were washed once gently with PBS followed by fixation with ice-cold methanol for 10 min at room temperature. After being washed three times with PBS, they were mounted on to glass microscope slides with Mowiol. The microscope slides were then stored in the dark at room temperature overnight to allow the coverslips to dry, then stored at 4 • C for later use.
Cells were imaged on LSM510 (Carl Zeiss) equipped with a TCSPC module (SPC-830, Becker & Hickl GmbH), to determine the fluorescence lifetime and consequently the amount of FRET. The donor is CFP with the excitation wavelength range of 350 ∼ 500 nm and the emission wavelength range of 450 ∼ 600 nm. The acceptor is YFP. The sample is scanned pixel by pixel by a femtosecond Ti:Sapphire laser (Chameleon, Coherent) with an average output laser power of 3.8 W at 800 nm, as a two-photon excitation source to reduce cellular damage. The laser power is controlled with two polarizers. The repetition rate is 80 MHz with illuminating duration <200 fs. The emitted fluorescence signal from the donor is collected through a 63× water-immersion objective lens (N.A. = 1.0), a 480 ∼ 520 nm bandpass filter, and transferred into a photomultiplier tube (PMT) detector. The FLIM scanning was performed in a dark room containing the microscope. A set of experimental data (256 × 256 pixels, M = 256, T = 10 ns) was collected over an exposure period of up to 15 min. The IRF is obtained from the measurement of dried urea [(NH 2 ) 2 CO] [61].

Average Lifetime Images With LSM, CMM, and IEM
Figures 7A,C show the τ I and τ A images of the data evaluated by LSM-2 with an execution time of 410 s. The lifetime images were evaluated on Matlab R2016a, 64-bit with the Intel(R) Celeron(R) CPU (2950M @ 2 GHz) with 20923 pixels above an intensity threshold. Figures 7B,D show the τ I and τ A images evaluated by CMM and IEM with execution times of 0.25 and 92.3 s, respectively. IEM can be further accelerated to 0.6 s per image with histogram classification methods (we will report the details soon), as shown in Figure 7E. Although Fast-IEM causes a small bias in some pixels, the mean square error is acceptable with 0.005 ns 2 . The color bar represents lifetimes and the pixel brightness represents photon counts. The Figures 7F,G are histograms of τ I and τ A , respectively. Although the histogram of τ I with CMM deviates slightly from the one with LSM-2, CMM is 1,800-fold faster than LSM-2. If T/τ i > 4, the bias of τ I with CMM would become smaller. The τ A images are almost the same with IEM and LSM-2, whereas IEM and Fast-IEM are much faster than LSM-2.
Since the FRET efficiency E has a linear relationship with the average lifetimes as shown in Equation (5), Figures 7A-E can also be used to represent E images with the color bar in the range of 0 ∼ 100%. As we mentioned in Introduction, it is straightforward to obtain E images from τ A images, so that Figures 7C-E are proper E images. If τ I images are misused for E, the results would be different, as shown in Figures 7A,B, leading to a different biological story.

Visualization of Multi-Exponential Decays
With τ A /τ I τ I and τ A can not only access the essential parameters in FRET and dynamic quenching processes but also indicate the positions where multi-exponential decays occur. As mentioned previously, a fluorescence signal can be approximated by a bi-exponential decay, so that the ratio of τ I and τ A can be expressed as where R = τ 1 /τ 2 . The distribution of τ A /τ I (Figure 8) shows that when R ≃ 1 or q 1 ≃ 0 or 1, τ A /τ I ≃ 1. With a decrease of R or an increase of q 1 , τ A /τ I decreases. Therefore, the ranges of q 1 and R of a pixel can be determined by τ A /τ I . To present the multi-exponential decay visualization performance of τ A /τ I , the τ I and τ A images evaluated by LSM-2,  Figures 9A,B, were used to generate the τ A /τ I image as shown in Figure 9C. The histograms of τ I and τ A and the phasor plot are shown in Figures 9D,E. Figure 9F shows the possible range of q 1 and R of the selected pixels in Figure 9C. Figures 9C,F share the same color bar. Figure 9D shows that τ A has a broader lifetime dynamic range than τ I , which is consistent with the theoretical lines shown in Figure 5B. The τ A histogram shows two clusters with different peaks, whereas the τ I histogram only indicates a single merged group, meaning that there is no way to differentiate these two clusters. It is why using τ I to analyze samples with a strong FRET can be misleading.
The results of the selected pixels within different τ A /τ I ranges are shown in Figure 10, τ A /τ I = 0.2 ∼ 0.5, and Figure 11, τ A /τ I = 0.5 ∼ 1. For the pixels with τ A /τ I = 0.2 ∼ 0.5, the histograms clearly show that τ A is much smaller than τ I , which means the difference between τ 1 and τ 2 is significant. Figure 10F shows that the ranges of q 1 and R are approximately 0.5 ∼ 1 and 0 ∼ 0.2, respectively. For the pixels with τ A /τ I = 0.5 ∼ 1, τ A is closer to τ I , meaning the pixels have decays close to mono-exponential. 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 show the q 1 and R ranges more intuitively than phasor plots. τ A /τ I can be a useful tool to visualize the properties of the fluorescence decays within a lifetime image.

DISCUSSION
In realistic samples, fluorescence signals always follow multiexponential decay models. However, extracting lifetime components with a traditional fitting method is a timeconsuming process. For some applications that require calculating FRET efficiency and accessing dynamic quenching behaviors, average lifetimes are satisfactory. Model-free lifetime determination algorithms can be used to evaluate average lifetimes directly, for instance, CMM and Phasor for intensityweighted average lifetimes τ I and IEM for amplitude-weighted average lifetimes τ A . Discussions of the influence of the model mismatch between the real signal and the model-based LDAs on τ I and τ A estimations suggest that a bi-exponential model can well-approximate a signal following a multiple-exponential model. The results of the Monte-Carlo simulations suggest that VPM and LSM based on a bi-exponential model can be used for applications requiring high accuracy. The constrained TRR and LM algorithms with proper initial conditions are supported for LSM to guarantee accuracy. In contrast, CMM and IEM are recommended for applications requiring high estimation speeds.   We also explained why τ I models can be misleading, and τ I and τ A models should be considered. Experimental data were used to compare the performances of LSM-2, CMM, and IEM for evaluating τ I and τ A images. Similar τ I and τ A images were generated, whereas CMM and IEM are much faster than LSM-2. The data were further analyzed with τ A /τ I , which is capable of indicating the possible ranges of the amplitude proportion of the short lifetime and the ratio of the short and long lifetimes. We believe τ A /τ I is a useful and intuitive tool for visualizing multi-exponential decays in a lifetime image.

DATA AVAILABILITY STATEMENT
The raw data supporting the conclusions of this article will be made available by the authors, without undue reservation.

AUTHOR CONTRIBUTIONS
YL conducted theoretical and experimental analysis and developed analysis tools. MS and MC conceived FRET experiments and prepared samples. SN and YC contributed to FRET-FLIM experiments. JT contributed to tool developments. DL initiated the research concept and supervised the project. All authors wrote and revised the paper.

APPENDIX
Derivation of τ CMM . Take the integration of t · h(t) and h(t), Dividing Equation (A1) by Equation (A2) gives Then,