## 1.

## Introduction

Fluorescence microscopy has become an indispensable tool in the study of tissues and cells, as well as in material science due to its high degree of specificity amid nonfluorescing material. In particular, the presence of a wide range of endogenous fluorophores^{1} in biological compounds has made this technique suitable for optically probing its complex cellular components. The steady state fluorescence intensity method is commonly employed in such studies, but is limited to distinguishing fluorophores with nonoverlapping spectra.^{2} The fluorescence lifetime technique, on the other hand, is able to overcome this limitation as it discriminates distinct fluorophores by characterizing the temporal response of the fluorescence emission. It is typically used to investigate the radiative and nonradiative decay rates of fluorophores on optical excitation. This promising technique has been exploited to monitor the subtle perturbation in cellular environment^{3, 4} and to determine the extent of protein-protein interactions in Förster resonance energy transfer (FRET) assays.^{5} Application of the fluorescence lifetime technique in cancer diagnosis has also been reported.^{6, 7} In particular, a recent work demonstrated the ability of the lifetime technique in detecting a tumor-targeted fluorescent marker with two distinct forms inside the tumor cells.^{8} This additional information provides a better insight into the uptake mechanism of tumor cells, and could possibly lead to an improved formulation of photosensitizer for photodynamic therapy and diagnosis.

The fluorescence lifetime can be measured with either modulated excitations in the frequency domain (FD) or pulsed excitation in the time domain (TD). In principle, both methods are equivalent and are related to each other by the Fourier transform. A comparison of the two approaches has recently been made to quantify their accuracy and signal-to-noise ratio^{9, 10} (SNR). Time-correlated single-photon counting (TCSPC), a TD technique, has been shown to offer a better SNR compared to the FD technique.^{9} The well-defined statistics in the single-photon counting measurement is also an advantage for data analysis.^{10} Despite the important advantages of the TCSPC technique, a successful discrimination of multiple fluorophores remains a complex issue. The temporal response of fluorescence emission is typically modeled as a summation of the multiple exponential decay function with discrete lifetimes. The decay parameters, comprising the fluorescence lifetimes and fractional contributions of each lifetime component, are typically estimated using the nonlinear least-squares (NLLS) method.^{11} This parametric estimation method, however, has been shown to be incapable of determining the actual number of fluorophores present in a sample.^{12} The NLLS method has been found to erroneously overfit the experimental data with a higher lifetime component number.^{2, 13} Global analysis^{14} of NLLS-fitted decay profiles is also widely used to improve the discrimination of closely spaced lifetimes.^{2} Despite the improved discrimination, the uncertainties in the extracted lifetimes are still significant.^{2} In most practical applications, a biexponential decay model is assumed and its use is often justified by a least-squares deviation value of close to 1. A stretched exponential function has been proposed to represent the fluorescence decay profile in complex, heterogeneous biological samples with a single characteristic time constant and a heterogeneity parameter.^{13} These alternative parameters are shown to yield better contrast in the imaging of tissue samples exhibiting continuous lifetime distributions. However, the technique requires good-quality fluorescence decay measurements with low noise and sufficiently large sample points.^{6} Furthermore, the stretched exponential representation may not be suitable in applications requiring the accurate discrimination and determination of discrete lifetimes for analysis.^{5}

In this paper, we investigate the discrimination of discrete components in fluorescing samples using an expectation-maximization (EM) method with joint deconvolution to extract the fluorescence decay parameters in the multiexponential model. The optimal number of fluorescent components is selected using the Bayesian information criterion (BIC). Both simulation and TCSPC measurements were carried out to compare the performance between the NLLS and our proposed method. The results showed that a better lifetime prediction and discrimination were obtained with the EM-BIC method. In particular, the EM-BIC method is able to accurately identify the correct number of fluorescent components in samples where one or more components are weakly fluorescing.

## 2.

## Theory

## 2.1.

### NLLS Method

Fluorescence emission involves a transition from the lowest vibrational level of an excited singlet electronic state to any of the vibrational sublevels in the ground state.^{15} The emission rate of fluorescence decay is related to the rate of depopulation of the upper-state electrons. The fluorescence dynamics
$F\left(t\right)$
as a function of time
$t$
is conventionally described by the multiexponential equation^{16}

## Eq. 1

$F\left(t\right)=\sum _{j=1}^{m}{\alpha}_{j}\phantom{\rule{0.2em}{0ex}}\mathrm{exp}(-t\u2215{\tau}_{j}),$In the conventional NLLS method, the estimated parameters
${\alpha}_{k}$
and
${\tau}_{k}$
are iteratively adjusted until a best fit between the estimated decay
${R}_{c}\left({t}_{i}\right)$
and the experimental decay
$R\left({t}_{i}\right)$
is obtained. While a number of nonlinear minimization methods are available, the Marquartdt-Levenberg method is widely used as a benchmark for fitting comparison.^{16, 17} In the curve-fitting procedure, the reduced chi-squared error
${\chi}_{r}^{2}$
given by

## Eq. 3

${\chi}_{r}^{2}=\frac{1}{a-b}\sum _{i=1}^{n}\frac{{[R\left({t}_{i}\right)-{R}_{c}\left({t}_{i}\right)]}^{2}}{{\sigma}_{i}^{2}}$## 2.2.

### Proposed EM Algorithm with Joint Deconvolution

The EM algorithm is a maximum likelihood estimator that can be applied to parameter estimation of a mixture model.^{18, 19} The advantages of EM have been exploited in many different applications.^{20, 21, 22} It provides an iterative method to simplify the search for maximum likelihood via an expectation step and a maximization step. For our application, we define the probability density function of a normalized observed fluorescence process
$g\left(t\right)$
as

## Eq. 4

$g\left(t\right)=\sum _{j=1}^{m}\frac{{\beta}_{j}\phantom{\rule{0.2em}{0ex}}\mathrm{exp}(-t\u2215{\tau}_{j})}{{\tau}_{j}}\otimes h\left(t\right),$Using the EM algorithm, the $\mathbf{\beta}$ and $\mathbf{\tau}$ components of a multicomponent fluorescence decay profile with $N$ total photons can be determined based on Eqs. 6, 7, respectively:

## Eq. 7

${\tau}_{j}=\frac{1}{N{\beta}_{j}}\sum _{i=1}^{N}P\left({y}_{ij}\right|{t}_{i})({t}_{i}-\frac{{M}_{ij}}{{K}_{ij}}),$## Eq. 8

$P\left({y}_{ij}\right|{t}_{i})=\frac{({\beta}_{j}{K}_{ij}\u2215{\tau}_{j})\mathrm{exp}(-{t}_{i}\u2215{\tau}_{j})}{\sum _{r=1}^{m}({\beta}_{r}{K}_{ir}\u2215{\tau}_{r})\mathrm{exp}(-{t}_{i}\u2215{\tau}_{r})},$## Eq. 9

${K}_{ij}={\int}_{0}^{{t}_{i}}h\left({t}^{\prime}\right)\mathrm{exp}\left(\frac{{t}^{\prime}}{{\tau}_{j}}\right)\mathrm{d}{t}^{\prime},$## Eq. 10

${M}_{ij}={\int}_{0}^{{t}_{i}}{t}^{\prime}h\left({t}^{\prime}\right)\mathrm{exp}\left(\frac{{t}^{\prime}}{{\tau}_{j}}\right)\mathrm{d}{t}^{\prime}.$Following parameter estimation, BIC (Ref. 23) is used to determine the optimal number of parameters to represent the given fluorescence decay profile. Its purpose is to reasonably grade the goodness of fit by penalizing the log-likelihood function with the model complexity. The BIC parameter $\lambda $ for the multiexponential model is given by

where## 3.

## Materials and Methods

## 3.1.

### Simulation

To compare the performance of the conventional NLLS and our proposed methods, theoretical lifetime decay profiles with different total photon count, lifetime combination, and contribution factor were generated. A Monte Carlo model was used to simulate biexponential decay profiles in MATLAB^{™}. The instrumental noise, including the electronic jitter of the TCSPC system and temporal width of laser pulses, was derived from the probability density function of a measured IRF in our measurement system and added to the simulated decay profiles. All simulated decay profiles consist of 3750 data points with temporal resolution of
$0.01\phantom{\rule{0.3em}{0ex}}\mathrm{ns}$
, giving a time span of
$37.5\phantom{\rule{0.3em}{0ex}}\mathrm{ns}$
.

Biexponential decay profiles with widely separated lifetimes of ${\tau}_{1}=1\phantom{\rule{0.3em}{0ex}}\mathrm{ns}$ and ${\tau}_{2}=4\phantom{\rule{0.3em}{0ex}}\mathrm{ns}$ were first generated at varying contribution factors ( ${\beta}_{1}=0.1$ , 0.3, 0.5, 0.7, and 0.9) for a total photon count of ${10}^{5}$ . To evaluate the performance of the methods on closely spaced lifetime parameters, biexponential decay profiles with ${\tau}_{1}=2\phantom{\rule{0.3em}{0ex}}\mathrm{ns}$ and ${\tau}_{2}=4\phantom{\rule{0.3em}{0ex}}\mathrm{ns}$ were generated at varying contribution factors ( ${\beta}_{1}=0.1$ , 0.5, and 0.9) for a total photon count of ${10}^{5}$ . Both methods were also evaluated using decay profiles with low photon counts. The biexponential decays with ${\tau}_{1}=1\phantom{\rule{0.3em}{0ex}}\mathrm{ns}$ , ${\tau}_{2}=4\phantom{\rule{0.3em}{0ex}}\mathrm{ns}$ , and ${\beta}_{1}=0.9$ were simulated for total photon counts ranging from $5\times {10}^{3}$ to ${10}^{5}$ . In all cases, 10 decay profiles were simulated for each parameter set. To quantify the accuracy of model selection and extracted lifetime parameters for each method, the 10 simulated decay profiles with random noise characteristics from each parameter set were fitted with the NLLS and EM-BIC methods. The model for each simulated decay profile is chosen based on the selection criterion in each method, i.e., a ${\chi}_{r}^{2}$ value closest to unity for the NLLS method and the largest BIC parameter $\lambda $ for the EM-BIC method. The model selection for a simulated decay profile is accurate if the chosen value of $m$ matches the true number of lifetime components in the parameter set. The model selection accuracy of each method is determined by taking the ratio of the number of simulated decay profiles that has its model correctly chosen to the total number of simulated decay profiles in a parameter set.

The NLLS method was implemented in MATLAB^{™} using the optimization function based on the Marquartdt-Levenberg method. The convergence criterion for the NLLS fitting is chosen to be a least-squares change of less than
${10}^{-8}$
. The EM method was also implemented in MATLAB^{™} and the convergence criterion was set to a change in the BIC parameter of less than
${10}^{-8}$
.

## 3.2.

### Experimental Technique

Experimental measurement of the fluorescence lifetime was carried out to verify the simulation study. Two commonly used fluorescent probes for biological studies, fluorescein and acridine orange, were selected owing to their highly overlapped emission spectrum from $500\phantom{\rule{0.3em}{0ex}}\text{to}\phantom{\rule{0.3em}{0ex}}650\phantom{\rule{0.3em}{0ex}}\mathrm{nm}$ . Fluorescein and acridine orange from Sigma Aldrich were dissolved in phosphate-buffered saline (PBS) to stabilize their pH at 7.4. For the measurement of fluorescence lifetime, fluorescein and acridine orange solutions were prepared at concentration of $0.01\phantom{\rule{0.3em}{0ex}}\mathrm{mM}$ and $50\phantom{\rule{0.3em}{0ex}}\mu \mathrm{M}$ , respectively.

A pulsed laser diode (LDH400, Picoquant) with an excitation wavelength centered at $400\phantom{\rule{0.3em}{0ex}}\mathrm{nm}$ was used to excite the fluorophore solutions. The laser diode was biased and pulsed at $20\phantom{\rule{0.3em}{0ex}}\mathrm{MHz}$ using a laser driver (PDL 800-B, Picoquant). The period of the optical pulse was selected such that it was sufficiently longer than the measured range of fluorescence lifetimes. To minimize the detection of spontaneous emission from the light source, the laser pulses were first transmitted through a $400\text{-}\mathrm{nm}$ bandpass filter before illuminating the samples.

The fluorophore solutions in quartz cuvettes were placed in a right-angled geometry for the lifetime measurement. Fluorescence emission from the solution is directed to a monochromator and detected by a microchannel plate photomultiplier tube (PMT; R3809U-51, Hamamatsu). The fluorescence emission emanating from fluorescein and acridine orange were measured at a wavelength of
$520\phantom{\rule{0.3em}{0ex}}\mathrm{nm}$
. The output current pulses of the PMT were amplified by a wideband preamplifier and converted into fast logic pulses before feeding into a constant fraction discriminator (CFD). The delayed reference pulses from the laser driver and the signal pulses from the CFD were then fed into a time analyzer (PTA 9308, Ortec) to record the photon arrival time. In TCSPC measurement, a single fluorescent photon is detected for every pulsed excitation. By repeating the single-photon-counting process a sufficiently large number of times, a histogram of photon arrival times was constructed to reveal the characteristic exponential decay of the fluorescence emission.^{24} The IRF of the measurement system was determined by measuring the excitation photon at
$400\phantom{\rule{0.3em}{0ex}}\mathrm{nm}$
from a scattering suspension of Ludox colloidal silica (Sigma-Aldrich).

## 4.

## Results

## 4.1.

### Simulation

The same initial guess values of the decay parameters $\mathbf{\beta}$ and $\mathbf{\tau}$ were used for both methods. All simulated decay profiles were fitted for lifetime component number $m$ of 1 to 4. In the NLLS method, the simulated biexponential decay profiles were fitted to Eq. 2 and the respective contribution factors were calculated based on Eq. 5. The optimal lifetime component number is obtained when ${\chi}_{r}^{2}$ is closest to unity. For the proposed EM-BIC method, Eqs. 6, 7, 8, 9, 10, 11, 12 were used to deduce the decay parameters. The best-fitting results are selected based on the largest $\lambda $ value. The typical results of the fitting procedure and model selection using the NLLS and EM-BIC methods are illustrated in Tables 1, 2 , respectively. The simulated decay profile consists of $1\times {10}^{5}$ photons originating from two fluorescent components ( ${\tau}_{1}=1\phantom{\rule{0.3em}{0ex}}\mathrm{ns}$ and ${\tau}_{2}=4\phantom{\rule{0.3em}{0ex}}\mathrm{ns}$ ) with equal contribution factors. In this example, the lifetime component number was erroneously determined as $m=4$ by the NLLS method, while the EM-BIC method correctly determines the lifetime component number.

## Table 1

Fluorescence lifetime parameters extracted from the biexponential decay profile simulated with τ1=1ns , τ2=4ns , and β1=0.5 using conventional NLLS method.

M | χr2 | β1 | β2 | β3 | β4 | τ1 (ns) | τ2 (ns) | τ3 (ns) | τ4 (ns) |
---|---|---|---|---|---|---|---|---|---|

1 | 4.2213 | 1.00 | — | — | — | 1.54 | — | — | — |

2 | 0.9699 | 0.51 | 0.49 | — | — | 1.02 | 4.09 | — | — |

3 | 0.9710 | 0.26 | 0.26 | 0.48 | — | 1.01 | 1.03 | 4.10 | — |

4 | 0.9717 | 0.26 | 0.26 | 0.24 | 0.24 | 1.02 | 1.02 | 4.06 | 4.13 |

## Table 2

Fluorescence lifetime parameters extracted from the bi-exponential decay profile simulated with τ1=1ns , τ2=4ns , and β1=0.5 using the EM-BIC method.

m | BIC Parameter, λ | β1 | β2 | β3 | β4 | τ1 (ns) | τ2 (ns) | τ3 (ns) | τ4 (ns) |
---|---|---|---|---|---|---|---|---|---|

1 | $-6.62143\times {10}^{5}$ | 1.00 | — | — | — | 2.38 | — | — | — |

2 | $-\mathbf{6.57859}\times {\mathbf{10}}^{\mathbf{6}}$ | 0.50 | 0.50 | — | — | 0.99 | 4.01 | — | — |

3 | $-6.57871\times {10}^{6}$ | 0.25 | 0.25 | 0.50 | — | 0.90 | 1.10 | 4.01 | — |

4 | $-6.57883\times {10}^{6}$ | 0.25 | 0.25 | 0.25 | 0.25 | 0.90 | 1.09 | 3.09 | 4.11 |

In the first study, we investigate the performance of the EM-BIC method on the simulated decay profiles comprising two fluorescent components with widely separated fluorescence lifetimes ( ${\tau}_{1}=1\phantom{\rule{0.3em}{0ex}}\mathrm{ns}$ and ${\tau}_{2}=4\phantom{\rule{0.3em}{0ex}}\mathrm{ns}$ ) and having a sufficiently large photon counts of $1\times {10}^{5}$ . Using these data, the two methods were compared for values of ${\beta}_{1}$ ranging from 0.1 to 0.9. The parameters predicted at a lifetime component number of $m=2$ are shown in Tables 3, 4 . Ten decay profiles were generated for each parameter set and fitted with each method to determine the corresponding mean and standard deviation of the extracted lifetimes. The model selection accuracy for each method is calculated from the proportion of the 10 simulated decay profiles for each parameter set that has its component number $m$ correctly determined by the model. The biexponential lifetime values and model selection accuracy determined from both methods are also plotted in Fig. 1 for comparison.

## Table 3

Biexponential lifetime parameters extracted from decay profiles simulated with widely separated lifetimes ( τ1=1ns , τ2=4ns , and 1×105 photons) using the conventional NLLS method.

β1 | τ1 (ns) | τ2 (ns) | Model Selection Accuracy (%) |
---|---|---|---|

0.1 | $1.02\pm 0.31$ | $3.96\pm 0.05$ | 20 |

0.3 | $0.99\pm 0.03$ | $4.00\pm 0.06$ | 10 |

0.5 | $1.00\pm 0.02$ | $4.03\pm 0.08$ | 30 |

0.7 | $1.00\pm 0.02$ | $4.01\pm 0.10$ | 50 |

0.9 | $1.00\pm 0.01$ | $4.26\pm 0.42$ | 60 |

## Table 4

Biexponential lifetime parameters extracted from decay profiles simulated with widely separated lifetimes ( τ1=1ns , τ2=4ns , and 1×105 photons) using the EM-BIC method.

β1 | τ1 (ns) | τ2 (ns) | Model selection accuracy (%) |
---|---|---|---|

0.1 | $1.02\pm 0.01$ | $4.00\pm 0.01$ | 100 |

0.3 | $0.99\pm 0.03$ | $4.01\pm 0.02$ | 100 |

0.5 | $0.99\pm 0.01$ | $4.03\pm 0.01$ | 100 |

0.7 | $0.99\pm 0.01$ | $4.00\pm 0.04$ | 100 |

0.9 | $1.00\pm 0.01$ | $4.00\pm 0.14$ | 100 |

The results show that the simulated decay profiles were fitted equally well with $m=2$ using both methods. However, there is a noticeable spread in the value of ${\tau}_{1}$ extracted with the NLLS method for ${\beta}_{1}=0.1$ . Similarly, we can see that the standard deviation of the estimated ${\tau}_{2}$ increases with ${\beta}_{1}$ for both fitting methods. This is attributed to the reduced photon number from the lifetime component whose fractional contribution is decreased, which resulted in a poorer estimate of the respective lifetime value. In particular, we note from Fig. 1 that the deviation is less significant when the EM-BIC method is used. As seen in Fig. 1a, the conventional NLLS method exhibits a decreasing accuracy in model selection for low values of ${\beta}_{1}$ even though the fluorescence lifetimes were relatively well estimated at $m=2$ . The NLLS method only correctly identifies the lifetime component number in 1 out of the 10 decay profiles (10%) simulated for the parameter set with ${\beta}_{1}=0.3$ . By contrast the EM-BIC method succeeds in selecting the correct lifetime component number for all the simulated decay profiles, as depicted in Fig. 1b. The results clearly indicate that the EM-BIC method outperforms the conventional NLLS method in determining the correct number of fluorescent components.

The next simulation study compares the two methods using decay profiles simulated with closely spaced lifetimes ( ${\tau}_{1}=2\phantom{\rule{0.3em}{0ex}}\mathrm{ns}$ and ${\tau}_{2}=4\phantom{\rule{0.3em}{0ex}}\mathrm{ns}$ ), a sufficiently large photon count of $1\times {10}^{5}$ and ${\beta}_{1}$ values of 0.1, 0.5, and 0.9. The lifetime values and model selection accuracy determined from the conventional NLLS method and our EM-BIC method are given in Tables 5, 6 , respectively. The results are also plotted in Fig. 2 for comparison. The average biexponential lifetimes predicted by both methods are in good agreement with the true values. However, the extracted lifetime value of the component that has a reduced photon count, i.e., a lower $\beta $ value, exhibits an increased standard deviation. The model selection accuracy also drops drastically with the use of the ${\chi}_{r}^{2}$ to resolve the two closely separated fluorescence lifetime components, with a maximum accuracy of only 40% for ${\beta}_{1}=0.1$ [Fig. 2a].

## Table 5

Fluorescence lifetime parameters extracted from the decay profiles simulated with closely spaced lifetimes ( τ1=2ns , τ2=4ns , and 1×105 photons) using the conventional NLLS method.

β1 | τ1 (ns) | τ2 (ns) | Model Selection Accuracy (%) |
---|---|---|---|

0.1 | $1.96\pm 0.33$ | $4.00\pm 0.05$ | 40 |

0.5 | $2.01\pm 0.03$ | $4.03\pm 0.06$ | 20 |

0.9 | $2.00\pm 0.02$ | $4.08\pm 0.09$ | 10 |

## Table 6

Fluorescence lifetime parameters extracted from decay profiles simulated with closely spaced lifetimes ( τ1=2ns , τ2=4ns , and 1×105 photons) using the EM-BIC method.

β1 | τ1 (ns) | τ2 (ns) | Model Selection Accuracy (%) |
---|---|---|---|

0.1 | $1.99\pm 0.07$ | $4.00\pm 0.02$ | 100 |

0.5 | $1.99\pm 0.02$ | $4.01\pm 0.03$ | 100 |

0.9 | $1.99\pm 0.01$ | $4.10\pm 0.04$ | 100 |

Last, we compare the performance of the conventional NLLS and our EM-BIC methods with varying total photon count in the simulated biexponential decay profile. Widely separated lifetime components with ${\tau}_{1}=1\phantom{\rule{0.3em}{0ex}}\mathrm{ns}$ and ${\tau}_{2}=4\phantom{\rule{0.3em}{0ex}}\mathrm{ns}$ were used to generate the decay profiles for various total photon count. The value of ${\beta}_{1}$ is set at 0.9 such that the simulated decay profiles are dominated by the shorter lifetime component. As shown in Tables 7, 8 , the shorter lifetime component ${\tau}_{1}$ was accurately extracted by both methods even at a very low photon count of $9\times {10}^{3}$ . By contrast, the extracted value of ${\tau}_{2}$ is less accurate and deviates significantly from the true value when the NLLS method is used. As depicted in Fig. 3a, the NLLS method gave a poor estimation of ${\tau}_{2}$ ranging from $3.01\phantom{\rule{0.3em}{0ex}}\text{to}\phantom{\rule{0.3em}{0ex}}4.73\phantom{\rule{0.3em}{0ex}}\mathrm{ns}$ for a total photon count of $1\times {10}^{4}$ . When the total photon count reduces to $5\times {10}^{3}$ , the uncertainty in ${\tau}_{2}$ becomes even larger and the estimate of ${\tau}_{2}$ ranges between 1.95 and $6.13\phantom{\rule{0.3em}{0ex}}\mathrm{ns}$ . At the lowest total photon count, the NLLS method was also unable to correctly determine the lifetime component number in any of the simulated decay profiles. By comparison, the EM-BIC method was able to determine the correct lifetime component number and the extracted lifetime values are reasonably close to the true values [see Table 8 and Fig. 3b].

## Table 7

Fluorescence lifetime parameters extracted from decay profiles simulated with varying total photon count ( τ1=1ns , τ2=4ns , and β1=0.9 ) using the conventional NLLS method.

Total PhotonCount (×104) | τ1 (ns) | τ2 (ns) | Model Selection Accuracy (%) |
---|---|---|---|

0.5 | $0.96\pm 0.08$ | $4.04\pm 2.09$ | 0 |

1 | $1.00\pm 0.03$ | $3.87\pm 0.86$ | 0 |

2 | $0.99\pm 0.03$ | $3.71\pm 0.73$ | 10 |

5 | $1.01\pm 0.02$ | $4.36\pm 0.40$ | 50 |

10 | $1.00\pm 0.01$ | $4.22\pm 0.33$ | 60 |

## Table 8

Fluorescence lifetime parameters extracted from decay profiles simulated with varying total photon count ( τ1=1ns , τ2=4ns , and β1=0.9 ) using the EM-BIC method.

Total PhotonCount (×104) | τ1 (ns) | τ2 (ns) | Model Selection Accuracy (%) |
---|---|---|---|

0.5 | $0.97\pm 0.02$ | $3.90\pm 0.44$ | 100 |

1 | $0.99\pm 0.02$ | $4.02\pm 0.21$ | 100 |

2 | $1.00\pm 0.01$ | $4.01\pm 0.19$ | 100 |

5 | $1.01\pm 0.01$ | $4.06\pm 0.11$ | 100 |

10 | $1.00\pm 0.01$ | $4.05\pm 0.16$ | 100 |

## 4.2.

### Experiment

The measured IRF and lifetime decay profiles of fluorescein and acridine orange for a total photon count of
$2\times {10}^{5}$
are shown in Fig. 4
. The full width at half maximum (FWHM) of the measured IRF is
$86\phantom{\rule{0.3em}{0ex}}\mathrm{ps}$
. We can see from Fig. 4 that the fluorescence lifetime of fluorescein is longer than that of acridine orange. Using the NLLS method, the fluorescence lifetime of acridine orange is determined to be
$1.84\phantom{\rule{0.3em}{0ex}}\mathrm{ns}$
with a
${\chi}_{r}^{2}$
of 1.03, while the estimated lifetime of fluorescein is
$4.18\phantom{\rule{0.3em}{0ex}}\mathrm{ns}$
with a
${\chi}_{r}^{2}$
of 1.09. The extracted lifetime values are in good agreement with the values reported in the literature^{25} and those from EM-BIC analysis.

Three biexponential decay profiles (A to C) with varying contributions from fluorescein and acridine orange were chosen to compare the performance and accuracy of the NLLS and EM-BIC methods. The respective contribution factors from each component and the total photon count are listed in Table 9 . The first biexponential decay profile A comprises fluorescence measured from the two fluorophores with equal contribution factors. Decay profile B is dominated by fluorescein, whereas decay profile C has a dominant contribution from acridine orange.

## Table 9

Fractional contribution and total photon count of the experimental biexponential decay profiles of fluorescein and acridine orange.

Decay profile | βacridineorange | βfluorescein | Total Photon Count (×104) |
---|---|---|---|

A | $1\u22152$ | $1\u22152$ | 40 |

B | $1\u221511$ | $10\u221511$ | 5.5 |

C | $10\u221511$ | $1\u221511$ | 5.5 |

The same initial guess values of the decay parameters $\mathbf{\beta}$ and $\mathbf{\tau}$ were used for both methods. All decay profiles were evaluated with each method using lifetime component numbers 1 to 3. The results of the parameter extraction and model selection using the NLLS and EM-BIC methods are shown in Tables 10, 11 , respectively. A comparison between the true lifetime parameters of each experimental decay profile and those determined from both methods with $m=2$ are also shown in Fig. 5 .

## Table 10

Fluorescence lifetime parameters extracted from the experimental biexponential decay profiles using the conventional NLLS method.

Decay Profile | m | χr2 | β1 | β2 | β3 | τ1 (ns) | τ2 (ns) | τ3 (ns) |
---|---|---|---|---|---|---|---|---|

1 | 3.6831 | 1.00 | — | — | 2.57 | — | — | |

A | 2 | 1.1064 | 0.50 | 0.50 | — | 1.85 | 4.22 | — |

3 | 1.1307 | 0.27 | 0.26 | 0.47 | 1.88 | 1.90 | 4.37 | |

1 | 0.9767 | 1.00 | — | — | 3.83 | — | — | |

B | 2 | 0.9945 | 0.12 | 0.88 | — | 2.30 | 4.19 | — |

3 | 0.9961 | 0.14 | 0.14 | 0.72 | 2.78 | 2.78 | 4.47 | |

1 | 1.1381 | 1.00 | — | — | 1.92 | — | — | |

C | 2 | 1.1307 | 0.73 | 0.27 | — | 1.68 | 3.21 | — |

3 | 1.1462 | 0.09 | 0.46 | 0.46 | 0.90 | 2.11 | 2.11 |

## Table 11

Fluorescence lifetime parameters extracted from the experimental biexponential decay profiles using the EM-BIC method.

Decay Profile | m | BIC Parameter, λ | β1 | β2 | β3 | τ1 (ns) | τ2 (ns) | τ3 (ns) |
---|---|---|---|---|---|---|---|---|

1 | $-2.70588\times {10}^{6}$ | 1.00 | — | — | 3.03 | — | — | |

A | 2 | $-\mathbf{2.70188}\times {\mathbf{10}}^{\mathbf{6}}$ | 0.49 | 0.51 | — | 1.80 | 4.26 | — |

3 | $-2.70189\times {10}^{6}$ | 0.25 | 0.25 | 0.51 | 1.75 | 1.85 | 4.27 | |

1 | $-3.8594\times {10}^{5}$ | 1.00 | — | — | 3.98 | — | — | |

B | 2 | $-\mathbf{3.85917}\times {\mathbf{10}}^{\mathbf{5}}$ | 0.10 | 0.90 | — | 1.88 | 4.21 | — |

3 | $-3.85928\times {10}^{5}$ | 0.05 | 0.05 | 0.90 | 1.86 | 1.94 | 4.21 | |

1 | $-3.52785\times {10}^{5}$ | 1.00 | — | — | 2.07 | — | — | |

C | 2 | $-\mathbf{3.52469}\times {\mathbf{10}}^{\mathbf{5}}$ | 0.89 | 0.11 | — | 1.80 | 4.66 | — |

3 | $-3.52479\times {10}^{5}$ | 0.05 | 0.05 | 0.89 | 4.60 | 4.72 | 1.80 |

The fitting results of profile A for both methods predicted a lifetime component number of $m=2$ , indicating the presence of two fluorescent components. The extracted lifetimes from both methods are in excellent agreement with the true lifetime parameters of profile A [Fig. 5a]. The goodness of fit for both methods becomes degraded at $m=3$ , although the corresponding fluorescence lifetimes are still reasonably close to the true values. For decay profile B, a lifetime component number of $m=3$ is predicted with the NLLS method. The value of ${\tau}_{1}(={\tau}_{2}=2.78\phantom{\rule{0.3em}{0ex}}\mathrm{ns})$ extracted from the NLLS method is obviously different from the true values. In contrast, the EM-BIC method predicted the correct lifetime component number of 2, with extracted lifetime values that are in good agreement with the true values, as shown in Fig. 5b. The NLLS analysis of profile C underestimated the lifetime values despite selecting the correct lifetime component number. On the other hand, the EM-BIC method accurately predicted the correct values for profile C [Fig. 5c].

## 5.

## Discussion

A comparison of the conventional NLLS method and our EM-BIC method using simulated decay profiles showed that the NLLS method is unable to accurately determine the number of fluorescent components and their corresponding lifetime values. The EM-BIC method consistently outperforms the NLLS method, particularly when the simulated biexponential decay profile contains a weakly contributing component or has a low total photon count.

Experimental measurements using two fluorophores to generate biexponential decay profiles to validate our EM-BIC method gave similar results as the simulation study. The experimental comparison indicates that the NLLS approach fails when there is a weakly fluorescing component that contributes a low photon count to the decay profile, i.e., profiles B and C. On the other hand, the EM-BIC method gave a more reliable lifetime prediction and discrimination even at low photon counts. The EM-BIC method can improve the accuracy of lifetime discrimination in fluorescence lifetime microscopy^{26, 27} to give better-quality lifetime images.

These results suggest that a unity-approaching
${\chi}_{r}^{2}$
in the NLLS method does not necessarily predict the correct lifetime component number and lifetime values. The failure of the NLLS method to discriminate the lifetime components is attributed to the poor model selection criterion with
${\chi}_{r}^{2}$
when the photon count is low. This can be explained by the assumptions of the NLLS method. In TCSPC experiments, the photon statistics are described by the Poisson distribution and can be approximated to the Gaussian distribution when the photon count is sufficiently large. The validity of the NLLS method holds only if the photon count in each time bin is independent and follows the Gaussian distribution.^{2, 28} A weakly fluorescing component results in a low photon contribution to the decay profile. Consequently the photon statistics are no longer Gaussian in the multiexponential decay profile. This scenario can commonly appear if the weakly fluorescing component has a long fluorescence lifetime. On the other hand, the EM-BIC method uses the correct photon statistics and thus has a more reliable lifetime prediction and discrimination even at low photon counts.

## 6.

## Conclusion

In this paper, an EM method with joint deconvolution was formulated to estimate the lifetime parameters, and the BIC was used to determine the number of fluorescing components in a lifetime decay profile. The EM-BIC method was shown to be more reliable in predicting and discriminating the fluorescence lifetime components in multiexponential fluorescence decay profiles. Unlike the conventional NLLS method, the EM-BIC method is able to correctly determine the number of lifetime components in the decay profiles and accurately extract the corresponding lifetime values, even when the photon count in the decay profile is low. This was attributed to the use of a more accurate description of the photon statistics in the EM-BIC method. The results indicate that the EM-BIC model is particularly useful for accurate fluorescence lifetime determination in samples with multiple fluorophores that are weakly fluorescing or lifetime measurements with a low photon count.

## Appendices

### Appendix

In the EM algorithm, the probability density function of a normalized observed fluorescence process $g\left(t\right)$ can be described as

## Eq. 13

$g\left(t\right)=\sum _{j=1}^{m}\frac{{\beta}_{j}\phantom{\rule{0.2em}{0ex}}\mathrm{exp}(-t\u2215{\tau}_{j})}{{\tau}_{j}}\otimes h\left(t\right)=\sum _{j=1}^{m}\frac{{\beta}_{j}{L}_{j}(t,{\tau}_{j})\mathrm{exp}(-t\u2215{\tau}_{j})}{{\tau}_{j}},$## Eq. 14

${L}_{j}(t,{\tau}_{j})={\int}_{0}^{t}h\left(\alpha \right)\mathrm{exp}\left(\frac{\alpha}{{\tau}_{j}}\right)\mathrm{d}\alpha .$## Eq. 15

$p\left(\mathbf{T}\right|\mathbf{\theta})=\prod _{i=1}^{n}\prod _{j=1}^{m}{\left[\frac{{\beta}_{j}{K}_{ij}}{{\tau}_{j}}\phantom{\rule{0.2em}{0ex}}\mathrm{exp}\left(\frac{-{t}_{i}}{{\tau}_{j}}\right)\right]}^{{y}_{ij}}$The aim of the EM algorithm is to simplify the maximization of the likelihood function via the expectation step and the maximization step. In the expectation step, the $Q$ function is defined as the conditional mean of log-likelihood function as follows:

## Eq. 16

$Q\left({\mathbf{\theta}}_{k}\right|{\mathbf{\theta}}_{k-1})=E\left[p\left(\mathbf{T}\right|{\mathbf{\theta}}_{k})\right|\mathbf{t},{\mathbf{\theta}}_{k-1}]=\sum _{i=1}^{n}\sum _{j=1}^{m}P\left({y}_{ij}\right|{t}_{i}),$Here $P\left({y}_{ij}\right|{t}_{i})$ is the probability that ${t}_{i}$ is generated by the $j$ ’th component and is given by

## Eq. 17

$P\left({y}_{ij}\right|{t}_{i})=\frac{({\beta}_{j}{K}_{ij}\u2215{\tau}_{j})\mathrm{exp}(-{t}_{i}\u2215{\tau}_{j})}{\sum _{r=1}^{m}({\beta}_{r}{K}_{ir}\u2215{\tau}_{r})\mathrm{exp}(-{t}_{i}\u2215{\tau}_{r})},$## Eq. 18

${K}_{ij}={\int}_{0}^{{t}_{i}}h\left({t}^{\prime}\right)\mathrm{exp}\left(\frac{{t}^{\prime}}{{\tau}_{j}}\right)\mathrm{d}{t}^{\prime}.$## Eq. 20

${\tau}_{j}=\frac{1}{n{\beta}_{j}}\sum _{i=1}^{n}P\left({y}_{ij}\right|{t}_{i})({t}_{i}-\frac{{M}_{ij}}{{K}_{ij}}),$## Eq. 21

${M}_{ij}={\int}_{0}^{{t}_{i}}{t}^{\prime}h\left({t}^{\prime}\right)\mathrm{exp}\left(\frac{{t}^{\prime}}{{\tau}_{j}}\right)\mathrm{d}{t}^{\prime}.$## Acknowledgments

This research was supported in part by a research grant from the Singapore Cancer Syndicate (Grant No. SCS-BU0052) under the Agency for Science, Technology and Research $\left({\mathrm{A}}^{*}\text{STAR}\right)$ and a manpower cofunding (Grant No. M48040149) from the Nanyang Technological University, Singapore.

## References

**,” Photochem. Photobiol. Sci., 3 (1), 127 –131 (2004). https://doi.org/10.1039/b306129a 1474-905X Google Scholar**

*Fluorescence lifetime imaging (FLIM) of rhodamine 123 in living cells***,” Opt. Express, 14 (10), 4412 –4426 (2006). https://doi.org/10.1364/OE.14.004412 1094-4087 Google Scholar**

*Time-resolved optical imaging provides a molecular snapshot of altered metabolic function in living human cancer cell models***,” Biosens. Bioelectron., 24 (3), 383 –390 (2008). https://doi.org/10.1016/j.bios.2008.04.015 0956-5663 Google Scholar**

*Time-resolved Forster-resonance-energytransfer DNA assay on an active CMOS microarray***,” J. Pathol., 199 (3), 309 –317 (2003). https://doi.org/10.1002/path.1286 0022-3417 Google Scholar**

*Flourescence lifetime imaging of unstained tissues: early results in human breast cancer***,” J. Biomed. Opt., 11 (2), 021004 (2006). https://doi.org/10.1117/1.2186045 1083-3668 Google Scholar**

*Laguerre-based method for analysis of time-resolved fluorescence data: application to in-vivo characterization and diagnosis of atherosclerotic lesions***,” Photochem. Photobiol., 76 (6), 686 –694 (2002). https://doi.org/10.1562/0031-8655(2002)076<0686:IOAPBT>2.0.CO;2 0031-8655 Google Scholar**

*Internalization of aggregated photosensitizers by tumor cells: Subcellular time-resolved fluorescence spectroscopy on derivatives of pyropheophorbide-a ethers and chlorin e6 under femtosecond one- and two-photon excitation***,” J. Biomed. Opt., 8 (3), 381 –390 (2003). https://doi.org/10.1117/1.1586704 1083-3668 Google Scholar**

*Fluorescence lifetime imaging for the two-photon microscope: time-domain and frequency-domain methods***,” Anal. Chem., 79 (5), 2137 –2149 (2007). https://doi.org/10.1021/ac062160k 0003-2700 Google Scholar**

*Fluorescence lifetime standards for time and frequency domain fluorescence spectroscopy***,” Anal. Biochem., 59 (2), 583 –598 (1974). https://doi.org/10.1016/0003-2697(74)90312-1 0003-2697 Google Scholar**

*On the analysis of fluorescence decay kinetics by the method of least-squares***,” Chem. Phys. Lett., 120 (4–5), 455 –459 (1985). https://doi.org/10.1016/0009-2614(85)85640-2 0009-2614 Google Scholar**

*A fallacy in the interpretation of fluorescence decay parameters***,” Biophys. J., 81 (3), 1265 –1274 (2001). https://doi.org/10.1016/S0006-3495(01)75784-0 0006-3495 Google Scholar**

*Application of the stretched exponential function to fluorescence lifetime imaging***,” Chem. Phys. Lett., 120 (4–5), 466 –472 (1985). https://doi.org/10.1016/0009-2614(85)85642-6 0009-2614 Google Scholar**

*Global analysis of fluorescence decay surfaces: excited-state reactions***,” J. Membr. Biol., 61 378 –425 (1979). 0022-2631 Google Scholar**

*Time-resolved fluorescence measurements***,” Opt. Commun., 134 (1–6), 371 –378 (1997). https://doi.org/10.1016/S0030-4018(96)00384-7 0030-4018 Google Scholar**

*Fast fitting of multi-exponential decay curves***,” J. R. Stat. Soc. Ser. B (Methodol.), 39 (1), 1 –38 (1977). 0035-9246 Google Scholar**

*Maximum likelihood from incomplete data via the EM algorithm***,” 753 –760 (2005). Google Scholar**

*Expectation maximization algorithms for conditional likelihoods***,” 198 –207 (2004). Google Scholar**

*Recursive EM algorithm for finite mixture models with application to Internet traffic modeling***,” IEEE Trans. Signal Process., 54 (4), 1289 –1303 (2006). https://doi.org/10.1109/TSP.2006.870586 1053-587X Google Scholar**

*Unsupervised learning of parsimonious mixtures on large spaces with integrated feature and component selection***,” J. Phys. Chem. B, 109 (1), 617 –628 (2005). https://doi.org/10.1021/jp0467548 1089-5647 Google Scholar**

*Detection of intensity change points in time-resolved single-molecule measurements***,” Curr. Sci., 80 (9), 1135 –1144 (2001). 0011-3891 Google Scholar**

*Model selection—an overview***,” Proc. SPIE, 4259 102 –109 (2001). https://doi.org/10.1117/12.432487 0277-786X Google Scholar**

*Time-domain measurement of fluorescence lifetime variation with pH***,” J. Phys. D: Appl. Phys., 35 (9), R61 –R76 (2002). https://doi.org/10.1088/0022-3727/35/9/201 0022-3727 Google Scholar**

*Time-resolved fluorescence imaging in biology and medicine***,” Photochem. Photobiol. Sci., 3 (8), 795 –801 (2004). https://doi.org/10.1039/b316456j 1474-905X Google Scholar**

*Time domain fluorescence lifetime imaging applied to biological tissue***,” Anal. Biochem., 206 (2), 215 –225 (1992). https://doi.org/10.1016/0003-2697(92)90356-C 0003-2697 Google Scholar**

*Why, when, and how biochemists should use least squares*