OSA's Digital Library

Virtual Journal for Biomedical Optics

Virtual Journal for Biomedical Optics

| EXPLORING THE INTERFACE OF LIGHT AND BIOMEDICINE

  • Editors: Andrew Dunn and Anthony Durkin
  • Vol. 8, Iss. 9 — Oct. 2, 2013
« Show journal navigation

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

Jeongtae Kim, Jiyeong Seok, Hwiin Lee, and Minyung Lee  »View Author Affiliations


Optics Express, Vol. 21, Issue 17, pp. 20240-20253 (2013)
http://dx.doi.org/10.1364/OE.21.020240


View Full Text Article

Acrobat PDF (840 KB)





Browse Journals / Lookup Meetings

Browse by Journal and Year


   


Lookup Conference Papers

Close Browse Journals / Lookup Meetings

Article Tools

Share
Citations

Abstract

We investigated the penalized maximum likelihood estimation of lifetime and amplitude images for fluorescence lifetime imaging microscopy. The proposed method penalizes large variations in the lifetimes and amplitudes in the spatial domain to reduces noise in the images, which is a serious problem in the conventional maximum likelihood estimation method. For an effective optimization of the objective function, we applied an optimization transfer method that is based on a separable surrogate function. Simulations show that the proposed method outperforms the conventional MLE method in terms of the estimation accuracy, and the proposed method yielded less noisy images in real experiments.

© 2013 OSA

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

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

, 2

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

] because the lifetime of a fluorophore is dependent on the molecular environment [3

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

]. FLIM has been used to study, for example, the intracellular pH distribution from an estimated lifetime map [3

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

,4

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]

] and fluorescence resonance energy transfer in a cell [5

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]

]. 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 χ2 or the Poisson likelihood function [1

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

]. It was reported that maximum likelihood estimation (MLE) based on maximizing the Poisson likelihood function outperforms other methods for TCSPC-based FLIM [6

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]

8

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]

]. Alternatively, the lifetime can be estimated in the frequency domain using intensity-modulated light, the associated phase shifts and signal demodulation [1

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

]. 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

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]

], 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

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

12

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]

]. 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

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]

] but since the exposure time is usually limited due to photobleaching and imaging speed, such an attempt may not be successful [5

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]

].

To improve the performance of the lifetime estimation without increasing the exposure time, global analysis-based methods have been studied [5

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]

, 13

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]

]. A central assumption is that the lifetimes of all the fluorophores in an image are the same [5

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]

, 13

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]

]. 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

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]

, 13

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

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]

]. 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

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

, 15

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]

].

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

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]

24

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]

]. 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

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]

23

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

] 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

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]

23

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

]. 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 (xi, yj) at a time instance tk by the Poisson random variable, g(tk, xi, yj), as follows [1

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

]:
g(tk,xi,yj)=Poisson{λ(tk;θi,j)},k=0,1,.,K1,
(1)
where the mean of the Poisson random variable is defined as
λ(tk;θi,j)=h(tk)*f(tk;θi,j).
(2)
The impulse response function (IRF) of the measurement system is h(t), and f(t) is the multi-exponential decay function at (xi, yj) defined by
f(t;θi,j)=p=0P1Ap(xi,yj)etτp(xi,yj),
(3)
where τp(xi, yj) is the p-th decay constant at (xi, yj), and Ap(xi, yj) is the associated amplitude. We define the parameter vector at (xi, yj) as
θi,j=[A0(xi,yj),A1(xi,yj),..,AP1(xi,yj),τ0(xi,yj),Δτ1(xi,yj),..,ΔτP1(xi,yj)].
(4)
Using the parameter vector defined in Eq. (4), the p-th decay time constant can be represented by
τp(xi,yj)=τ0(xi,yj)+p=1P1Δτp(xi,yj),p=1,2,..,P1.
(5)
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 = {gi,j|i = 0,,N − 1, j = 1,,M − 1}, where gi,j = [g(t0, xi, yj),., g(tK−1, xi, yj)]. 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)C},
(6)
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)=(i,j)CLi,j(θi,j;gij),
(7)
where the log-likelihood function of the parameter vector θi,j from the measurement gij is
Li,j(θi,j;gij)=k=0K1λ(tk;θi,j)+g(tk,xi,yj)logλ(tk;θi,j)logg(tk,xi,yj)!.
(8)
The MLE of Θ from the measured data g is given by
Θ^=argmaxΘ0L(Θ;g).
(9)
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:
θ^i,j=argmaxθi,j>0Li,j(θi,j;gi,j),(i,j)C.
(10)
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(Θ),
(11)
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(Θ)=p=02P1(i,j)Cψ([DijhΘ]p)+ψ([DijvΘ]p),
(12)
where the change of the p-th component at the (i, j)-th pixel location in the horizontal direction is measured by
ψ([DijhΘ]p)={ϕ([θi+1,jθi,j]p),if(i+1,j)C0,otherwise,
(13)
and that in the vertical direction is measured by
ψ([DijvΘ]p)={ϕ([θi,j+1θi,j]p),if(i,j+1)C0,otherwise,
(14)
where the 2P × (2P × 𝒞) matrices Dijh and Dijv 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

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

]. We instead adopt an edge preserving penalty function such as the TV penalty function, which is defined as follows [23

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

]:
ϕ(x)=x2+ε2,
(15)
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

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

]:
Φs(θ;θk)Φ(θ),
(16)
Φs(θk;θk)=Φ(θk).
(17)
That is, the following sequence θk monotonically converges to the minimizer of the objective function Φ(θ):
θk+1=argminθΦs(θ;θk).
(18)
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

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

] to design the quadratic surrogate functions of the penalty functions and then apply De Pierre’s approach [24

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]

] to construct the quadratic surrogate functions separately as follows [21

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]

23

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

]:
ϕ([θi+1,jθi,j]p)q([DijhΘ]p;[DijhΘk]p)=q(12{2[θi+1,jθi+1,jk]p}+12{2[θi,jθi,jk]p}+[DijhΘk]p;[DijhΘk]p)12q(2[θi+1,jθi+1,jk]p+[DijhΘk]p;[DijhΘk]p)+12q(2[θi,jθi,jk]p+[DijhΘk]p;[DijhΘk]p),
(19)
where q(a;b)=ϕ(b)+ϕ˙(b)(ab)+12c(ab)2 and φ̇ is the derivative of φ [23

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

]. The curvature c that guarantees the first inequality of a quadratic or TV penalty function can be determined using Huber’s method [25

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

]. The last inequality in Eq. (19) is due to the convex inequality [24

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]

]. Therefore, for the p-th component, the following inequality holds:
(i,j)Cψ([DijhΘ]p)(i,j)CRh([θi,j]p;Θpk),
(20)
where,
Rh([θi,j]p;Θk)={12q(2[θi,jθi,jk]p+[Di1,jhΘk]p;[Di1,jhΘk]p)+12q(2[θi,jθi,jk]p+[Di,jhΘk]p;[Di,jhΘk]p),(i1,j)C,(i+1,j)C12q(2[θi,jθi,jk]p+[Di1,jhΘk]p;[Di1,jhΘk]p),(i1,j)C,(i+1,j)C12q(2[θi,jθi,jk]p+[Di,jhΘk]p;[Di,jhΘk]p),(i1,j)C,(i+1,j)C0,(i+1,j)C,(i1,j)C.
It is relatively straightforward to design a similar surrogate function Rv([θi,j]p; Θk) for ψ([DjvΘ]p)=ψ([θi,j+1θi,j]p). Based on Eq. (20), the separable quadratic surrogate function of the penalty function, Rs(Θ), satisfies the following inequality:
R(Θ)Rs(Θ;Θk)=(i,jC)Rhv(θi,j;Θk),
(21)
where
Rhv(θi,j;Θk)=p=02P1Rh([θi,j]p;Θk)+Rv([θi,j]p;Θk).
(22)
Note that Rk) = Rskk). 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,jk+1=argminθi,j>0Lij(θi,j;gi,j)+βRs(θi,j;Θk).
(23)
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 (τ1 and τ2) and amplitudes (A1 and A2) 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)

table-icon
View This Table
| View All Tables

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 (τ1, τ2, A1, A2) 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) τ1, (b) τ2, (c) A1, and (d) A2.
Fig. 2 Estimated lifetime and amplitude images obtained with PMLEQ: (a) τ1, (b) τ2, (c) A1, and (d) A2.
Fig. 3 Estimated lifetime and amplitude images obtained with PMLETV: (a) τ1, (b) τ2, (c) A1, and (d) A2.
Fig. 4 Estimated lifetime and amplitude images obtained using MLE with binning of 3×3 window: (a) τ1, (b) τ2, (c) A1, and (d) A2.

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

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]

]. 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

table-icon
View This Table
| View All Tables

We also evaluated the similarity between the estimated and true images by computing the correlation coefficients of the four estimated images (τ1, τ2, A1, A2). 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

table-icon
View This Table
| View All Tables

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

table-icon
View This Table
| View All Tables

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.

3.2.2. cell data

Figure 6 shows the two estimated lifetimes (τ1 and τ2) and associated amplitudes (A1 and A2) obtained using the MLE method. The MLE method yielded very noisy τ1 and τ2 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

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

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]

]. 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) τ1, (b) τ2, (c) A1, (d) A2.
Fig. 7 Estimated lifetime and amplitude images from the HeLa cell using the PMLEQ: (a) τ1, (b) τ2, (c) A1, and (d) A2.
Fig. 8 Estimated lifetime and amplitude images from the HeLa cell using the PMLETV: (a) τ1, (b) τ2, (c) A1, and (d) A2.

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]

OCIS Codes
(100.3190) Image processing : Inverse problems
(180.2520) Microscopy : Fluorescence microscopy
(300.6280) Spectroscopy : Spectroscopy, fluorescence and luminescence

ToC Category:
Microscopy

History
Original Manuscript: June 26, 2013
Revised Manuscript: August 12, 2013
Manuscript Accepted: August 13, 2013
Published: August 21, 2013

Virtual Issues
Vol. 8, Iss. 9 Virtual Journal for Biomedical Optics

Citation
Jeongtae Kim, Jiyeong Seok, Hwiin Lee, and Minyung Lee, "Penalized maximum likelihood estimation of lifetime and amplitude images from multi-exponentially decaying fluorescence signals," Opt. Express 21, 20240-20253 (2013)
http://www.opticsinfobase.org/vjbo/abstract.cfm?URI=oe-21-17-20240


Sort:  Author  |  Year  |  Journal  |  Reset  

References

  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. Express21, 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. A20, 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. Express17, 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. Express16, 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 A52A, 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]

Cited By

Alert me when this paper is cited

OSA is able to provide readers links to articles that cite this paper by participating in CrossRef's Cited-By Linking service. CrossRef includes content from more than 3000 publishers and societies. In addition to listing OSA journal articles that cite this paper, citing articles from other participating publishers will also be listed.


« Previous Article  |  Next Article »

OSA is a member of CrossRef.

CrossCheck Deposited