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 fluorophores1 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 environment3, 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 ratio9, 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 analysis14 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.
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 as a function of time is conventionally described by the multiexponential equation16is the number of fluorescence lifetime components, and and are the decay amplitude and fluorescence lifetime of the ’th component, respectively. The observed fluorescence decay contains the instrumental effect and is given by the convolution of with the instrumental response function (IRF) written as
In the conventional NLLS method, the estimated parameters and are iteratively adjusted until a best fit between the estimated decay and the experimental decay 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 given byis the number of data points, is the number of estimated parameters, and is the standard deviation of the photon count in the ’th time bin. In TCSPC measurements, the photon statistics follow the Poisson distribution, and is approximated by the square root of . The value is expected to approach unity when the estimated function matches closely with the data.
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 asis unity. Here of the ’th component is calculated from the and vector components using following expression:
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 for the multiexponential model is given by
Materials and Methods
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 , giving a time span of .
Biexponential decay profiles with widely separated lifetimes of and were first generated at varying contribution factors ( , 0.3, 0.5, 0.7, and 0.9) for a total photon count of . To evaluate the performance of the methods on closely spaced lifetime parameters, biexponential decay profiles with and were generated at varying contribution factors ( , 0.5, and 0.9) for a total photon count of . Both methods were also evaluated using decay profiles with low photon counts. The biexponential decays with , , and were simulated for total photon counts ranging from to . 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 value closest to unity for the NLLS method and the largest BIC parameter for the EM-BIC method. The model selection for a simulated decay profile is accurate if the chosen value of 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 . The EM method was also implemented in MATLAB™ and the convergence criterion was set to a change in the BIC parameter of less than .
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 . 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 and , respectively.
A pulsed laser diode (LDH400, Picoquant) with an excitation wavelength centered at was used to excite the fluorophore solutions. The laser diode was biased and pulsed at 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 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 . 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 from a scattering suspension of Ludox colloidal silica (Sigma-Aldrich).
The same initial guess values of the decay parameters and were used for both methods. All simulated decay profiles were fitted for lifetime component number 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 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 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 photons originating from two fluorescent components ( and ) with equal contribution factors. In this example, the lifetime component number was erroneously determined as by the NLLS method, while the EM-BIC method correctly determines the lifetime component number.
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)|
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)|
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 ( and ) and having a sufficiently large photon counts of . Using these data, the two methods were compared for values of ranging from 0.1 to 0.9. The parameters predicted at a lifetime component number of 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 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.
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 (%)|
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 (%)|
The results show that the simulated decay profiles were fitted equally well with using both methods. However, there is a noticeable spread in the value of extracted with the NLLS method for . Similarly, we can see that the standard deviation of the estimated increases with 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 even though the fluorescence lifetimes were relatively well estimated at . 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 . 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 ( and ), a sufficiently large photon count of and 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 value, exhibits an increased standard deviation. The model selection accuracy also drops drastically with the use of the to resolve the two closely separated fluorescence lifetime components, with a maximum accuracy of only 40% for [Fig. 2a].
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 (%)|
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 (%)|
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 and were used to generate the decay profiles for various total photon count. The value of 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 was accurately extracted by both methods even at a very low photon count of . By contrast, the extracted value of 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 ranging from for a total photon count of . When the total photon count reduces to , the uncertainty in becomes even larger and the estimate of ranges between 1.95 and . 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].
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 (%)|
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 (%)|
The measured IRF and lifetime decay profiles of fluorescein and acridine orange for a total photon count of are shown in Fig. 4 . The full width at half maximum (FWHM) of the measured IRF is . 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 with a of 1.03, while the estimated lifetime of fluorescein is with a of 1.09. The extracted lifetime values are in good agreement with the values reported in the literature25 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.
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)|
The same initial guess values of the decay parameters and 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 are also shown in Fig. 5 .
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)|
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)|
The fitting results of profile A for both methods predicted a lifetime component number of , 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 , although the corresponding fluorescence lifetimes are still reasonably close to the true values. For decay profile B, a lifetime component number of is predicted with the NLLS method. The value of 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].
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 microscopy26, 27 to give better-quality lifetime images.
These results suggest that a unity-approaching 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 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.
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.
In the EM algorithm, the probability density function of a normalized observed fluorescence process can be described asbe an indicator variable defined as follows:13 as is the complete data, is the parameter set, and is the total number of recorded photons.
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 function is defined as the conditional mean of log-likelihood function as follows:and are the current and previously estimated parameters, respectively.
Here is the probability that is generated by the ’th component and is given byand by maximizing the derived function subjected to the constraint . With the use of Lagrange multipliers, the solutions are and are assumed to be constants, where is taken as the old estimate. This approximation is valid since . These expectation and maximization steps are then repeated until convergence is achieved.
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 and a manpower cofunding (Grant No. M48040149) from the Nanyang Technological University, Singapore.