Fluorescence lifetime discrimination using expectation-maximization algorithm with joint deconvolution

C. Y. Fu; Beng-Koon Ng; Sirajudeen Gulam Razul

doi:10.1117/1.3258835

1 November 2009 Fluorescence lifetime discrimination using expectation-maximization algorithm with joint deconvolution

C. Y. Fu, Beng-Koon Ng, Sirajudeen Gulam Razul

Author Affiliations +

Journal of Biomedical Optics, Vol. 14, Issue 6, 064009 (November 2009). https://doi.org/10.1117/1.3258835

Abstract

The fluorescence lifetime technique offers an effective way to resolve fluorescent components with overlapping emission spectra. The presence of multiple fluorescent components in biological compounds can hamper their discrimination. The conventional method based on the nonlinear least-squares technique is unable to consistently determine the correct number of fluorescent components in a fluorescence decay profile. This can limit the applications of the fluorescence lifetime technique in biological assays and diagnoses where more than one fluorescent component is typically encountered. We describe the use of an expectation-maximization (EM) method with joint deconvolution to estimate the fluorescence decay parameters, and the Bayesian information criterion (BIC) to accurately determine the number of fluorescent components. A comprehensive simulation and experimental study is carried out to compare the performance and accuracy of the proposed method. The results show that the EM-BIC method is able to accurately identify the correct number of fluorescent components in samples with weakly fluorescing components.

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¹ 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.² 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.⁵ 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.⁸ 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.⁹ The well-defined statistics in the single-photon counting measurement is also an advantage for data analysis.¹⁰ 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.¹¹ This parametric estimation method, however, has been shown to be incapable of determining the actual number of fluorophores present in a sample.¹² The NLLS method has been found to erroneously overfit the experimental data with a higher lifetime component number.^{2, 13} Global analysis¹⁴ of NLLS-fitted decay profiles is also widely used to improve the discrimination of closely spaced lifetimes.² Despite the improved discrimination, the uncertainties in the extracted lifetimes are still significant.² 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.¹³ 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.⁶ Furthermore, the stretched exponential representation may not be suitable in applications requiring the accurate discrimination and determination of discrete lifetimes for analysis.⁵

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.¹⁵ The emission rate of fluorescence decay is related to the rate of depopulation of the upper-state electrons. The fluorescence dynamics $F (t)$ as a function of time $t$ is conventionally described by the multiexponential equation¹⁶

Eq. 1

F (t) = \sum_{j = 1}^{m} α_{j} \exp (- t ∕ τ_{j}),

where

m

is the number of fluorescence lifetime components, and

α_{j}

and

τ_{j}

are the decay amplitude and fluorescence lifetime of the

j

’th component, respectively. The observed fluorescence decay

R (t)

contains the instrumental effect and is given by the convolution of

F (t)

with the instrumental response function (IRF)

h (t)

written as

Eq. 2

R (t) = F (t) \otimes h (t) .

The IRF contains information about the temporal profile of the laser pulses and the response time characteristics of the measurement system.

In the conventional NLLS method, the estimated parameters $α_{k}$ and $τ_{k}$ are iteratively adjusted until a best fit between the estimated decay $R_{c} (t_{i})$ and the experimental decay $R (t_{i})$ 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 $χ_{r}^{2}$ given by

Eq. 3

χ_{r}^{2} = \frac{1}{a - b} \sum_{i = 1}^{n} \frac{{[R (t_{i}) - R_{c} (t_{i})]}^{2}}{σ_{i}^{2}}

is commonly used to gauge the goodness of fit. Here,

a

is the number of data points,

b

is the number of estimated parameters, and

σ_{i}

is the standard deviation of the photon count in the

i

’th time bin. In TCSPC measurements, the photon statistics follow the Poisson distribution, and

σ_{i}

is approximated by the square root of

R (t_{i})

. The

χ_{r}^{2}

value is expected to approach unity when the estimated function matches closely with the data.

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 (t)$ as

Eq. 4

g (t) = \sum_{j = 1}^{m} \frac{β_{j} \exp (- t ∕ τ_{j})}{τ_{j}} \otimes h (t),

where the sum of all fractional contributions

\sum_{j = 1}^{m} β_{j}

is unity. Here

β_{j}

of the

j

’th component is calculated from the

α

and

τ

vector components using following expression:

Eq. 5

β_{j} = \frac{α_{j} τ_{j}}{\sum_{k = 1}^{m} α_{k} τ_{k}} .

Using the EM algorithm, the $β$ and $τ$ components of a multicomponent fluorescence decay profile with $N$ total photons can be determined based on Eqs. 6, 7, respectively:

Eq. 6

β_{j} = \frac{1}{N} \sum_{i = 1}^{N} P (y_{i j} | t_{i}),

Eq. 7

τ_{j} = \frac{1}{N β_{j}} \sum_{i = 1}^{N} P (y_{i j} | t_{i}) (t_{i} - \frac{M_{i j}}{K_{i j}}),

where

Eq. 8

P (y_{i j} | t_{i}) = \frac{(β_{j} K_{i j} ∕ τ_{j}) \exp (- t_{i} ∕ τ_{j})}{\sum_{r = 1}^{m} (β_{r} K_{i r} ∕ τ_{r}) \exp (- t_{i} ∕ τ_{r})},

Eq. 9

K_{i j} = \int_{0}^{t_{i}} h (t^{'}) \exp (\frac{t^{'}}{τ_{j}}) d t^{'},

Eq. 10

M_{i j} = \int_{0}^{t_{i}} t^{'} h (t^{'}) \exp (\frac{t^{'}}{τ_{j}}) d t^{'} .

The derivations of the equations used in the EM algorithm are given in the Appendix.

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

Eq. 11

λ = \log [L (β, τ; t)] - \frac{m}{2} \log_{e} (N),

where

Eq. 12

L (β, τ; t) = \prod_{i = 1}^{N} g (t_{i} | β, τ) .

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 ns$ , giving a time span of $37.5 ns$ .

Biexponential decay profiles with widely separated lifetimes of $τ_{1} = 1 ns$ and $τ_{2} = 4 ns$ were first generated at varying contribution factors ( $β_{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 $τ_{1} = 2 ns$ and $τ_{2} = 4 ns$ were generated at varying contribution factors ( $β_{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 $τ_{1} = 1 ns$ , $τ_{2} = 4 ns$ , and $β_{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 $χ_{r}^{2}$ 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 $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 to 650 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 mM$ and $50 μ M$ , respectively.

A pulsed laser diode (LDH400, Picoquant) with an excitation wavelength centered at $400 nm$ was used to excite the fluorophore solutions. The laser diode was biased and pulsed at $20 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 - 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 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.²⁴ The IRF of the measurement system was determined by measuring the excitation photon at $400 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 $β$ and $τ$ 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 $χ_{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 $λ$ 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 ( $τ_{1} = 1 ns$ and $τ_{2} = 4 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

The optimum lifetime component number has a χr2 that is closest to unity (shown in bold).

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	$- 6.57859 \times 10^{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

The optimum lifetime component number has the largest λ value (shown in bold).

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 ( $τ_{1} = 1 ns$ and $τ_{2} = 4 ns$ ) and having a sufficiently large photon counts of $1 \times 10^{5}$ . Using these data, the two methods were compared for values of $β_{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.

Fig. 1

Biexponential fluorescence lifetime values (–●–,–◼–) extracted from decay profiles simulated with widely separated lifetimes of $τ_{1} = 1 ns$ and $τ_{2} = 4 ns$ (dashed lines), and a photon count of $1 \times 10^{5}$ using (a) the conventional NLLS method and (b) the EM-BIC method. The model selection accuracy (–▲–) specifies the proportion of 10 decay profiles simulated for each parameter set that has its component number $m$ correctly determined by a model.

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

The model selection accuracy specifies the proportion of 10 decay profiles simulated for each parameter set that has its component number correctly determined by the model.

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 model selection accuracy specifies the proportion of 10 decay profiles simulated for each parameter set that has its component number correctly determined by the model.

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 $τ_{1}$ extracted with the NLLS method for $β_{1} = 0.1$ . Similarly, we can see that the standard deviation of the estimated $τ_{2}$ increases with $β_{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 $β_{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 $β_{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 ( $τ_{1} = 2 ns$ and $τ_{2} = 4 ns$ ), a sufficiently large photon count of $1 \times 10^{5}$ and $β_{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 $β$ value, exhibits an increased standard deviation. The model selection accuracy also drops drastically with the use of the $χ_{r}^{2}$ to resolve the two closely separated fluorescence lifetime components, with a maximum accuracy of only 40% for $β_{1} = 0.1$ [Fig. 2a].

Fig. 2

Biexponential fluorescence lifetime values (–●–,–◼–) extracted from decay profiles simulated with closely spaced lifetimes of $τ_{1} = 2 ns$ and $τ_{2} = 4 ns$ (dashed lines), and a photon count of $1 \times 10^{5}$ using (a) the conventional NLLS method and (b) the EM-BIC method. The model selection accuracy (–▲–) specifies the proportion of 10 decay profiles simulated for each parameter set that has its component number $m$ correctly determined by a model.

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

The model selection accuracy specifies the proportion of 10 decay profiles simulated for each parameter set that has its component number correctly determined by the model.

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

The model selection accuracy specifies the proportion of 10 decay profiles simulated for each parameter set that has its component number correctly determined by the model.

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 $τ_{1} = 1 ns$ and $τ_{2} = 4 ns$ were used to generate the decay profiles for various total photon count. The value of $β_{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 $τ_{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 $τ_{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 $τ_{2}$ ranging from $3.01 to 4.73 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 $τ_{2}$ becomes even larger and the estimate of $τ_{2}$ ranges between 1.95 and $6.13 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].

Fig. 3

Biexponential fluorescence lifetime values (–●–,–◼–) extracted from decay profiles simulated with widely separated lifetimes of $τ_{1} = 1 ns$ and $τ_{2} = 4 ns$ (dashed lines), and $β_{1} = 0.9$ using (a) the conventional NLLS method and (b) the EM-BIC method. The model selection accuracy (–▲–) specifies the proportion of 10 decay profiles simulated for each parameter set that has its component number $m$ correctly determined by a model.

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

The model selection accuracy specifies the proportion of 10 decay profiles simulated for each parameter set that has its component number correctly determined by the model.

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

The model selection accuracy specifies the proportion of 10 decay profiles simulated for each parameter set that has its component number correctly determined by the model.

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 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 ns$ with a $χ_{r}^{2}$ of 1.03, while the estimated lifetime of fluorescein is $4.18 ns$ with a $χ_{r}^{2}$ of 1.09. The extracted lifetime values are in good agreement with the values reported in the literature²⁵ and those from EM-BIC analysis.

Fig. 4

Normalized plots of (a) the system IRF, (b) the fluorescence decay of acridine orange, and (c) the fluorescence decay of fluorescein measured with a total photon count of $2 \times 10^{5}$ .

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 ∕ 2$	$1 ∕ 2$	40
B	$1 ∕ 11$	$10 ∕ 11$	5.5
C	$10 ∕ 11$	$1 ∕ 11$	5.5

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 $m = 2$ are also shown in Fig. 5 .

Fig. 5

Fluorescence lifetime parameters of experimental decay profiles (a) A, (b) B, and (c) C extracted using the conventional NLLS method (◼) and the EM-BIC method (△). The true fluorescence lifetime parameters (×) are also plotted for comparison.

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

The optimum lifetime component number has χr2 that is closest to unity (shown in bold).

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	$- 2.70188 \times 10^{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	$- 3.85917 \times 10^{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	$- 3.52469 \times 10^{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 optimum lifetime component number has the largest λ value (shown in bold).

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 $τ_{1} (= τ_{2} = 2.78 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 $χ_{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 $χ_{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 (t)$ can be described as

Eq. 13

g (t) = \sum_{j = 1}^{m} \frac{β_{j} \exp (- t ∕ τ_{j})}{τ_{j}} \otimes h (t) = \sum_{j = 1}^{m} \frac{β_{j} L_{j} (t, τ_{j}) \exp (- t ∕ τ_{j})}{τ_{j}},

where

Eq. 14

L_{j} (t, τ_{j}) = \int_{0}^{t} h (α) \exp (\frac{α}{τ_{j}}) d α .

Let

y_{i j}

be an indicator variable defined as follows:

y_{i j} = {\begin{cases} 1 & if t_{i} is generated by j th component \\ 0 & otherwise . \end{cases}

Then the likelihood of the complete data can be derived from Eq. 13 as

Eq. 15

p (T | θ) = \prod_{i = 1}^{n} \prod_{j = 1}^{m} {[\frac{β_{j} K_{i j}}{τ_{j}} \exp (\frac{- t_{i}}{τ_{j}})]}^{y_{i j}}

where

T = [t, y_{i j}]

is the complete data,

θ = [β, τ]

is the parameter set, and

n

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 $Q$ function is defined as the conditional mean of log-likelihood function as follows:

Eq. 16

Q (θ_{k} | θ_{k - 1}) = E [p (T | θ_{k}) | t, θ_{k - 1}] = \sum_{i = 1}^{n} \sum_{j = 1}^{m} P (y_{i j} | t_{i}),

where

θ_{k}

and

θ_{k - 1}

are the current and previously estimated parameters, respectively.

Here $P (y_{i j} | t_{i})$ is the probability that $t_{i}$ is generated by the $j$ ’th component and is given by

Eq. 17

P (y_{i j} | t_{i}) = \frac{(β_{j} K_{i j} ∕ τ_{j}) \exp (- t_{i} ∕ τ_{j})}{\sum_{r = 1}^{m} (β_{r} K_{i r} ∕ τ_{r}) \exp (- t_{i} ∕ τ_{r})},

where

Eq. 18

K_{i j} = \int_{0}^{t_{i}} h (t^{'}) \exp (\frac{t^{'}}{τ_{j}}) d t^{'} .

In the maximization step, the EM algorithm updates the estimated parameters

β_{j}

and

τ_{j}

by maximizing the derived

Q

function subjected to the constraint

\sum_{j = 1}^{m} β_{j} = 1

. With the use of Lagrange multipliers, the solutions are

Eq. 19

β_{j} = \frac{1}{n} \sum_{i = 1}^{n} P (y_{i j} | t_{i}),

Eq. 20

τ_{j} = \frac{1}{n β_{j}} \sum_{i = 1}^{n} P (y_{i j} | t_{i}) (t_{i} - \frac{M_{i j}}{K_{i j}}),

where

Eq. 21

M_{i j} = \int_{0}^{t_{i}} t^{'} h (t^{'}) \exp (\frac{t^{'}}{τ_{j}}) d t^{'} .

In the preceding derivation,

K_{i j}

and

M_{i j}

are assumed to be constants, where

τ_{j}

is taken as the old estimate. This approximation is valid since

M_{i j} ∕ K_{i j} < t_{i}

. These expectation and maximization steps are then repeated until convergence is achieved.

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 $(A^{*} STAR)$ and a manpower cofunding (Grant No. M48040149) from the Nanyang Technological University, Singapore.

References

1.

M.-A. Mycek and B. W. Pogue, Handbook of Biomedical Fluorescence, Marcel Dekker, New York (2003). Google Scholar

2.

J. R. Lakowicz, Principles of Fluorescence Spectroscopy, Kluwer Academic/Plenum, New York (1999). Google Scholar

3.

H. Schneckenburger, K. Stock, M. Lyttek, W. S. L. Strauss, and R. Sailer, “Fluorescence lifetime imaging (FLIM) of rhodamine 123 in living cells,” Photochem. Photobiol. Sci., 3 (1), 127 –131 (2004). https://doi.org/10.1039/b306129a 1474-905X Google Scholar

4.

D. Sud, W. Zhong, D. G. Beer, and M. A. Mycek, “Time-resolved optical imaging provides a molecular snapshot of altered metabolic function in living human cancer cell models,” Opt. Express, 14 (10), 4412 –4426 (2006). https://doi.org/10.1364/OE.14.004412 1094-4087 Google Scholar

5.

D. E. Schwartz, P. Gong, and K. L. Shepard, “Time-resolved Forster-resonance-energytransfer DNA assay on an active CMOS microarray,” Biosens. Bioelectron., 24 (3), 383 –390 (2008). https://doi.org/10.1016/j.bios.2008.04.015 0956-5663 Google Scholar

6.

P. J. Tadrous, J. Siegel, P. M. W. French, S. Shousha, E. N. Lalani, and G. W. H. Stamp, “Flourescence lifetime imaging of unstained tissues: early results in human breast cancer,” J. Pathol., 199 (3), 309 –317 (2003). https://doi.org/10.1002/path.1286 0022-3417 Google Scholar

7.

J. A. Jo, Q. Fang, T. Papaioannou, J. D. Baker, A. H. Dorafshar, T. Reil, J. H. Qiao, M. C. Fishbein, J. A. Freischlag, and L. Marcu, “Laguerre-based method for analysis of time-resolved fluorescence data: application to in-vivo characterization and diagnosis of atherosclerotic lesions,” J. Biomed. Opt., 11 (2), 021004 (2006). https://doi.org/10.1117/1.2186045 1083-3668 Google Scholar

8.

L. Kelbauskas and W. Dietel, “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,” 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

9.

E. Gratton, S. Breusegem, J. Sutin, Q. Ruan, and N. Barry, “Fluorescence lifetime imaging for the two-photon microscope: time-domain and frequency-domain methods,” J. Biomed. Opt., 8 (3), 381 –390 (2003). https://doi.org/10.1117/1.1586704 1083-3668 Google Scholar

10.

N. Boens, W. Qin, N. Basaric, J. Hofkens, M. Ameloot, J. Pouget, J. P. Lefévre, B. Valeur, E. Gratton, M. VandeVen, N. D. Silva Jr, Y. Engelborghs, K. Willaert, A. Sillen, G. Rumbles, D. Phillips, A. J. W. G. Visser, A. Van Hoek, J. R. Lakowicz, H. Malak, I. Gryczynski, A. G. Szabo, D. T. Krajcarski, N. Tamai, and A. Miura, “Fluorescence lifetime standards for time and frequency domain fluorescence spectroscopy,” Anal. Chem., 79 (5), 2137 –2149 (2007). https://doi.org/10.1021/ac062160k 0003-2700 Google Scholar

11.

A. Grinvald and I. Z. Steinberg, “On the analysis of fluorescence decay kinetics by the method of least-squares,” Anal. Biochem., 59 (2), 583 –598 (1974). https://doi.org/10.1016/0003-2697(74)90312-1 0003-2697 Google Scholar

12.

D. R. James and W. R. Ware, “A fallacy in the interpretation of fluorescence decay parameters,” Chem. Phys. Lett., 120 (4–5), 455 –459 (1985). https://doi.org/10.1016/0009-2614(85)85640-2 0009-2614 Google Scholar

13.

K. C. Benny Lee, J. Siegel, S. E. D. Webb, S. Lévêque-Fort, M. J. Cole, R. Jones, K. Dowling, M. J. Lever, and P. M. W. French, “Application of the stretched exponential function to fluorescence lifetime imaging,” Biophys. J., 81 (3), 1265 –1274 (2001). https://doi.org/10.1016/S0006-3495(01)75784-0 0006-3495 Google Scholar

14.

J. M. Beechem, M. Ameloot, and L. Brand, “Global analysis of fluorescence decay surfaces: excited-state reactions,” Chem. Phys. Lett., 120 (4–5), 466 –472 (1985). https://doi.org/10.1016/0009-2614(85)85642-6 0009-2614 Google Scholar

15.

B. Valeur, Wiley InterScience (Online service), Molecular Fluorescence Principles and Applications, Wiley-VCH, New York (2001). Google Scholar

16.

M. G. Badea and L. Brand, “Time-resolved fluorescence measurements,” J. Membr. Biol., 61 378 –425 (1979). 0022-2631 Google Scholar

17.

J. Enderlein and R. Erdmann, “Fast fitting of multi-exponential decay curves,” Opt. Commun., 134 (1–6), 371 –378 (1997). https://doi.org/10.1016/S0030-4018(96)00384-7 0030-4018 Google Scholar

18.

A. P. Dempster, N. M. Laird, and D. B. Rubin, “Maximum likelihood from incomplete data via the EM algorithm,” J. R. Stat. Soc. Ser. B (Methodol.), 39 (1), 1 –38 (1977). 0035-9246 Google Scholar

19.

J. Salojarvi, K. Puolamaki, and S. Kaski, “Expectation maximization algorithms for conditional likelihoods,” 753 –760 (2005). Google Scholar

20.

Z. Liu, J. Almhana, V. Choulakian, and R. McGorman, “Recursive EM algorithm for finite mixture models with application to Internet traffic modeling,” 198 –207 (2004). Google Scholar

21.

M. W. Graham and D. J. Miller, “Unsupervised learning of parsimonious mixtures on large spaces with integrated feature and component selection,” IEEE Trans. Signal Process., 54 (4), 1289 –1303 (2006). https://doi.org/10.1109/TSP.2006.870586 1053-587X Google Scholar

22.

L. P. Watkins and H. Yang, “Detection of intensity change points in time-resolved single-molecule measurements,” J. Phys. Chem. B, 109 (1), 617 –628 (2005). https://doi.org/10.1021/jp0467548 1089-5647 Google Scholar

23.

J. K. Ghosh and T. Samanta, “Model selection—an overview,” Curr. Sci., 80 (9), 1135 –1144 (2001). 0011-3891 Google Scholar

24.

W. Becker, Advanced Time-Correlated Single Photon Counting Techniques, Springer, Berlin, New York (2005). Google Scholar

25.

A. G. Ryder, S. Power, T. J. Glynn, and J. J. Morrison, “Time-domain measurement of fluorescence lifetime variation with pH,” Proc. SPIE, 4259 102 –109 (2001). https://doi.org/10.1117/12.432487 0277-786X Google Scholar

26.

R. Cubeddu, D. Comelli, C. D’Andrea, P. Taroni, and G. Valentini, “Time-resolved fluorescence imaging in biology and medicine,” J. Phys. D: Appl. Phys., 35 (9), R61 –R76 (2002). https://doi.org/10.1088/0022-3727/35/9/201 0022-3727 Google Scholar

27.

D. Elson, J. Requejo-Isidro, I. Munro, F. Reavell, J. Siegel, K. Suhling, P. Tadrous, R. Benninger, P. Lanigan, J. McGinty, C. Talbot, B. Treanor, S. Webb, A. Sandison, A. Wallace, D. Davis, J. Lever, M. Neil, D. Philips, G. Stamp, and P. French, “Time domain fluorescence lifetime imaging applied to biological tissue,” Photochem. Photobiol. Sci., 3 (8), 795 –801 (2004). https://doi.org/10.1039/b316456j 1474-905X Google Scholar

28.

M. L. Johnson, “Why, when, and how biochemists should use least squares,” Anal. Biochem., 206 (2), 215 –225 (1992). https://doi.org/10.1016/0003-2697(92)90356-C 0003-2697 Google Scholar

Citation Download Citation

C. Y. Fu, Beng-Koon Ng, and Sirajudeen Gulam Razul "Fluorescence lifetime discrimination using expectation-maximization algorithm with joint deconvolution," Journal of Biomedical Optics 14(6), 064009 (1 November 2009). https://doi.org/10.1117/1.3258835

Published: 1 November 2009

Access the abstract

JOURNAL ARTICLE
10 PAGES

DOWNLOAD PAPER SAVE TO MY LIBRARY

GET CITATION

CITATIONS

Cited by 15 scholarly publications and 1 patent.

Explore citations on Lens.org

RIGHTS & PERMISSIONS

Get copyright permission Get copyright permission on Copyright Marketplace

KEYWORDS

Luminescence

Photon counting

Expectation maximization algorithms

Computer simulations

Deconvolution

Monte Carlo methods

Biological research

1.

Introduction

2.

Theory

2.1.

NLLS Method

Eq. 1

Eq. 2

Eq. 3

2.2.

Proposed EM Algorithm with Joint Deconvolution

Eq. 4

Eq. 5

Eq. 6

Eq. 7

Eq. 8

Eq. 9

Eq. 10

Eq. 11

Eq. 12

3.

Materials and Methods

3.1.

Simulation

3.2.

Experimental Technique

4.

Results

4.1.

Simulation

Table 1

Table 2

Fig. 1

Table 3

Table 4

Fig. 2

Table 5

Table 6

Fig. 3

Table 7

Table 8

4.2.

Experiment

Fig. 4

Table 9

Fig. 5

Table 10

Table 11

5.

Discussion

6.

Conclusion

Appendices

Appendix

Eq. 13

Eq. 14

Eq. 15

Eq. 16

Eq. 17

Eq. 18

Eq. 19

Eq. 20

Eq. 21

Acknowledgments

References

Show All Keywords

Keywords/Phrases

Search In:

Publication Years