Estimation of multiexponential fluorescence decay parameters using compressive sensing

Sejung Yang; Joohyun Lee; Youmin Lee; Minyung Lee; Byung-Uk Lee

doi:10.1117/1.JBO.20.9.096003

3 September 2015 Estimation of multiexponential fluorescence decay parameters using compressive sensing

Sejung Yang, Joohyun Lee, Youmin Lee, Minyung Lee, Byung-Uk Lee

Author Affiliations +

Journal of Biomedical Optics, Vol. 20, Issue 9, 096003 (September 2015). https://doi.org/10.1117/1.JBO.20.9.096003

Abstract

Fluorescence lifetime imaging microscopy (FLIM) is a microscopic imaging technique to present an image of fluorophore lifetimes. It circumvents the problems of typical imaging methods such as intensity attenuation from depth since a lifetime is independent of the excitation intensity or fluorophore concentration. The lifetime is estimated from the time sequence of photon counts observed with signal-dependent noise, which has a Poisson distribution. Conventional methods usually estimate single or biexponential decay parameters. However, a lifetime component has a distribution or width, because the lifetime depends on macromolecular conformation or inhomogeneity. We present a novel algorithm based on a sparse representation which can estimate the distribution of lifetime. We verify the enhanced performance through simulations and experiments.

1. Introduction

Fluorescence lifetime imaging microscopy (FLIM) is an imaging technique observing the decreasing rate of emitted photon counts with the lapse of time, after applying a laser beam to a cell which is labeled with fluorophores. It provides essential information for cell analysis and can be used for intracellular measurements of fluorescence resonance energy transfer (FRET).¹ The FLIM acquisition methods can be categorized into time-domain and frequency-domain methods. Time-domain methods, such as time correlated single-photon counting (TCSPC) techniques, are more robust to scattering and can measure short lifetimes better than frequency-domain methods.²^,³ FLIM can be obtained using analog mean delay, which is intensity weighted mean delay time; however, it cannot estimate multiple lifetimes; it measures the intensity weighted average lifetime.⁴ Recently, many researchers have used FLIM with multiphoton confocal microscopy to noninvasively monitor changes in metabolism. Skala et al. have observed cellular nicotinamide adenine dinucleotide and flavin adenine dinucleotide lifetimes at the earliest stages of cancer development.⁵ Roberts et al. have investigated coexisting drugs and their metabolites in the skin using FLIM techniques.⁶ Deka et al. have also used FLIM to measure the metabolic rate of cells for in vivo wound healing diagnosis.⁷

Accurate estimation of fluorescence lifetime is crucial to the FLIM imaging. A simple estimation method is a nonlinear least-squares (NLLS) fitting. The NLLS algorithm is a conventional method to estimate lifetimes by minimizing the weighted sum of the squared residuals.⁸^–¹¹ However, this method is not effective in low counts, since it does not adopt a Poisson noise probability density function of the photon counting process. FRET has donor and acceptor fluorophores; however, NLLS does not accurately estimate multiexponential decays. Pearson and Neyman’s $χ^{2}$ -based estimators have been proposed to overcome the Gaussian noise assumption.¹² Sasaki and Masuhara proposed a convolved autoregressive model for fitting multiexponential decay curves.¹³ Apanasovich and Novikov also proposed a fitting method, which is mathematically equivalent to Sasaki’s and Masuhara method.¹⁴ An improved version was presented by Enderlein and Erdmann, which estimates the optimal fit parameters for multiexponential decay curves by examining an error function.¹⁵

Recently, the maximum-likelihood estimator (MLE) based on the Poisson-likelihood function has proven to be superior in terms of unbiasedness and efficiency from theoretical analysis and computer simulations.¹⁶ Although MLE is an excellent estimator to fit the multiexponential decay parameters, it is sensitive to the initial parameter values. Thus, the choice of the numerical optimization algorithm is significantly important.

The number of fluorescent lifetime components is quite sparse; therefore, a compressed sensing algorithm is suitable for the estimation. We observed that a sparse Poisson intensity reconstruction algorithm (SPIRAL),¹⁷ which is a penalized maximum Poisson-likelihood estimator, is effective for FLIM. SPIRAL has been proposed for the estimation of image intensity from limited photon counts such as tomographic reconstruction in nuclear medicine. Using this algorithm, a unique solution to the minimization can be obtained.

This paper proposes a novel algorithm for multiexponential fluorescence decay parameters that can estimate distribution of the lifetime in FLIM. To the best of our knowledge, compressed sensing has not been applied to FLIM yet. We evaluate the performance of our algorithm from simulations and TCSPC data from confocal microscopy, which suffer from very low photon counts. Experimental results demonstrate the efficacy of the proposed approach.

2. Prior Work

Before proposing a new method for estimating lifetime based on a Poisson distribution with a regularizing function, a brief review of prior work for conventional lifetime estimation is provided. We especially discuss two approaches: modified Neyman’s nonlinear least squares (MN-NLLS) and fast fitting, which are embedded in the widely used TCSPC software and SPCImage.¹⁸

2.1.

Modified Neyman’s Nonlinear Least-Squares Method

NLLS algorithm is used in over-determined systems and this algorithm gives the solution that minimizes the sum of squared error, which is the difference between the observed data $y = {[y_{1}, y_{2}, \dots, y_{N}]}^{T}$ and the fitting model $g = {[g (t_{1}), g (t_{2}), \dots, g (t_{N})]}^{T}$ .⁸^–¹⁰ The objective function is as follows:

Eq. (1)

Φ_{NLLS} = \sum_{i = 1}^{N} {[g (t_{i}) - y_{i}]}^{2} .

The NLLS algorithm is a widely used method to obtain the distribution of fluorescence lifetime, since it is simple and fast. However, NLLS assumes that the measured data are acquired under additive Gaussian noise with equal variance. It is not suitable for photon counts which follow a Poisson distribution. Thus, it is less accurate for fluorescence lifetime estimation, as demonstrated by Kim and Seok.¹⁶

Neyman’s $χ^{2}$ function has been introduced using the property of the Poisson process; therefore, the variance is proportional to the intensity. However, the original Neyman’s $χ^{2}$ function is not defined when $y_{i}$ is zero.¹⁶ As a remedy, a modified Neyman’s NLLS method is proposed, which is obtained by substituting the variance in the denominator with $\max (y_{i}, 1)$ . This objective function can be written as

Eq. (2)

Φ_{MN - NLLS} = \sum_{i = 1}^{N} \frac{{[g (t_{i}) - y_{i}]}^{2}}{\max (y_{i}, 1)} .

Observation of multiexponential decay is modeled as

Eq. (3)

g (t_{i}) = [\sum_{j = 1}^{M} f_{j} \exp (- t_{i} / τ_{j})] * h_{i},

where

f = {[f_{1}, f_{2}, \dots, f_{M}]}^{T}

is the amplitude of each lifetime component corresponding to

τ = [τ_{1}, τ_{2}, \dots, τ_{M}]

, and

h_{i}

stands for the instrument response function

h (t_{i})

. The observation model vector

g

can be represented as

A f

, where the matrix

A

consists of the convolution of the exponential decay functions with each lifetime component

τ = [τ_{1}, τ_{2}, \dots, τ_{M}]

at measurement time

t = [t_{1}, t_{2}, \dots, t_{N}]

,

Eq. (4)

A = {[\begin{matrix} h_{i} * [\exp (- t_{1} / τ_{1}) & \dots & \exp (- t_{N} / τ_{1})] \\ ⋮ & ⋱ & ⋮ \\ ⋮ & ⋱ & ⋮ \\ h_{i} * [\exp (- t_{1} / τ_{M}) & \dots & \exp (- t_{N} / τ_{M})] \end{matrix}]}^{T} .

2.2.

Fast Fitting Using a Convolved Autoregressive Model

For estimating the decay parameters of photon emission curves convolved with an instrumental response function, a convolved autoregressive model was proposed.¹³^–¹⁵ The discrete number of photons observed in a finite number of time channels $[t_{j} - Δ t / 2, t_{j} + Δ t / 2]$ is given by

Eq. (5)

y_{k} = \sum_{i = 0}^{\infty} \sum_{j = 1}^{M} f_{j} \exp (- i Δ t / τ_{j}) h_{k - i},

where

τ_{j}

and

f_{j}

are the decay constant and the amplitude of the

j

’th component decay.

y_{k}

is the number of observed photon counts, and

h_{k}

is the instrumental impulse response function for the

k

’th time channel. The

z

-transform of Eq. (5) yields

Eq. (6)

\tilde{y} (z) = \sum_{j = 1}^{M} f_{j} \tilde{h} (z) {(1 - z_{j} / z)}^{- 1},

where

z_{j} = \exp (- Δ t / τ_{j})

and

\tilde{h} (z)

is the

z

transform of the estimated instrumental impulse response function

h_{k}

.

Equation (6) can be rearranged as follows:

Eq. (7)

\tilde{y} (z) \prod_{j = 1}^{M} (1 - z_{j} / z) = \sum_{j = 1}^{M} f_{j} \tilde{h} (z) \prod_{i \neq j}^{M} (1 - z_{i} / z) .

The products can be expanded into sums by defining new coefficients $b_{j}$ and $c_{j}$ :

Eq. (8)

\prod_{j = 1}^{M} (1 - z_{j} / z) = 1 - \sum_{j = 1}^{M} \frac{b_{j}}{z^{j}}, \sum_{j = 1}^{M} f_{j} \prod_{i \neq j}^{M} (1 - z_{i} / z) = \sum_{j = 0}^{M - 1} \frac{c_{j}}{z^{j}} .

After rearranging and expanding, Eq. (6) can be rewritten in the following form:

Eq. (9)

\tilde{y} (z) = \sum_{j = 1}^{M} \tilde{y} (z) \frac{b_{j}}{z^{j}} + \sum_{j = 0}^{M - 1} \tilde{h} (z) \frac{c_{j}}{z^{j}} .

The inverse $z$ -transform gives the desired autoregressive convolution:

Eq. (10)

y_{k} = \sum_{j = 1}^{M} b_{j} y_{k - j} + \sum_{j = 0}^{M - 1} c_{j} h_{k - j} .

It can be written in a compact matrix form.

Eq. (11)

y = D p, with y = {[y_{M + 1}, y_{M + 2}, \dots, y_{N}]}^{T}, p = {[b_{1}, b_{2}, \dots, b_{M}, c_{0}, c_{1}, \dots, c_{M - 1}]}^{T}, D = [\begin{matrix} y_{M} & y_{M - 1} & \dots & y_{1} \\ y_{M + 1} & y_{M} & \dots & y_{2} \\ ⋮ & ⋮ & \dots & ⋮ \\ y_{N - 1} & y_{N - 2} & \dots & y_{N - M} \end{matrix} \begin{matrix} h_{M + 1} & h_{M} & \dots & h_{2} \\ h_{M + 2} & h_{M + 1} & \dots & h_{3} \\ ⋮ & ⋮ & \dots & ⋮ \\ h_{n} & h_{N - 1} & \dots & h_{N - M + 1} \end{matrix}],

where

N

is the number of time sampled data points and

M

is the number of lifetime components.

The general form of Eq. (11) is as follows:

Eq. (12)

y = D p + ϵ,

where

ϵ

is a random vector of residuals. The optimal values for the vectors

p

and

y

can be estimated.

Eq. (13)

\hat{y} = D \hat{p},

where

\hat{y}

and

\hat{p}

are the estimated versions of

y

and

p

, respectively. From this notation, we have

ϵ = y - \hat{y}

. The covariance matrix of the vector

ϵ

is

Eq. (14)

cov (ϵ) = ⟨ (ϵ - ⟨ ϵ ⟩) {(ϵ - ⟨ ϵ ⟩)}^{T} ⟩ = V,

where

⟨ ϵ ⟩

denotes the expectation of

ϵ

. The MLE solution

p

is given by

Eq. (15)

p = {(D^{T} V^{- 1} D)}^{- 1} D^{T} V^{- 1} y,

which enables the calculation of the distribution of multiexponential components

f_{j}

from Eq. (8).

3. Proposed Method

3.1.

Problem Formulation

Data collected by counting photons hitting a detector are inherently noisy due to low count levels. Harmany et al. developed SPIRAL for photon limited imaging in the spatial domain such as emission tomography.¹⁷ The algorithm is reported to outperform the state-of-the-art emission tomography algorithms. We found that SPIRAL can be adopted for multiexponential lifetime estimation in the time domain. Under a Poisson process, the observation is modeled as

Eq. (16)

y \sim Poisson ({A f}^{*}),

where

f^{*} = [f_{1}^{*}, f_{2}^{*}, \dots, f_{M}^{*}]

is the amplitude of fluorescence lifetimes to be estimated, the matrix

A

is given in Eq. (3), and

y

is a length

N

vector of the observed photon counts;

f^{*}

and

A

represented the image intensity and a projection matrix for measurements, respectively, in the original SPIRAL.¹⁷

3.2.

Poisson Log-Likelihood

The sensing matrix $A$ represents multiexponential decay in FLIM: each column represents an exponential decay of a single lifetime component, while the row direction is the time axis. Thus, the observed photon count is modeled as a linear superposition of sparse lifetime components $f^{*}$ . In addition, we need to include the effect of impulse response function of the microscopy. Therefore, the sensing matrix $A$ is the same as Eq. (3).

In the Poisson model, the probability of the observation vector is given by

Eq. (17)

p (y | A f) = Π_{i = 1}^{N} \frac{{(e_{i}^{T} {A f}^{*})}^{y_{i}}}{y_{i}!} \exp (- e_{i}^{T} {A f}^{*}),

where

e_{i}

is the

i

’th canonical basis unit vector. Therefore, the negative Poisson log-likelihood is given by

Eq. (18)

- \log [p (y | A f)] = - \log [Π_{i = 1}^{N} \frac{{(e_{i}^{T} A f)}^{y_{i}}}{y_{i}!} \exp (- e_{i}^{T} A f)] = 1^{T} A f - \sum_{i = 1}^{N} \log \frac{{(e_{i}^{T} A f)}^{y_{i}}}{y_{i}!},

where 1 is an

n

-vector of ones. Since

y_{i}

is the constant with respect to

f

, we can ignore the

\log (y_{i}!)

term during optimization. A small parameter

β = 10^{- 10}

is added to avoid the singularity problem of a log function when

f

is zero. The maximum-likelihood estimation for the measurement modeled in Eq. (16) is determined by the minimizing the following equation:

Eq. (19)

F (f) = 1^{T} A f - \sum_{i = 1}^{N} y_{i} \log (e_{i}^{T} A f + β) .

3.3.

Objective Function

We use an objective function which is a sum of the negative Poisson log-likelihood function and $L 1$ -norm regularization term. By adding the $L 1$ -norm penalty term to the objective function, we can obtain a sparse solution $f^{*}$ .¹⁷^,¹⁹^–²³ Exponential decay of the fluorescence can be modeled as a sum of a finite number of lifetime components; therefore, sparse representation is an excellent choice.

Function $F (f)$ is the negative Poisson log-likelihood function of the observed photon counts $y$ from the Poisson process. We adopt a penalized negative Poisson log-likelihood objective function with a non-negativity constraint:

Eq. (20)

Φ (f) = F (f) + ω {‖ f ‖}_{1} f^{*} = argmin Φ (f), subject to f \geq 0,

where

ω

is the weighting factor for the regularization term. A sequential quadratic approximation to

F (f)

is employed for SPIRAL, and the global convergence and uniqueness of the solution are established under a mild set of assumptions.¹⁷

The SPIRAL solves the penalized Poisson log-likelihood function; therefore, it can find multiple exponential decay time constants and amplitudes that are sparse and non-negative.

4. Numerical Simulations

To validate the proposed algorithm, a simulation study is conducted in order to estimate multifluorescence lifetimes and to compare these results to existing methods. A biexponential decay model for a simulation is as follows:

Eq. (21)

y_{i} = [a_{1} \exp (- \frac{t_{i}}{τ_{1}}) + a_{2} \exp (- \frac{t_{i}}{τ_{2}})] * h_{i},

where

a

and

τ

denote the photon count and lifetime, respectively.

The simulation uses two fluorescence components with lifetimes of 2 and 6 ns. To adjust the Poisson noise variance, we employed a scaling factor $α$ :

Eq. (22)

x = Poisson {α y} .

The variance and mean are proportional to $α$ ; therefore, signal power is proportional to $α^{2}$ , while the signal-to-noise power ratio is proportional to $α$ . As $α$ decreases, the signal-to-noise ratio decreases with $α$ . The observed time interval is from 0 to 19.1 ns and the sampling interval is 0.1 ns. The number of lifetime components for sensing matrix $A$ is 196, which covers from 32 ps to 32 ns in logscale. Since the conventional MLE is sensitive to initial values, we set the initial parameter values to the true values, while the initial values are set to zero for other methods. Conventional methods need to fix the number of lifetime components; therefore, they find the best fit after increasing the lifetime component from one to three in most cases. On the other hand, we obtain the distribution of lifetime components using the proposed algorithm. From the distribution, the lifetime is estimated as the amplitude weighted mean of the distribution in each contiguous group and the corresponding intensity is the sum of the intensity components of the group.

Table 1 shows the simulation results of lifetime estimation using the MN-NLLS method,⁸^–¹⁰ fast fitting method,¹⁵ MLE method initialized with ground truth,¹⁶ and our algorithm. Mean and standard deviation are as shown in Table 1 as calculated from 500 realizations of the Poisson random variable. The accuracy of algorithms decreases as the scaling Poisson variance factor $α$ decreases. As shown in Table 1, the MN-NLLS and fast fitting methods are quite sensitive to noise: the estimated mean value significantly deviates from true lifetime, especially when $α = 0.2$ . In comparison, the proposed method is robust to noise: the relative error for the 2 ns component is 33%, 37.5%, 7%, and 7.5% for NM-NLLS, fast fitting, MLE, and the proposed method. The proposed method provides similar results the as MLE method initialized with the ground truth, since the proposed algorithm is based on MLE. However, the MLE method is not practical since the initial solution must be the ground truth. Otherwise, it falls in a local minimum, which is quite different from the desired solution.

Table 1

Simulation results of fluorescence lifetime estimation using modified Neyman’s nonlinear least squares (MN-NLLS) algorithm,12 fast fitting,15 maximum-likelihood method,16 and our algorithm (true values of lifetimes: 2, 6 ns and true values of amplitudes: 100, 100).

α	True lifetimes and amplitudes	MN-NLLS12: mean (standard deviation)	Fast fitting15: mean (standard deviation)	MLE with ground truth16: mean (standard deviation)	The proposed algorithm: mean (standard deviation)
no noise	$τ_{1} = 2 ns$	2.00	1.99	2.00	2.00
	$τ_{2} = 6 ns$	6.00	5.91	6.00	6.00
	$a_{1} = 100$	100.00	101.58	100.00	100.00
	$a_{2} = 100$	100.00	98.42	100.00	100.00
1	$τ_{1} = 2 ns$	1.77 (0.33)	2.30 (0.61)	2.02 (0.24)	2.02 (0.26)
	$τ_{2} = 6 ns$	5.33 (0.21)	6.32 (1.15)	6.03 (0.18)	6.02 (0.20)
	$α a_{1} = 100$	82.02 (12.03)	111.05 (28.53)	101.14 (8.50)	100.28 (9.27)
	$α a_{2} = 100$	117.65 (12.73)	88.93 (27.93)	909.05 (8.71)	99.08 (9.58)
0.5	$τ_{1} = 2 ns$	2.33 (0.81)	2.36 (0.77)	2.02 (0.34)	2.00 (0.34)
	$τ_{2} = 6 ns$	6.10 (1.33)	6.46 (1.63)	6.03 (0.26)	6.03 (0.26)
	$α a_{1} = 50$	55.54 (16.82)	57.12 (19.08)	50.81 (5.95)	50.85 (5.49)
	$α a_{2} = 50$	42.28 (19.95)	45.31 (19.58)	49.25 (6.24)	49.52 (6.12)
0.2	$τ_{1} = 2 ns$	2.66 (0.69)	2.75 (1.02)	1.86 (0.75)	1.85 (0.67)
	$τ_{2} = 6 ns$	7.71 (1.75)	6.68 (2.64)	5.98 (0.51)	6.01 (0.48)
	$α a_{1} = 20$	29.73 (5.77)	20.51 (7.67)	20.21 (4.68)	20.92 (3.72)
	$α a_{2} = 20$	9.14 (6.90)	15.76 (9.43)	20.60 (4.91)	20.24 (4.53)

Typical execution times of various estimation methods in Table 1 are approximately as follows: MN-NLLS method 0.41 s, fast fitting method 0.72 s, maximum-likelihood method 2.10 s, and the proposed method 1.25 s. The compared algorithms were implemented by MATLAB®. The execution times were obtained with an Intel Core i7-3770 processor running at 3.40 GHz.

Lifetime depends on the macromolecular conformation and also heterogeneity; therefore, the lifetime component has a distribution or width depending on the conformation. Since the proposed method can estimate the distribution of the lifetime, it results in a more accurate estimation for cell where the number of lifetime components is not known.

5. Experimental Results

We compare the estimated lifetime of MN-NLLS and fast fitting algorithms to that of the proposed algorithm with FLIM data. A Nikon TE2000-S microscope and Picoquant LDH-P-C-440M laser were used for FLIM imaging. The laser operates at a wavelength of 442 nm and a repetition rate of 20 MHz. A lens for the microscope has a numerical aperture 1.4, $60 \times$ . The FLIM image scanner is atomic force microscopy controller by PSIA and the fluorometric detector is the R3809U-07 model by Hamamatsu Photonics. A SPC-830 board by Becker-Hickl is used for data processing. A 473-nm long path cut-off filter (Semrock) was used in front of the detector to block the excitation beam. Photon counts were continuously monitored by a photon counter (EG&G Ortec). The thin film containing dihexyloxacarbocyanine iodide (DOCI) and Coumarin 343 dyes were prepared on a cover glass by the spin-coating method. FITC- $A β 11$ -25 and FITC- $A β 11$ -25-DABMI to the HeLa cells were incubated for 20 h.

First, we compare the resolving capabilities with a double exponential decay: we coated a mixture of DOCI and Coumarin 343, whose measured lifetimes are 1.59 and 2.82 ns respectively, on polyacrylic acid polymer. As a measure of quality, we use Neyman’s $χ^{2}$ , defined as

Eq. (23)

χ^{2} = \frac{1}{N - p} \sum_{i = 1}^{N} \frac{{(y_{i} - g_{i})}^{2}}{y_{i}},

for

y_{i}

larger than a threshold. In our experiments, the threshold is set to 5. The numbers of the fit parameters and data points are

p

and

N

. The lifetime estimation is accurate when

χ^{2}

is close to one.

To experiment with the biexponential decay, we mixed DOCI and Coumarin 343, which have different lifetimes. TCSPC collected data were analyzed with MN-NLLS, the fast fitting method, and the proposed algorithm, and results are summarized in Table 2. We selected six image positions for brevity: the rest of the image is quite similar to the shown data. The proposed method can more accurately estimate double exponential fluorescence lifetimes than others. The MN-NLLS method estimated double exponential components at two positions among the six samples and the fast fitting method was able to resolve the double exponential fluorescence lifetimes at only one position. The performance of our algorithm is satisfactory in terms of $χ^{2}$ , which is close to one. The size of the acquired FLIM image is $248 \times 254 pixels$ and the total time sampling is 256 points with time interval 85.9 ps. The proposed method could resolve the biexponential decay for 81% of them, whereas the MN-NLLS and fast fitting methods detected the biexponential for 37% and 17%, respectively. Even the proposed SPIRAL-based algorithm could not detect the biexponential decay for the all data points due to low photon counts, but it is still much better than the conventional algorithms.

Table 2

Comparison of estimated fluorescence lifetimes of DOCI and Coumarin 343 fluorescence lifetime imaging microscopy (FLIM) data from MN-NLLS, fast fitting, and the proposed algorithm at arbitrarily selected positions.

Pixel position, initial photon count	MN-NLLS12: lifetime/amplitude (χ2)	Fast fitting15: lifetime/amplitude (χ2)	The proposed method: lifetime/amplitude (χ2)
(100, 50) Photon count: 70	$1.54 ns / 65.07$ , $4.62 ns / 2.31$ (1.34)	$1.77 ns / 85.90$ (0.71)	$1.85 ns / 67.65$ , $5.06 ns / 1.10$ (0.96)
(100, 100) Photon count: 74	$1.00 ns / 58.06$ , $2.34 ns / 27.39$ (1.43)	$1.76 ns / 120.88$ (0.83)	$1.67 ns / 66.26$ , $3.31 ns / 6.52$ (1.14)
(100, 200) Photon count: 107	$1.85 ns / 98.92$ (0.73)	$1.94 ns / 190.31$ (0.62)	$1.77 ns / 101.23$ , $3.23 ns / 9.79$ (0.94)
(200, 50) Photon count: 73	$1.79 ns / 58.26$ (1.23)	$1.98 ns / 74.97$ (1.06)	$1.85 ns / 67.26$ (1.05)
(200, 100) Photon count: 78	$1.67 ns / 74.04$ (1.20)	$2.07 ns / 115.04$ (0.88)	$1.78 ns / 54.47$ , $2.78 ns / 16.39$ (1.13)
(200, 200) Photon count: 122	$1.76 ns / 99.65$ (1.29)	$1.44 ns / 82.87$ , $2.74 ns / 47.60$ (1.05)	$1.41 ns / 93.32$ , $2.91 ns / 30.55$ (1.08)

Figure 1 shows the lifetime distribution at the (100, 100) pixel position obtained from the proposed method. We can observe two distinct peaks in the distribution. Each estimated lifetime is the amplitude weighted mean of the components and the intensity is the sum of the constituents. Figure 2 presents the observed photon counts at the (100, 50) pixel position of the model acquired from TCSPC and the exponential decay model from estimated lifetime components over time. The black solid line represents the observed photon counts. The red line with a circle, the blue line with a plus, and the green line with a diamond come from the calculation of the biexponential decay estimated from the proposed method, MN-NLLS, and fast fitting method, respectively. The proposed method accurately estimates double exponents; thus, the fitted exponential decay correctly matches the observed data. While the curves from fast fitting and MN-NLLS closely follow the proposed method for the initial stage for the time interval between 0 and 4 ns, these results deviate from the measured photon counts and the proposed method after roughly 4 ns. The fast fitting results diminish faster than those of the MN-NLLS method.

Fig. 1

Lifetime distribution at (100, 100) pixel position from the proposed algorithm. Two distinct lifetime components, 1.67 and 3.31 ns, are resolved.

Fig. 2

Measured and estimated results at (100, 50) pixel position in logarithmic scale; black solid line: measured photon counts, red-dotted line: the estimated result from our proposed method, blue-dashed line: the result from modified Neyman’s nonlinear least squares (MN-NLLS), green dash-dotted line: the result from fast fitting.

Now we demonstrate the performance of the SPIRAL-based algorithm on FLIM-based FRET, which is a powerful tool that provides information about the vicinity between donor- and acceptor-labeled molecules.²⁴ FRET occurs when donor and acceptor molecules are placed within a distance that is close ( $< 10 nm$ ) enough to transfer the excited state energy from an initially excited donor to an acceptor.²⁵ After the energy transfer between a donor and an acceptor, the lifetime of the fluorophore is shortened. We adopted a fluorescence amyloid ( $A β$ ) probe, in which fluoresceine isothiocyanate (FITC) is a donor at the end of the $N$ -terminus of $A β 11$ -25 and 4-dimethylaminophenylazophenyl-4’-maleimide (DABMI) is an acceptor at the end of the C-terminus. FRET provides us not only the conformational information of $A β$ , but also the localization of the peptide in the cellular environment. We observed HeLa cells that were grown in Dulbecco’s modified Eagles medium supplemented with 10% (v/v) fetal bovine serum and 1% penicillin. For the confocal and fluorescence lifetime imaging, the cells were plated in a 35-mm glass bottom dish (MatTek Corporation) at low density and grown to 80% confluency on culture dishes at 37°C in a humidified 5% ${CO}_{2}$ atmosphere incubator. At first, we analyze a FLIM image without the DABMI acceptor to observe the data without FRET. The fluorescence of the FITC is known to follow a double exponential decay.²⁶^,²⁷ As shown in Table 3, the proposed algorithm estimates double exponential lifetimes better than MN-NLLS and fast fitting algorithms. Table 3 shows the estimation results at six typical pixel data: MN-NLLS and fast fitting methods cannot resolve biexponential lifetimes for half of the pixels, while the presented method estimates biexponential decay at five positions. The size of the acquired FLIM image is $256 \times 256 pixels$ . The proposed method could resolve the biexponential decay for 97% of them, whereas the MN-NLLS and fast fitting methods detected biexponential for 56% and 68%, respectively. Please note that our method does not estimate just two lifetime components; it estimates the distribution of a lifetime, which is a more accurate model for the fluorescence decay.

Table 3

Comparison of estimated fluorescence lifetimes of fluoresceine isothiocyanate (FITC) FLIM data for HeLa cells from MN-NLLS, fast fitting, and the proposed algorithm at arbitrary positions.

Pixel position, initial photon count	MN-NLLS12: lifetime/amplitude (χ2)	Fast fitting15: lifetime/amplitude (χ2)	The proposed method: lifetime/amplitude (χ2)
(140, 120) Photon count: 14	$3.13 ns / 8.60$ (2.18)	$0.90 ns / 5.09$ , $3.85 ns / 31.19$ (1.28)	$1.42 ns / 5.33$ , $3.56 ns / 11.15$ (1.12)
(150, 120) Photon count: 21	$3.00 ns / 18.75$ (1.43)	$2.67 ns / 67.42$ (0.83)	$1.47 ns / 2.07$ , $3.10 ns / 20.90$ (0.86)
(160, 120) Photon count: 23	$2.97 ns / 23.87$ (1.02)	$0.93 ns / 15.74$ , $3.53 ns / 56.68$ (0.82)	$2.96 ns / 26.52$ (1.02)
(170, 120) Photon count: 22	$1.88 ns / 19.01$ , $8.57 ns / 1.57$ (1.57)	$2.82 ns / 52.06$ (1.03)	$2.33 ns / 16.04$ , $6.25 ns / 3.67$ (1.05)
(180, 120) Photon count: 18	$3.45 ns / 9.14$ (2.21)	$2.97 ns / 41.99$ (1.08)	$1.47 ns / 4.33$ , $3.2 ns / 13.82$ (1.07)
(140, 130) Photon count: 11	$2.22 ns / 8.68$ , $5.46 ns / 2.44$ (2.47)	$3.48 ns / 43.36$ (1.14)	$3.31 ns / 10.63$ , $5.43 ns / 1.83$ (1.05)
(150, 130) Photon count: 22	$1.04 ns / 19.34$ , $4.14 ns / 9.30$ (1.25)	$2.40 ns / 49.18$ , $7.98 ns / 13.83$ (0.69)	$1.64 ns / 12.32$ , $4.40 ns / 10.18$ (0.96)
(160, 130) Photon count: 25	$1.44 ns / 26.77$ , $4.62 ns / 8.33$ (1.63)	$2.61 ns / 46.18$ (1.20)	$1.33 ns / 13.96$ , $3.43 ns / 20.67$ (1.14)
(170, 130) Photon count: 21	$2.87 ns / 15.38$ (1.07)	$4.41 ns / 16.39$ , $9.62 ns / 47.51$ (0.80)	$1.70 ns / 10.64$ , $4.39 ns / 9.29$ (1.05)
(180, 130) Photon count: 19	$1.39 ns / 15.75$ , $5.07 ns / 4.47$ (1.64)	$2.50 ns / 19.61$ , $5.91 ns / 17.66$ (0.87)	$1.33 ns / 5.73$ , $3.55 ns / 13.35$ (1.13)

Figure 3 shows the FLIM images and distributions of lifetimes without and with an acceptor in a HeLa cell estimated by SPC software and the proposed method. Figures 3(a) and 3(c) are FITC- $A β 11$ -25 donor FLIM images analyzed for the same TCSPC photon counts: Fig. 3(a) was obtained using SPC software and Fig 3(c) from the proposed method. Similarly, Figs. 3(e) and 3(g) are FITC- $A β 11$ -25-DABMI donor–acceptor FLIM images analyzed by the two algorithms, where Fig. 3(e) is from the SPC software and Fig. 3(g) is from the proposed method. To generate these FLIM images with multiple lifetimes, we used the intensity weighted mean lifetime at each pixel and it is coded with pseudo color. The left column shows FLIM images and the right column corresponds to lifetime distributions. As stated above, the pseudo color of the FLIM image in Figs. 3(e) and 3(g), which show shortened lifetimes due to FRET, is more yellowish than Figs. 3(a) and 3(c), respectively. SPC software shows a single lifetime parameter: it shows a peak at 2.5 ns in Fig. 3(b) and the peak position changes to 1.8 ns in Fig. 3(f) after FRET. The same TCSPC data are analyzed by SPIRAL: it estimates two distinctive peaks at 1.8 and 3.4 ns as shown in Fig. 3(d). The amplitude of the 1.8 ns component increases after FRET, while the amplitude of the 3.4 ns component decreases, as observed in Fig. 3(h). The distribution of the lifetime after FRET shows that energy transfer dominantly occurs at 1.8 ns. Since a HeLa cell has biological structure, it cannot be uniformly stained. Moreover, the reaction of FITC differs depending on the local environment such as pH distribution or strain. We believe that the spatial nonuniformity that resulted from the proposed method shows the spatial variation of the structure and the local environment of the cell, while the conventional SPC method cannot resolve the fine details. Thus, the ability to show the fine structure of the HeLa cell reflects an improvement in the quality of the proposed imaging technique. Also please note that the SPIRAL analysis of Figs. 3(c) and 3(g) delineates cell boundary much better than the corresponding images of Figs. 3(a) and 3(e) from the SPC software. Although Fig. 3(c) is the image composed by the data obtained before FRET occurs, the data still have double lifetimes due to the characteristics of FITC and the lifetime in each local region varies significantly depending on pH concentration.¹ Therefore, local homogeneous regions in Fig. 3(c) imply similar pH concentrations in the region. Since Fig. 3(g) is an image after FRET occurs, Fig. 3(g) shows the homogeneous regions with blue or yellow pseudo colors. The bluish region is where FRET does not occur, while the yellowish region represents where FRET occurs, since the region with FRET has shorter lifetimes than the region without FRET. Reference 1 supports this phenomenon. The experimental results demonstrate the applicability of the proposed method to the multiexponential decay of fluorophores from low count cell data, where the maximum count is around 40 photons.

Fig. 3

Comparison of fluorescence lifetime imaging microscopy (FLIM) images and the distributions of intensity weighted lifetimes between the proposed method and SPC software (a) FITC- $A β$ FLIM image for HeLa cell using SPC software. (b) The corresponding distribution of lifetime for (a). (c) The same FITC- $A β$ FLIM image analyzed by sparse Poisson intensity reconstruction algorithm (SPIRAL). (d) The corresponding distribution of lifetime for (c). (e) FITC- $A β$ -DABMI FLIM image for HeLa cell using SPC software. (f) The corresponding distribution of lifetime for (e). (g) The same FITC- $A β$ -DABMI FLIM image analyzed by SPIRAL. (h) The corresponding distribution of lifetime for (g).

6. Conclusion

We presented a powerful method to estimate the distribution of fluorescence lifetime using Poisson noise modeling and $L 1$ -norm regularization. We use the sensing matrix $A$ based on the multiexponential decay of fluorophores and the instrument response function, together with $L 1$ -norm regularization for compressive sensing.

We verified that the presented algorithm shows a more robust estimation of lifetime than the maximum-likelihood estimator based on a Poisson noise from computer simulations and a more accurate distribution with real cell with low photon count images with FRET. SPIRAL has the merit of showing the width or distribution of lifetime components. We need extensive work for the theoretical analysis and practical evaluation of SPIRAL. For future work, we can incorporate a spatial prior in order to improve the quality of FLIM images.

Acknowledgments

This work was supported in part by the Basic Science Research Program through the National Research Foundation of Korea (NRF) funded by the Ministry of Science, ICT & Future Planning (NRF-2012R1A1A3013525) and by the Bio and Medical Technology Development Program through the National Research Foundation of Korea (NRF) funded by the Ministry of Education, Science and Technology (MEST) (No. 2012M3A9B6055379).

References

1.

J. R. Lakowicz, Principles of Fluorescence Spectroscopy, Springer, New York (2006). Google Scholar

2.

B. B. Collier and M. J. McShane, “Dynamic windowing algorithm for the fast and accurate determination of luminescence lifetimes,” Anal. Chem., 84 4725 –4731 (2012). http://dx.doi.org/10.1021/ac300023q ANCHAM 0003-2700 Google Scholar

3.

M. Wahl, Technical Note on Time-correlated Single Photon Counting, PicoQuant GmbH, Berlin (2009). Google Scholar

4.

Y. Won et al., “High-speed confocal fluorescence lifetime imaging microscopy (FLIM) with the analog mean delay (AMD) method,” Opt. Express, 19 (4), 3396 –3404 (2011). http://dx.doi.org/10.1364/OE.19.003396 OPEXFF 1094-4087 Google Scholar

5.

M. C. Skala et al., “In vivo multiphoton microscopy of NADH and FAD redox states, fluorescence lifetimes, and cellular morphology in precancerous epithelia,” Proc. Natl. Acad. Sci. U.S.A., 104 (49), 19494 –19499 (2007). http://dx.doi.org/10.1073/pnas.0708425104 Google Scholar

6.

M. S. Roberts et al., “Non-invasive imaging of skin physiology and percutaneous penetration using fluorescence spectral and lifetime imaging with multiphoton and confocal microscopy,” Eur. J. Pharm. Biopharm., 77 469 –488 (2011). http://dx.doi.org/10.1016/j.ejpb.2010.12.023 EJPBEL 0939-6411 Google Scholar

7.

G. Deka, W. W. Wu and F. J. Kao, “In vivo wound healing diagnosis with second harmonic and fluorescence lifetime imaging,” J. Biomed. Opt., 18 (6), 061222 (2013). http://dx.doi.org/10.1117/1.JBO.18.6.061222 JBOPFO 1083-3668 Google Scholar

8.

D. V. O’Connor and D. Phillips, Time-Correlated Single Photon Counting, Academic Press, London (1984). Google Scholar

9.

A. E. McKinnon, A. G. Szabo and D. R. Millar, “The deconvolution of photoluminescence data,” J. Phys. Chem., 81 (16), 1564 –1570 (1977). http://dx.doi.org/10.1021/j100531a009 JPCHAX 0022-3654 Google Scholar

10.

D. V. O’Connor, W. R. Ware and J. C. Andre, “Deconvolution of fluorescence decay curves. A critical comparison of techniques,” J. Phys. Chem., 83 1333 –1343 (1979). http://dx.doi.org/10.1021/j100473a019 JPCHAX 0022-3654 Google Scholar

11.

J. R. Knutson, J. M. Beechem and L. Brand, “Simultaneous analysis of multiple fluorescence decay curves: a global approach,” Chem. Phys. Lett., 102 (6), 501 –507 (1983). http://dx.doi.org/10.1016/0009-2614(83)87454-5 CHPLBC 0009-2614 Google Scholar

12.

T. Hauschild and M. Jentschel, “Comparison of maximum likelihood estimation and chi-square statistics applied to counting experiments,” Nucl. Instrum. Methods Phys. Res., Sect. A, 457 384 –401 (2001). http://dx.doi.org/10.1016/S0168-9002(00)00756-7 NIMAER 0168-9002 Google Scholar

13.

K. Sasaki and H. Masuhara, “Analysis of transient emission curves by a convolved autoregressive model,” Appl. Opt., 30 (8), 977 –980 (1991). http://dx.doi.org/10.1364/AO.30.000977 Google Scholar

14.

V. V. Apanasovich and E. G. Novikov, “Identification of lifetime characteristics in fluorescence experiments,” Proc. SPIE, 2388 378 –384 (1995). http://dx.doi.org/10.1117/12.208496 PSISDG 0277-786X Google Scholar

15.

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

16.

J. Kim and J. Seok, “Statistical properties of amplitude and decay parameter estimators for fluorescence lifetime imaging,” Opt. Express, 21 (5), 6061 –6074 (2013). http://dx.doi.org/10.1364/OE.21.006061 OPEXFF 1094-4087 Google Scholar

17.

Z. T. Harmany, R. F. Marcia and R. M. Willett, “This is SPIRAL-TAP: sparse Poisson intensity reconstruction algorithms-theory and practice,” IEEE Trans. Image Process., 21 (3), 1084 –1096 (2012). http://dx.doi.org/10.1109/TIP.2011.2168410 IIPRE4 1057-7149 Google Scholar

18.

The bh TCSPC Handbook, 5th ed.Becker & Hickl GmbH, Berlin (2012). Google Scholar

19.

D. L. Donoho, “Compressed sensing,” IEEE Trans. Inf. Theory, 52 (4), 1289 –1306 (2006). http://dx.doi.org/10.1109/TIT.2006.871582 IETTAW 0018-9448 Google Scholar

20.

E. J. Candes and T. Tao, “Decoding by linear programming,” IEEE Trans. Inf. Theory, 51 (12), 4203 –4215 (2006). http://dx.doi.org/10.1109/TIT.2005.858979 IETTAW 0018-9448 Google Scholar

21.

E. Candes, J. Romberg and T. Tao, “Near-optimal signal recovery from random projection: universal encoding from highly incomplete frequency information,” IEEE Trans. Inf. Theory, 52 (12), 5406 –5425 (2006). http://dx.doi.org/10.1109/TIT.2006.885507 IETTAW 0018-9448 Google Scholar

22.

R. Willett and M. Raginsky, “Performance bounds for compressed sensing with Poisson noise,” in Proc. IEEE Int. Symp. Information Theory, 174 –178 (2009). Google Scholar

23.

M. Raginsky and Z. T. Harmany, “Compressed sensing performance bounds under Poisson noise,” IEEE Trans. Signal Process., 58 (8), 3990 –4002 (2010). http://dx.doi.org/10.1109/TSP.2010.2049997 ITPRED 1053-587X Google Scholar

24.

N. S. Caron et al., “Using FLIM-FRET to measure conformational changes of transglutaminase type 2 in live cells,” PLoS One, 7 (8), e44159 (2012). http://dx.doi.org/10.1371/journal.pone.0044159 POLNCL 1932-6203 Google Scholar

25.

B. Valeur and M. N. Berberan-Santos, Molecular Fluorescence Principles and Applications, Wiley-VCH, (2012). Google Scholar

26.

C. A. Swindlehurst and Jr. E. W. Voss, “Fluorescence measurements of immune complexes of Mab 4-4-20 with isomeric haptens,” Biophys. J., 59 619 –628 (1991). http://dx.doi.org/10.1016/S0006-3495(91)82277-9 BIOJAU 0006-3495 Google Scholar

27.

S. Santra et al., “Fluorescence lifetime measurements to determine the core-shell nanostructure of FITC-doped silica nanoparticles: an optical approach to evaluate nanoparticle photostability,” J. Lumin., 117 75 –82 (2006). http://dx.doi.org/10.1016/j.jlumin.2005.04.008 JLUMA8 0022-2313 Google Scholar

Biography

Sejung Yang is a research assistant professor at the Ewha Womans University. She received her BS degree in electrical engineering from the Ewha Womans University in 2001, her MS degree in electrical engineering from the University of Southern California in 2005, and her PhD in electronics engineering from the Ewha Womans University in 2011. Her current research interests include image processing and biomedical imaging.

Joohyun Lee received her BS degree in electrical engineering from Ewha Womans University in 2010 and her MS degree in 2012. Currently, she is an associate research engineer at KT Institute of convergence technology. Her research interests include color image enhancement, denoising, and medical image processing.

Youmin Lee received her bachelor’s degree from the Department of Chemistry and Nanoscience, Ewha Womans University in South Korea in 2012, where she worked toward a master’s degree in the physical chemistry group. She received her MS degree in 2014. Her research was focused on investigating photophysics of fluorescence probes in microenvironments using time-resolved fluorescence.

Minyung Lee has been a professor in the Department of Chemistry and Nanoscience at Ewha Womans University. He completed his doctorate in physical chemistry at the University of Pennsylvania in 1987. He developed a couple of measurement techniques incorporating atomic force microscopy and fluorescence lifetime imaging microscopy (FLIM). His major research field has been biomolecular dynamics and current research interest includes investigations on the structural dynamics and aggregation processes of peptides and proteins, employing FLIM.

Byung-Uk Lee received his BS degree from Seoul National University in 1979, his MS degree from the KAIST in 1981, and his PhD in electrical engineering from Stanford University in 1991. He worked for Daewoo Electronics Co. from 1991 to 1995. He joined the Department of Electronics Engineering, Ewha Womans University in 1995, where he is currently a professor. His research interests include color image enhancement, denoising, and medical image processing.

CC BY: © The Authors. Published by SPIE under a Creative Commons Attribution 4.0 Unported License. Distribution or reproduction of this work in whole or in part requires full attribution of the original publication, including its DOI.

Citation Download Citation

Sejung Yang, Joohyun Lee, Youmin Lee, Minyung Lee, and Byung-Uk Lee "Estimation of multiexponential fluorescence decay parameters using compressive sensing," Journal of Biomedical Optics 20(9), 096003 (3 September 2015). https://doi.org/10.1117/1.JBO.20.9.096003

Published: 3 September 2015

Access the abstract

JOURNAL ARTICLE
9 PAGES

DOWNLOAD PAPER SAVE TO MY LIBRARY

GET CITATION

CITATIONS

Cited by 7 scholarly publications.

Explore citations on Lens.org

KEYWORDS

Fluorescence lifetime imaging

Photon counting

Luminescence

Fluorescence resonance energy transfer

Image analysis

Atrial fibrillation

Biological research

1.

Introduction

2.

Prior Work

2.1.

Modified Neyman’s Nonlinear Least-Squares Method

Eq. (1)

Eq. (2)

Eq. (3)

Eq. (4)

2.2.

Fast Fitting Using a Convolved Autoregressive Model

Eq. (5)

Eq. (6)

Eq. (7)

Eq. (8)

Eq. (9)

Eq. (10)

Eq. (11)

Eq. (12)

Eq. (13)

Eq. (14)

Eq. (15)

3.

Proposed Method

3.1.

Problem Formulation

Eq. (16)

3.2.

Poisson Log-Likelihood

Eq. (17)

Eq. (18)

Eq. (19)

3.3.

Objective Function

Eq. (20)

4.

Numerical Simulations

Eq. (21)

Eq. (22)

Table 1

5.

Experimental Results

Eq. (23)

Table 2

Fig. 1

Fig. 2

Table 3

Fig. 3

6.

Conclusion

Acknowledgments

References

Biography

Show All Keywords

Keywords/Phrases

Search In:

Publication Years