Penalized maximum likelihood estimation of lifetime and amplitude images from multi-exponentially decaying fluorescence signals

Jeongtae Kim; Jiyeong Seok; Hwiin Lee; Minyung Lee

doi:10.1364/OE.21.020240

1. Introduction

Fluorescence lifetime imaging microscopy (FLIM) is based on the excited-state decay properties of fluorophores and is useful for studying molecular environments in living cells and human tissues [1, 2] because the lifetime of a fluorophore is dependent on the molecular environment [3]. FLIM has been used to study, for example, the intracellular pH distribution from an estimated lifetime map [3,4] and fluorescence resonance energy transfer in a cell [5]. Since the lifetime is an indicator of the molecular environment, an accurate and precise estimation of the lifetime is crucial for applying FLIM to biological studies. The fluorescence lifetime can be estimated in the time domain using a time-correlated single photon counting (TCSPC) technique that fits an exponential decay model to measured photon counts by maximizing a similarity measure such as χ² or the Poisson likelihood function [1]. It was reported that maximum likelihood estimation (MLE) based on maximizing the Poisson likelihood function outperforms other methods for TCSPC-based FLIM [6–8]. Alternatively, the lifetime can be estimated in the frequency domain using intensity-modulated light, the associated phase shifts and signal demodulation [1]. For both methods, however, the precision of the estimated lifetime map is limited by the signal-to-noise ratio (SNR) of the measured photon counts [5], which suffer from Poisson noise. Note that the precision of unbiased lifetime estimation becomes worse as SNR decreases since the variance of the estimation is bounded by the Cramer-Rao bound, which becomes higher as SNR becomes lower [9–12]. As the SNR is proportional to the total number of photon counts, increasing the exposure time of the fluorophores to the laser source may improve the SNR [5] but since the exposure time is usually limited due to photobleaching and imaging speed, such an attempt may not be successful [5].

To improve the performance of the lifetime estimation without increasing the exposure time, global analysis-based methods have been studied [5, 13]. A central assumption is that the lifetimes of all the fluorophores in an image are the same [5, 13]. Based on this assumption, data sets from different locations are analyzed simultaneously to estimate the lifetime, while the fluorophore concentration at each pixel location is allowed to be different [5, 13, 14]. The method is successful in improving the accuracy of the lifetime estimation provided the assumption is valid, but this is not always true; for example, when the lifetime is an indicator of pH, the lifetime at each pixel in an image will not necessarily be the same [3, 15].

A different approach is based on image restoration techniques such as denoising and deconvolution [2, 16, 17]. Since the variance in the lifetime estimation becomes larger as the SNR of the fluorescence decreases, an estimated lifetime map from a low-SNR signal is usually noisy. Therefore, denoising of the estimated lifetime image may decrease the variance in the lifetime estimation. Previous investigations have applied the iterative constrained Tikhonov–Miller algorithm to restore the Fourier coefficients in a frequency-domain lifetime estimation [16] and total variation (TV)-based denoising methods for time-gated FLIM [2,17]. However, the former method was applied to frequency domain FLIM, while the latter method was investigated for time-gated FLIM. These methods are also unable to handle data that should be excluded from the lifetime estimation, i.e., the background signal and dark pixels, where the total number of photon counts is very small [18]. To the best of our knowledge, no investigations applying image restoration techniques to TCSPC-based lifetime estimation have been published, even though such investigations are essential as TCSPC-based lifetime estimation performs better than any other method, especially for low photon counts (i.e., low SNR cases) [3, 19, 20].

In this paper, we propose a penalized maximum likelihood estimation (PMLE) of the lifetimes and amplitudes for TCSPC-based FLIM. To overcome the noise problem of independent MLE of the lifetimes and amplitudes at each pixel, the proposed method estimates the lifetimes and amplitudes by maximizing an objective function that is the weighted sum of a negative Poisson likelihood function and a penalty function that penalizes large variations in the lifetimes and amplitudes in the spatial domain. A similar PMLE method has been shown to be effective for reducing noise in medical image reconstructions [21–24]. We calculate the variations in the lifetimes and amplitudes after excluding pixels for which the number of photon counts is below a certain threshold. We also introduce a novel parameterization to determine the lifetimes uniquely. Since the proposed objective function is not an additively separable function of parameters at each pixel, the optimization of the proposed objective function requires simultaneous optimization of a significantly large size parameter vector (i.e., lifetimes and amplitudes from every pixel), which is often intractable. To overcome this, we apply the principle of optimization transfer, which is based on a series of optimizations of a separable surrogate function [21–23] for the penalty function. It has been shown that the series of optimizations of the surrogate functions eventually converges to the minimizer of the original objective function [21–23]. We demonstrate here the effectiveness of the proposed method in comparison with the MLE through simulations and experiments.

2. Theory

2.1. Problem formulation

We model the number of photon counts detected at a location (x_i, y_j) at a time instance t_k by the Poisson random variable, g(t_k, x_i, y_j), as follows [1]:

g (t_{k}, x_{i}, y_{j}) = Poisson {λ (t_{k}; θ_{i, j})}, k = 0, 1, ., K - 1,

where the mean of the Poisson random variable is defined as

λ (t_{k}; θ_{i, j}) = h (t_{k}) * f (t_{k}; θ_{i, j}) .

The impulse response function (IRF) of the measurement system is h(t), and f(t) is the multi-exponential decay function at (x_i, y_j) defined by

f (t; θ_{i, j}) = \sum_{p = 0}^{P - 1} A_{p} (x_{i}, y_{j}) e^{- \frac{t}{τ_{p} (x_{i}, y_{j})}},

where τ_p(x_i, y_j) is the p-th decay constant at (x_i, y_j), and A_p(x_i, y_j) is the associated amplitude. We define the parameter vector at (x_i, y_j) as

θ_{i, j} = [A_{0} (x_{i}, y_{j}), A_{1} (x_{i}, y_{j}), . ., A_{P - 1} (x_{i}, y_{j}), τ_{0} (x_{i}, y_{j}), Δ τ_{1} (x_{i}, y_{j}), . ., Δ τ_{P - 1} (x_{i}, y_{j})] .

Using the parameter vector defined in Eq. (4), the p-th decay time constant can be represented by

τ_{p} (x_{i}, y_{j}) = τ_{0} (x_{i}, y_{j}) + \sum_{p = 1}^{P - 1} Δ τ_{p} (x_{i}, y_{j}), p = 1, 2, . ., P - 1 .

We use this parameterization to determine each decay time constant uniquely. The goal of FLIM imaging is to obtain an accurate estimation of the decay and amplitude constants θ_i,j at every location using the collected number of detected photon counts for an M × N image, g = {g_i,j|i = 0,,N − 1, j = 1,,M − 1}, where g_i,j = [g(t₀, x_i, y_j),., g(t_K−1, x_i, y_j)]. Because it is not possible to determine the lifetime accurately when the total number of detected photon counts is too small, we do not attempt to estimate lifetimes at such locations. In other words, we attempt to estimate the following parameter vector,

Θ = {θ_{i, j} | (i, j) \in C},

where C is a set that defines the indices of the pixel locations where the number of detected photon counts is larger than some threshold γ. The log-likelihood function of the entire parameter vector Θ from the entire measurement data g is defined as

L (Θ; g) = \sum_{(i, j) \in C} L_{i, j} (θ_{i, j}; g_{i j}),

where the log-likelihood function of the parameter vector θ_i,j from the measurement g_ij is

L_{i, j} (θ_{i, j}; g_{i j}) = \sum_{k = 0}^{K - 1} - λ (t_{k}; θ_{i, j}) + g (t_{k}, x_{i}, y_{j}) log λ (t_{k}; θ_{i, j}) - log g (t_{k}, x_{i}, y_{j})! .

The MLE of Θ from the measured data g is given by

\hat{Θ} = \underset{Θ \geq 0}{argmax} L (Θ; g) .

Note that the size of the parameter vector is 2P × 𝒞, where 𝒞 s the cardinality of the set C. Since L(Θ;g) is additively separable, the MLE of the parameter vector θ_i,j for each pixel location can be determined independently from other parameter vectors:

{\hat{θ}}_{i, j} = \underset{θ_{i, j} > 0}{argmax} L_{i, j} (θ_{i, j}; g_{i, j}), (i, j) \in C .

One problem of using such an MLE is the noise sensitivity, since the method does not consider spatial correlations between neighboring decay and amplitude constants.

2.2. Proposed method

To reduce the noise in the estimated lifetime and amplitude images using the MLE method, we propose a method that incorporates a priori knowledge: the amplitudes and lifetimes do not change rapidly in the spatial domain. The objective function of the proposed method Φ(Θ;g) consists of the weighted sum of the Poisson likelihood function and a penalty function R(Θ) that discourages large variations in the spatial domain:

Φ (Θ; g) = - L (Θ; g) + β R (Θ),

where β is a regularization parameter that determines the weight between the first data fidelity term and second penalty term. The penalty function measures the total roughness of the lifetime and associated amplitude maps and is defined as

R (Θ) = \sum_{p = 0}^{2 P - 1} \sum_{(i, j) \in C} ψ ({[D_{i j}^{h} Θ]}_{p}) + ψ ({[D_{i j}^{v} Θ]}_{p}),

where the change of the p-th component at the (i, j)-th pixel location in the horizontal direction is measured by

ψ ({[D_{i j}^{h} Θ]}_{p}) = {\begin{array}{l} ϕ ({[θ_{i + 1, j} - θ_{i, j}]}_{p}), & if & (i + 1, j) \in C \\ 0, & otherwise, \end{array}

and that in the vertical direction is measured by

ψ ({[D_{i j}^{v} Θ]}_{p}) = {\begin{array}{l} ϕ ({[θ_{i, j + 1} - θ_{i, j}]}_{p}), & if & (i, j + 1) \in C \\ 0, & otherwise, \end{array}

where the 2P × (2P × 𝒞) matrices

D_{i j}^{h}

and

D_{i j}^{v}

are horizontal and vertical parameter difference calculating matrices at the (i, j)-th pixel, respectively. A natural choice for φ (·) is a quadratic function; however, it is well known that quadratic penalty functions over-smooth edges in restored images [23]. We instead adopt an edge preserving penalty function such as the TV penalty function, which is defined as follows [23]:

ϕ (x) = \sqrt{x^{2} + ε^{2}},

where ε > 0 is a small positive constant to ensure that the penalty function is differentiable at x = 0. Since the functions defined in Eq. (13) and (14) are not additively separable functions of θ_i,j, minimization of the objective function defined in Eq. (11) requires the simultaneous optimization of the 2P × 𝒞 -size parameter vector Θ, which is often intractable. To overcome this difficulty, we apply an optimization transfer approach to determine the minimizer of the objective function. For the minimization of an objective function Φ(θ), the iterative minimization of a surrogate function Φ^s(θ; θ^k) (θ^k denotes the estimated value of θ at the k-th iteration) monotonically converges to the minimizer of the objective function Φ(θ) if the surrogate function satisfies the following two conditions [23]:

Φ^{s} (θ; θ^{k}) \geq Φ (θ),

Φ^{s} (θ^{k}; θ^{k}) = Φ (θ^{k}) .

That is, the following sequence θ^k monotonically converges to the minimizer of the objective function Φ(θ):

θ^{k + 1} = \underset{θ}{argmin} Φ^{s} (θ; θ^{k}) .

We design additively separable quadratic surrogate functions for the penalty functions defined in Eq. (13) and (14) so that we can estimate the parameters at each pixel independently. To do this, we first apply Huber’s method [25] to design the quadratic surrogate functions of the penalty functions and then apply De Pierre’s approach [24] to construct the quadratic surrogate functions separately as follows [21–23]:

\begin{array}{l} ϕ ({[θ_{i + 1, j} - θ_{i, j}]}_{p}) & \leq & q ({[D_{i j}^{h} Θ]}_{p}; {[D_{i j}^{h} Θ^{k}]}_{p}) \\ = & q (\frac{1}{2} {2 {[θ_{i + 1, j} - θ_{i + 1, j}^{k}]}_{p}} + \frac{1}{2} {- 2 {[θ_{i, j} - θ_{i, j}^{k}]}_{p}} + {[D_{i j}^{h} Θ^{k}]}_{p}; {[D_{i j}^{h} Θ^{k}]}_{p}) \\ \leq & \frac{1}{2} q (2 {[θ_{i + 1, j} - θ_{i + 1, j}^{k}]}_{p} + {[D_{i j}^{h} Θ^{k}]}_{p}; {[D_{i j}^{h} Θ^{k}]}_{p}) \\ + \frac{1}{2} q (- 2 {[θ_{i, j} - θ_{i, j}^{k}]}_{p} + {[D_{i j}^{h} Θ^{k}]}_{p}; {[D_{i j}^{h} Θ^{k}]}_{p}), \end{array}

where

q (a; b) = ϕ (b) + \dot{ϕ} (b) (a - b) + \frac{1}{2} c {(a - b)}^{2}

and φ̇ is the derivative of φ [23]. The curvature c that guarantees the first inequality of a quadratic or TV penalty function can be determined using Huber’s method [25]. The last inequality in Eq. (19) is due to the convex inequality [24]. Therefore, for the p-th component, the following inequality holds:

\sum_{(i, j) \in C} ψ ({[D_{i j}^{h} Θ]}_{p}) \leq \sum_{(i, j) \in C} R_{h} ({[θ_{i, j}]}_{p}; Θ_{p}^{k}),

where,

R_{h} ({[θ_{i, j}]}_{p}; Θ^{k}) = {\begin{array}{l} \frac{1}{2} q (2 {[θ_{i, j} - θ_{i, j}^{k}]}_{p} + {[D_{i - 1, j}^{h} Θ^{k}]}_{p}; {[D_{i - 1, j}^{h} Θ^{k}]}_{p}) \\ + \frac{1}{2} q (- 2 {[θ_{i, j} - θ_{i, j}^{k}]}_{p} + {[D_{i, j}^{h} Θ^{k}]}_{p}; {[D_{i, j}^{h} Θ^{k}]}_{p}), & (i - 1, j) \in C, (i + 1, j) \in C \\ \frac{1}{2} q (2 {[θ_{i, j} - θ_{i, j}^{k}]}_{p} + {[D_{i - 1, j}^{h} Θ^{k}]}_{p}; {[D_{i - 1, j}^{h} Θ^{k}]}_{p}), & (i - 1, j) \in C, (i + 1, j) \notin C \\ \frac{1}{2} q (- 2 {[θ_{i, j} - θ_{i, j}^{k}]}_{p} + {[D_{i, j}^{h} Θ^{k}]}_{p}; {[D_{i, j}^{h} Θ^{k}]}_{p}), & (i - 1, j) \notin C, (i + 1, j) \in C \\ 0, & (i + 1, j) \notin C, (i - 1, j) \notin C . \end{array}

It is relatively straightforward to design a similar surrogate function R_v([θ_i,j]_p; Θ^k) for

ψ ({[D_{j}^{v} Θ]}_{p}) = ψ ({[θ_{i, j + 1} - θ_{i, j}]}_{p})

. Based on Eq. (20), the separable quadratic surrogate function of the penalty function, R_s(Θ), satisfies the following inequality:

R (Θ) \leq R_{s} (Θ; Θ^{k}) = \sum_{(i, j \in C)} R_{h v} (θ_{i, j}; Θ^{k}),

where

R_{h v} (θ_{i, j}; Θ^{k}) = \sum_{p = 0}^{2 P - 1} R_{h} ({[θ_{i, j}]}_{p}; Θ^{k}) + R_{v} ({[θ_{i, j}]}_{p}; Θ^{k}) .

Note that R(Θ^k) = R_s(Θ^k;Θ^k). Since both the likelihood function and the surrogate function of the penalty function are additively separable, the minimizer of the weighted sum of the log-likelihood function and the additively separable surrogate function at the k-th iteration can be determined at each pixel location independently as follows:

θ_{i, j}^{k + 1} = \underset{θ_{i, j} > 0}{argmin} - L_{i j} (θ_{i, j}; g_{i, j}) + β R_{s} (θ_{i, j}; Θ^{k}) .

The minimization problem defined in Eq. (23) can be solved numerically. The minimizer sequence defined in Eq. (23) monotonically converges to the minimizer of the objective function defined in Eq. (11).

3. Results

To evaluate the performance of the proposed PMLE method and compare it with the conventional MLE method, we conducted simulation studies and real experiments. We implemented the MLE and the proposed PMLE method using the MATLAB (Mathworks, USA). In addition, we solve the minimization problems in Eq. (10) for the MLE and in Eq. (23) for the PMLE method using the fmincon function in the Optimization toolbox of the MATLAB. We run the methods on a workstation equipped with two Intel Xeon X5650 processors (2.67 GHz) and 96 GB memory.

3.1. Simulations

We conducted simulation studies using a 64×64 synthetic image. The synthetic image has an elliptical structure that contains a square and a circle. We synthesized a double exponentially decaying fluorescence signal at each pixel location that has different decay constants (τ₁ and τ₂) and amplitudes (A₁ and A₂) in different regions of the image and convolved these synthesized signals with an IRF function of [0.014, 0.042, 0.121, 0.273, 0.357, 0.194]. We then generated 261 Poisson random variables at each pixel, whose means were uniformly sampled values of the convolved decaying fluorescence signal. The true values of the lifetime, amplitude and the average photon counts in each region are summarized in Table 1.

Table 1. True values of the lifetimes, amplitudes and the average photon counts in the synthetic image (the lifetime is expressed in units of nanoseconds)

View Table | View all tables in this article

For estimating the lifetime and amplitude images from the synthesized noisy fluorescence signals (i.e., the Poisson random variables), we applied the MLE, the PMLE with a quadratic penalty function (PMLEQ), or the PMLE with a TV penalty function (PMLETV). Regularization parameters for (τ₁, τ₂, A₁, A₂) were (250, 100, 0.5, 1.5) in PMLEQ method and (8, 6, 0.4, 0.5) in PMLETV method. We also applied the MLE method with binning of the measured photon counts (i.e., summing the photon counts of neighboring pixels) from 3 × 3 neighboring pixels, which is used to reduce the effect of noise in the commercial software SPCImage (Becker & Hickl GmbH). The parameters Fig. 1 shows the estimated lifetime and amplitude images using the MLE method, which as expected yielded noisy lifetime and amplitude images due to Poisson noise. To reduce the effect of Poisson noise, we applied the PMLEQ and PMLETV methods. Figure 2 shows the estimated images obtained from the PMLEQ method, and as can be seen, the PMLEQ method was able to reduce the effect of noise. However, the quadratic penalty function blurred the edges in the estimated images. The blurring can be reduced by applying the PMLETV method, as shown in Fig. 3. Figure 4 shows estimated images using the MLE method with binning. As shown in the figure, although the binning strategy was able to reduce the effect of noise, it yielded images with artifacts, especially in the region where lifetimes in neighboring pixels are different. This is due to that photon counts data from different lifetimes are summed together.

Fig. 1 Estimated lifetime and amplitude images obtained with MLE: (a) τ₁, (b) τ₂, (c) A₁, and (d) A₂.

Download Full Size | PDF

Fig. 2 Estimated lifetime and amplitude images obtained with PMLEQ: (a) τ₁, (b) τ₂, (c) A₁, and (d) A₂.

Download Full Size | PDF

Fig. 3 Estimated lifetime and amplitude images obtained with PMLETV: (a) τ₁, (b) τ₂, (c) A₁, and (d) A₂.

Download Full Size | PDF

Fig. 4 Estimated lifetime and amplitude images obtained using MLE with binning of 3×3 window: (a) τ₁, (b) τ₂, (c) A₁, and (d) A₂.

Download Full Size | PDF

We repeated the estimation of the lifetime and amplitude images using the four methods for 10 times with different noise realizations and computed the bias and standard deviation (STD) of the lifetimes and amplitudes for each region using the pixel values in all the estimated images. For the MLE method with binning, we normalized the bias and the standard deviation by the number of pixels in the binning window. We summarize the biases and standard deviations of the lifetimes and amplitudes in each region in Table 2. The biases of the PMLETV method were quite similar to those of the MLE method (slightly larger in estimating amplitudes), whereas the standard deviations were significantly lower than the MLE method, indicating that the effect of random noise is greatly reduced by the PMLETV method. While the PMLEQ method also yielded lower standard deviations than the MLE method, the standard deviations of the PMLEQ method were larger than those of the PMLETV method. The standard deviations of the MLE method with binning were larger than those of the PMLE methods. Note that as the effect of penalty function increases, noise can be reduced at the price of resolution. The regularization parameter controls trade-off between noise and resolution [26]. Therefore, one should design the regularization parameter based on desired resolution and noise property. Automatic determination of the parameter is a challenging task, and we defer it to future study.

Table 2. Bias and STD of the estimated lifetimes and amplitudes in the ellipse, square, and circle

View Table | View all tables in this article

We also evaluated the similarity between the estimated and true images by computing the correlation coefficients of the four estimated images (τ₁, τ₂, A₁, A₂). Table 3 summarizes the results and reveals that the PMLE methods yield images that are more similar to the true images than those images obtained with the MLE method and the MLE method with binning. Among the two PMLE methods, the PMLETV method yielded better results than the PMLEQ method.

Table 3. Correlation between the true and estimated images

View Table | View all tables in this article

The disadvantage of the proposed method is longer computation time. We summarize the average computation time for estimating the lifetime and amplitude maps using the MLE, PM-LEQ and PMLETV methods in Table 4. Since we initialized lifetime and amplitude maps using the ones estimated by the MLE method, one should think that total computation times for the PMLE methods are approximately the same as the sums of the computation time for the MLE and for the PMLE methods. We expect that computation time for the PMLE methods can be reduced by using graphics processing unit and parallel computing technology as well as by the improvement of computer technology. Investigations on faster PMLE methods are deferred to future study.

Table 4. Average computation time in seconds

View Table | View all tables in this article

3.2. Real data

3.2.1. polymer data

To prove the effectiveness of the proposed method for real data, we conducted experiments using a polymer sample which contains two different dyes that have different lifetimes. The sample was prepared as follows. Coumarin 6 (C6) and 4-(4-methoxybenzylamino)-7-nitrobenzofurazan (MBD), tetrahydrofuran (THF), poly(methyl methacrylate) (PMMA, MW 120,000) and polyvinyl chloride (PVC, MW 97,000) were purchased from Sigma-Aldrich. C6 and PVC were dissolved in THF and MBD and PMMA were dissolved in THF. Each solution was dropped at contact on a cover glass and this sample was dried in an oven at 80°C for 24 hours. The light sources were two picosecond diode lasers operating at the 442 nm wavelength with the repetition rate of 20 MHz (Picoquant). A platform for sample excitation and fluorescence detection was an inverted confocal microscope (Nikon, TE2000-S) with an oil immersion objective lens (NA 1.4, ×60). A long pass cut-off filter was placed in front of the detector to block any remaining excitation beam. The total fluorescence was detected by a micro channel plate photomultiplier tube (Hamamatsu R3809U-51). The signal was amplified by a GHz preamplifier (EG&G 9306) and the signal was processed by a fast time-correlated single photon counting (TCSPC) board (Becker-Hickl, SPC-830). The x–y images consist of arrays of pixels (256 × 256), each of which contains the fluorescence decay curve. We obtained low and high photon counts by controlling the excitation laser intensities. Using the fluorescence decay curves from 80×80 sub-image region, we attempted to estimate lifetime and amplitude images. Figure 5(a) shows estimated lifetime image from high photon counts data (average number of photons per each pixel is 2686) using the MLE method while Fig. 5(b) does from low photon counts data (average number of photons per each pixel is 335), which is noisier than Fig. 5(a). The short and long lifetime components correspond to C6 and MBD, respectively. To acquire less noisy lifetime image, we applied the PMLEQ and the PMLETV method to the low photon counts data. The results are shown in Fig. 5(c) and Fig. 5(d). As shown in the figures, the estimated lifetime images using the PMLE methods look less noisier than the image using the MLE method. The correlation coefficients between the image from the high photon counts data (which serves as a ground truth image) and the three images from the low photon counts data are 0.959 (MLE), 0.983 (PMLEQ) and 0.977 (PMLETV), respectively. As shown in the results, the estimated lifetime images using the PMLE methods are more similar to the lifetime image from high photon counts data.

Fig. 5 Estimated lifetime and amplitude images from the polymer: (a) MLE from high photon counts, (b) MLE from low photon count, (c) PMLEQ from low photon counts, and (d) PMLETV from low photon counts.

Download Full Size | PDF

3.2.2. cell data

We also attempted to estimate the lifetime and amplitude images from a HeLa cell. The HeLa cell was stained with a fluorescein dye that was covalently attached to amyloid₁₁_–₂₅ (Fl-amyloid₁₁_–₂₅). The fluorescence lifetime of fluorescein is known to be pH-sensitive, and we know that the amyloid₁₁_–₂₅ is localized to specific positions in the cells. By incubating the HeLa cells with Fl-amyloid₁₁_–₂₅, it is thought that information about the probe distributions in the cells, which is related to the pH conditions of the localized states, will be obtained from FLIM images. The fluorescence decay curves were acquired using the same system as in the polymer experiment except for that a 467-nm-wavelength laser was used for this case. The average number of photon counts from pixels where the maximum number of photon count was greater than a given threshold (for this experiment 10) was 308. We applied the MLE, PMLEQ, and PMLETV methods to estimate the lifetime and amplitude maps.

Figure 6 shows the two estimated lifetimes (τ₁ and τ₂) and associated amplitudes (A₁ and A₂) obtained using the MLE method. The MLE method yielded very noisy τ₁ and τ₂ maps, thereby making it difficult to identify meaningful pH states. Figure 7 shows the estimated lifetime and amplitudes obtained using the PMLEQ method, and those obtained using the PMLETV method are shown in Fig. 8. We can see that the effects of noise have been reduced in the estimated images obtained using these two methods. Compared with the images obtained with the PMLEQ method, the lifetime maps estimated using the PMLETV method show more distinct features. Figure 8(a) and Fig. 8(b) show areas where the lifetimes are larger than the surrounding area, and we suspect that the pH states in this area are different based on the lifetime images. The variation of the fluorescence lifetime of the fluorocein as a function of pH is well known phenomenon [27, 28]. It is associated with the proton transfer of the dye in the electronically excited state. The existence of two lifetimes stems from a chemical origin in which the normal and proton-transferred species are in dynamic equilibrium. Although it is not possible to directly show that the estimated lifetime maps obtained using the PMLE methods are more similar to true lifetime maps since we do not know the true lifetimes in the HeLa cell, it is evident that the estimated images using the proposed method are less noisy than those obtained using the MLE method. Moreover, we believe that the PMLE methods provide more accurate results based on the simulated results and the assumption that the true pH states do not change abruptly in the spatial domain.

Fig. 6 Estimated lifetime and amplitude images from the HeLa cell using the MLE: (a) τ₁, (b) τ₂, (c) A₁, (d) A₂.

Download Full Size | PDF

Fig. 7 Estimated lifetime and amplitude images from the HeLa cell using the PMLEQ: (a) τ₁, (b) τ₂, (c) A₁, and (d) A₂.

Download Full Size | PDF

Fig. 8 Estimated lifetime and amplitude images from the HeLa cell using the PMLETV: (a) τ₁, (b) τ₂, (c) A₁, and (d) A₂.

Download Full Size | PDF

4. Conclusion

We have proposed a PMLE method for estimating the lifetime and amplitude images from measured fluorescence signals. Because the objective function in the proposed method considered the spatial correlation in the lifetime and amplitude images, the method was able to reduce the effects of noise more efficiently than the conventional MLE method. The objective function of the PMLE method, however, is not a separable function, and so we applied the optimization transfer principle to design a separable surrogate function to update the lifetimes and amplitudes at each pixel independently. The simulation results showed that the proposed PMLE method outperforms the conventional MLE method in terms of accuracy, and the proposed method yielded less noisy image in real experiments. We believe that the proposed method will be useful in estimating lifetime and amplitude images for FLIM, especially for cases in which the true lifetimes and amplitudes do not change rapidly in the spatial domain.

Acknowledgments

This research was supported by the Basic Science Research Program through the National Research Foundation of Korea (NRF) funded by the Ministry of Education, Science and Technology ( NRF-2012R1A1A2009370) and Ewha Global Top 5 Grant 2011 of Ewha Womans University.

References and links

1. J. R. Lakowicz, Principles of Fluorescence Spectroscopy (Kluwer Academic/Plenum, 1999). [CrossRef]

2. C. W. Chang and M.-A. Mycek, “Enhancing precision in time-domain fluorescence lifetime imaging,” J. Biomed. Opt. 15, 056013 (2010). [CrossRef] [PubMed]

3. W. Becker, “Fluorescence lifetime imaging techniques and applications,” J. Microsc. 247, 119–136 (2012). [CrossRef] [PubMed]

4. M. Kneen, J. Farinas, Y. Li, and A. Verkman, “Green fluorescent protein as a noninvasive intracellular ph indicator,” Biophys. J. 74, 1591–1599 (1998). [CrossRef] [PubMed]

5. P. J. Verveer, A. Squire, and P. I. Bastiaens, “Global analysis of fluorescence lifetime imaging microscopy data,” Biophys. J. 78, 2127–2137 (2000). [CrossRef] [PubMed]

6. Z. Bajzer, T. M. Therneau, J. C. Sharp, and F. G. Prendergast, “Maximum likelihood method for the analysis of time-resolved fluorescence decay curves,” Eur. Biophys. J. 20, 247–262 (1991). [CrossRef]

7. M. Kollner and J. Wolfrum, “How many photons are necessary for fluorescence-lifetime measurements?” Chem. Phys. Lett. 200, 199 –204 (1992). [CrossRef]

8. J. Kim and J. Seok, “Statistical properties of amplitude and decay parameter estimators for fluorescence lifetime imaging,” Opt. Express 21, 6061–6075 (2013). [CrossRef] [PubMed]

9. H. Cramer, Mathematical Methods of Statistics (PMS-9), Princeton Landmarks in Mathematics and Physics (Princeton University Press, 1999).

10. H. C. Gerritsen, M. A. H. Asselbergs, A. V. Agronskaia, and W. G. J. H. M. Van Sark, “Fluorescence lifetime imaging in scanning microscopes: acquisition speed, photon economy and lifetime resolution,” J. Microsc. 206, 218–224 (2002). [CrossRef] [PubMed]

11. L. P. Watkins and H. Yang, “Information bounds and optimal analysis of dynamic single molecule measurements,” Biophys. J. 86, 4015 – 4029 (2004). [CrossRef] [PubMed]

12. J. Philip and K. Carlsson, “Theoretical investigation of the signal-to-noise ratio in fluorescence lifetime imaging,” J. Opt. Soc. Am. A 20, 368–379 (2003). [CrossRef]

13. H. E. Grecco, P. Roda-Navarro, and P. J. Verveer, “Global analysis of time correlated single photon counting fret-flim data,” Opt. Express 17, 6493–6508 (2009). [CrossRef] [PubMed]

14. S. Pelet, M. J. R. Previte, L. H. Laiho, and P. T. C. So, “A fast global fitting algorithm for fluorescence lifetime imaging microscopy based on image segmentation.” Biophys. J. 87, 2807–2817 (2004). [CrossRef] [PubMed]

15. K. M. Hanson, M. J. Behne, N. P. Barry, T. M. Mauro, E. Gratton, and R. M. Clegg, “Two-photon fluorescence lifetime imaging of the skin stratum corneum ph gradient,” Biophys. J. 83, 1682–1690 (2002). [CrossRef] [PubMed]

16. A. Squire and P. I. H. Bastiaens, “Three dimensional image restoration in fluorescence lifetime imaging microscopy,” J. Microsc. 193, 36–49 (1999). [CrossRef] [PubMed]

17. D. Sud and M.-A. Mycek, “Image restoration for fluorescence lifetime imaging microscopy (FLIM),” Opt. Express 16, 19192–19200 (2008). [CrossRef]

18. M. Heilemann, D. P. Herten, R. Heintzmann, C. Cremer, C. Mller, P. Tinnefeld, K. D. Weston, J. Wolfrum, and M. Sauer, “High-resolution colocalization of single dye molecules by fluorescence lifetime imaging microscopy,” Anal. Chem. 74, 3511–3517 (2002). [CrossRef] [PubMed]

19. 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). [CrossRef] [PubMed]

20. 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, 381–390 (2003). [CrossRef] [PubMed]

21. J. Fessler and A. Hero, “Penalized maximum-likelihood image reconstruction using space-alternating generalized em algorithms,” IEEE Trans. Image Process. , 4, 1417–1429 (1995). [CrossRef] [PubMed]

22. J.-H. Chang, J. Anderson, and J. Votaw, “Regularized image reconstruction algorithms for positron emission tomography,” IEEE Trans. Med. Imag. , 23, 1165 – 1175 (2004). [CrossRef]

23. J. Fessler, “Image reconstruction: Algorithms and analysis,” Online preprint of book in preparation.

24. A. De Pierro, “A modified expectation maximization algorithm for penalized likelihood estimation in emission tomography,” IEEE Trans. Med. Imag. , 14, 132–137 (1995). [CrossRef]

25. P. Huber, Robust Statistics (Wiley, 1974).

26. J. Fessler and W. Rogers, “Spatial resolution properties of penalized-likelihood image reconstruction: space-invariant tomographs,” IEEE Trans. Image Process. , 5, 1346–1358 (1996). [CrossRef] [PubMed]

27. H.-J. Lin, P. Herman, and J. R. Lakowicz, “Fluorescence lifetime-resolved ph imaging of living cells,” Cytometry Part A 52A, 77–89 (2003). [CrossRef]

28. C. Hille, M. Berg, L. Bressel, D. Munzke, P. Primus, H.-G. Lahmannsraben, and C. Dosche, “Time-domain fluorescence lifetime imaging for intracellular ph sensing in living tissues,” Anal. Bioanal. Chem. 391, 1871–1879 (2008). [CrossRef] [PubMed]

		τ₁		τ₂		A₁		A₂

		Bias	STD	Bias	STD	Bias	STD	Bias	STD
MLE	ellipse	0.01	0.20	0.05	0.34	0.70	3.69	−0.09	2.16
	square	0.02	0.39	0.20	0.85	0.74	3.98	−0.43	4.28
	circle	0.01	0.24	0.04	0.30	0.70	4.12	−0.21	3.38

MLE (binnig)	ellipse	0.02	0.11	0.08	0.27	−0.57	2.63	−0.28	1.45
	square	−0.08	0.18	−0.09	0.22	−0.66	1.60	−0.91	1.47
	circle	0.05	0.12	0.06	0.16	−0.19	2.33	−2.33	1.56

PMLEQ	ellipse	0.03	0.07	0.11	0.15	0.93	0.84	−0.19	0.65
	square	−0.27	0.12	−0.36	0.22	−0.19	0.65	0.45	0.65
	circle	0.01	0.05	−0.01	0.10	−1.67	1.28	−1.91	1.24

PMLETV	ellipse	−0.01	0.04	0.03	0.07	0.68	0.44	−0.17	0.29
	square	−0.07	0.04	−0.18	0.08	−1.34	0.53	1.23	0.36
	circle	−0.01	0.04	−0.06	0.06	−1.16	0.72	0.46	0.33

	τ₁	τ₂	A₁	A₂
MLE	0.799	0.853	0.658	0.780
MLE (binning)	0.921	0.916	0.802	0.875
PMLEQ	0.946	0.952	0.941	0.960
PMLETV	0.991	0.993	0.978	0.994

		τ₁		τ₂		A₁		A₂

		Bias	STD	Bias	STD	Bias	STD	Bias	STD
MLE	ellipse	0.01	0.20	0.05	0.34	0.70	3.69	−0.09	2.16
	square	0.02	0.39	0.20	0.85	0.74	3.98	−0.43	4.28
	circle	0.01	0.24	0.04	0.30	0.70	4.12	−0.21	3.38

MLE (binnig)	ellipse	0.02	0.11	0.08	0.27	−0.57	2.63	−0.28	1.45
	square	−0.08	0.18	−0.09	0.22	−0.66	1.60	−0.91	1.47
	circle	0.05	0.12	0.06	0.16	−0.19	2.33	−2.33	1.56

PMLEQ	ellipse	0.03	0.07	0.11	0.15	0.93	0.84	−0.19	0.65
	square	−0.27	0.12	−0.36	0.22	−0.19	0.65	0.45	0.65
	circle	0.01	0.05	−0.01	0.10	−1.67	1.28	−1.91	1.24

PMLETV	ellipse	−0.01	0.04	0.03	0.07	0.68	0.44	−0.17	0.29
	square	−0.07	0.04	−0.18	0.08	−1.34	0.53	1.23	0.36
	circle	−0.01	0.04	−0.06	0.06	−1.16	0.72	0.46	0.33

	τ₁	τ₂	A₁	A₂
MLE	0.799	0.853	0.658	0.780
MLE (binning)	0.921	0.916	0.802	0.875
PMLEQ	0.946	0.952	0.941	0.960
PMLETV	0.991	0.993	0.978	0.994

Penalized maximum likelihood estimation of lifetime and amplitude images from multi-exponentially decaying fluorescence signals

Abstract

1. Introduction

2. Theory

2.1. Problem formulation

2.2. Proposed method

3. Results

3.1. Simulations

3.2. Real data

3.2.1. polymer data

3.2.2. cell data

4. Conclusion

Acknowledgments

References and links

Cited By

Figures (8)

Tables (4)

Equations (24)

Optics Express