Blind unmixing of spectrally resolved lifetime images

Farzad Fereidouni; Gerhard A. Blab; Hans C. Gerritsen

doi:10.1117/1.JBO.18.8.086006

12 August 2013 Blind unmixing of spectrally resolved lifetime images

Farzad Fereidouni, Gerhard A. Blab, Hans C. Gerritsen

Author Affiliations +

Journal of Biomedical Optics, Vol. 18, Issue 8, 086006 (August 2013). https://doi.org/10.1117/1.JBO.18.8.086006

Abstract

A method, is presented for blind unmixing spectrally resolved fluorescence lifetime images. The method is based on the combined analysis of spectral and lifetime phasors and allows unmixing of up to three components without any prior knowledge. Fractional intensities, spectra and decay curves of the individual components can be extracted with this new technique. The reliability and sensitivity are investigated and the possibility of extending the method to unmix more components is discussed. The method is evaluated on mixtures of fluorescent dyes and labeled cells.

1. Introduction

Unmixing of spectral and lifetime images has been extensively used to study the distribution and colocalization of fluorophores in specimens. In spectral or lifetime unmixing, the total spectrum or decay is decomposed into the separate components. Applications of unmixing can be found in biology, chemistry, and earth sciences and numerous unmixing algorithms have been developed. Among the techniques developed for finding the individual components, linear unmixing serves as a standard and straightforward technique. It assumes a linear superposition of different components¹ but this approach requires knowledge of reliable references for the individual components. However, the recording of reference spectra or decays is not always trivial; they may be subject to changes induced by, e.g., solvents. This complication especially occurs when one is dealing with endogenous fluorophores, in which case recording of the reference spectra or decays is often not feasible at all.

Recently, novel techniques including spectral phasor analyses (SPA),² non-negative matrix factorization (NMF),³ and principal component analysis (PCA)⁴ have been developed that overcome this problem. PCA is based on setting up a perpendicular coordinate system, called principal components, and maximization of the variances of the data points along each of the perpendicular coordinates. The method is mathematically complicated and has been successfully used, for instance, for unmixing of Flavin autofluorescence components.⁵ NMF uses constraints including non-negative spectra and intensity fractions and the minimization of a cost function to unmix the spectral images. The application of NMF for analyzing spectral lifetime images has also been reported.⁶ In this technique, a rough estimation of spectra of constituent components is still necessary as an input and the final result is dependent on the initial input. In SPA, unmixing is achieved by employing reference phasors and exploiting the linear behavior of the phasors. The reference phasors can be estimated from the recorded spectral image itself. This is performed by extracting the maximum emission wavelength of pure spectra from the average spectrum (of a selected region). The choice of the reference points is subject to errors and the accuracy of the unmixing depends on the estimate of the reference points by the user.

Efforts have been made to employ extra information, like excitation spectra or lifetimes, to alleviate the requirement of reference data and to increase the robustness of the unmixing. The parallel factor analysis⁷ was introduced to blindly unmix up to five components without prior knowledge; here, a filter wheel was employed to record fluorescence excitation images and increase the amount of data available for unmixing. This method suffers from comparatively long acquisition times and a fairly complicated acquisition process.

Schlachter et al.⁸ employed a series of lifetime images recorded at different excitation wavelengths to enhance the unmixing process. They successfully demonstrated the unmixing of two components without prior knowledge about the components. This method is, however, limited to dyes exhibiting single exponential decays, and the unmixing is performed based on a global lifetime estimation using a phasor approach.

A combined multimodal analysis has the potential to more reliably unmix the fluorophores than with a single technique. Moreover, the simultaneous recording of different properties of the fluorophore with a single detector, for instance, emission wavelength and fluorescence lifetime, produces additional information from the same number of photons compared to the recording of only a single property. This potentially opens the door to reliable results with fewer numbers of detected photons. This requires, however, that the fluorophores to be unmixed have differences in multiples of the imaging modalities. A good example is the presence of a free and bound state of NADH. The emission maxima of free and bound NADH differ by 20 nm⁹ and their lifetimes are 0.3 to 0.5 ns and 1.6 to 3.7 ns¹⁰ for free and bound NADH, respectively.

In frequency-domain lifetime imaging systems, two different modulation frequencies are needed to extract the lifetimes of a two-component mixture; this complicates the image acquisition and results in longer acquisition times. Hanley et al.¹¹ used frequency-domain spectral lifetime imaging and employed a phasor approach to unmix two components and their lifetimes using only a single-modulation frequency for the excitation light. The spectral information removes the need for the extra recording at other harmonics. In a more recent publication,¹² it is shown that the unmixing process becomes more accurate when the spectral and lifetime information are combined. It is possible to separate the signals from Rhodamine and fluorescein even when the fractional intensity of the Rhodamine is below 2%.

In this article, we report a new technique which combines spectral and lifetime data to unmix the fractional intensities, emission spectra and decays of pure components. This method relies on the use of spectral phasors and lifetime phasors and does not require any reference data for unmixing. Every spectrum or decay curve is transformed into a single point on the phasor plot; this reduces the unmixing procedure to find the reference phasor points of pure components. The method is tested on digitally added images of single dye solutions to generate a two- and a three-component system. We show that the unmixing can be performed with average photon counts as low as $400 counts / pixel$ . Moreover, images with spectra separated by only 15 nm and lifetimes separated by 700 ps can be unmixed. Finally, the method is also tested on images of labeled cells.

1.1.

Temporal and Spectral Phasor Analysis

Our novel unmixing method is based on the use of phasors. Previously, phasor analysis has been successfully employed for independently analyzing spectral and lifetime images.²^,¹³ In this approach, the decay curve or the spectrum are Fourier transformed and mapped as a single point in a phasor plot. Here, the $x$ axis is used to map the real part and the $y$ axis is used to map the imaginary part of the first harmonic of the Fourier transformation. This approach facilitates the analysis of lifetime and spectral images by visualization of the data in a two-dimensional (2-D) plot. The position of the phasor depends on the shape of the spectrum or the decay curve. There is a one to one correlation between phasor points and pixels with a particular spectrum or lifetime; this allows for simple and straightforward image segmentation. Figure 1(a) shows reference phasor curves for a time-gated system with different numbers of gates ( $K$ ) and the same total measurement window ( $T$ ); the lifetime range runs from $τ / T = 0.01$ to $τ / T = 5$ . The phasors of artificial Gaussian-shaped spectra with variable peak position $k_{0}$ and variable width $w$ are shown in Fig. 1(b). Reference curves are shown for different numbers of spectral detection channels, $K = 8$ and $K = 128$ . The phasor curves for zero width ( $w = 0$ , corresponding to a delta function) show the maximum extent of the phasors. Nonzero width spectra fall inside these curves.

Fig. 1

(a) Time-gated phasors for different (normalized) lifetimes and numbers of time gates ( $K$ ); (b) spectral phasor reference curves for different numbers of spectral channels and spectral widths.

Importantly, phasors show a linear behavior. The phasor of a linear superposition of two spectra or decay curves is a linear superposition of the phasors of the pure components. This allows direct unmixing; the relative position of the phasor of a mixture with respect to the pure phasors yields the fractional intensities of each component. In a binary mixture, the phasors fall on a line connecting the phasors of the pure components.

The relative distance of the total phasor to the reference points is determined by the fraction of that component. This is shown for both temporal and spectral phasors in Fig. 2(a) and 2(b). In the case of a three-component system, shown in Fig. 2(c) and 2(d), the total phasor falls inside a triangle whose vertices are made up by the phasors of the pure components. The fractional intensities can be obtained by calculating the area spanned by the total phasor and two vertices and dividing this by the total area of the triangle formed by the phasors of the pure components. Both spectral and temporal phasors have been used for unmixing. In temporal phasors, unmixing of two components is realized simply by fitting a line through all phasor points; the intersections of the line with the reference semicircle yield the reference phasor points (lifetimes) of the pure components. This, however, only works in the case of single exponential decays. Determining the reference points of multiexponential decays is not possible in this way and the reference points need to be assigned manually.

Fig. 2

The linear superposition of two components using (a) spectral phasor and (b) temporal phasor. For a three-component system, the linear superposition falls inside the triangle made up by the phasors of the pure components (c and d).

In the case of spectral phasors, the reference points can be estimated from the spectral image. First, the emission maxima of the pure components are estimated. This is accomplished by selecting regions in the phasor plot and showing the average spectra of the pixels that correspond to these regions. The average spectra show the emission maxima that can serve as reference values. Next, lines can be drawn in the phasor plot for the emission maxima and a range of spectral widths. By estimating straight lines through the phasor cloud, the position of suitable reference phasor points can be estimated. This is subject to errors when there are no pixels in the image with pure components and no unique triangle in the phasor cloud can be determined. In this article, we will introduce a new technique for finding the reference points in the phasor diagrams which is based on combining the information from lifetimes and spectra. The correlation between temporal and spectral phasors for a system with two and three components will be derived. This is the base for extracting all reference points, estimates of spectra and lifetimes and performs blind unmixing.

2. Methods

2.1.

Theory

Unmixing methods need a model to describe how different components are superimposed at the pixel level to estimate the fractional intensities and the pure spectra. In the absence of (interaction) effects such as solvatochromic shifts, linear superposition can be assumed. The normalized intensity $s$ of a pixel with position index $j$ in wavelength channel $n$ and time bin $m$ is related to the “fraction” $α$ of fluorophore number $k$ by:

Eq. (1)

s_{n m j} = \sum_{k} a_{n k} d_{m k} α_{k j},

where

a_{n k}

and

d_{m k}

are the fractional contribution to the intensity of spectrum

k

at wavelength channel

n

and the decay

k

at time bin

m

, respectively. The decays do not need to be single exponential and can take on any shape. Here, we assume the integrals of the spectra (

a_{n 1}, a_{n 2}, \dots

) and decay (

d_{m 1}, d_{m 2}, \dots

) to be normalized to one. Therefore, we can write:

Eq. (2)

s_{m j} = \sum_{n = 0}^{N - 1} \sum_{k} a_{n k} d_{m k} α_{k j} = \sum_{k} d_{m k} α_{k j},

Eq. (3)

s_{n j} = \sum_{m = 0}^{M - 1} \sum_{k} a_{n k} d_{m k} α_{k j} = \sum_{k} a_{n k} α_{k j},

where

N

and

M

are the number of spectral channels and time gates, respectively. The Fourier transformation of Eqs. (1)–(3) can be written as follows:

Eq. (4)

S_{j}^{A D} = \sum_{m = 0}^{M - 1} \sum_{n = 0}^{N - 1} s_{n m j} e^{i \frac{2 π}{N} q (n + \frac{1}{2})} e^{i \frac{2 π}{M} q (m + \frac{1}{2})} = \sum_{k} A_{k} D_{k} α_{k j},

Eq. (5)

S_{j}^{A} = \sum_{n = 0}^{N - 1} s_{n j} e^{i \frac{2 π}{N} q (n + \frac{1}{2})} = \sum_{k} A_{k} α_{k j},

Eq. (6)

S_{j}^{D} = \sum_{m = 0}^{M - 1} s_{m j} e^{i \frac{2 π}{M} q (m + \frac{1}{2})} = \sum_{k} D_{k} α_{k j},

where

A_{k} = \sum_{n = 1}^{N} a_{n k} e^{i (2 π / N) q [n + (1 / 2)]}

and

D_{k} = \sum_{m = 1}^{M} d_{m k} e^{i (2 π / M) q [m + (1 / 2)]}

are the discrete Fourier transforms of

a_{k}

and

d_{k}

, and

q

is the harmonic number. The Fourier transformation of Eq. (1) is basically a two dimensional transformation where one dimension is time and the other is wavelength.

For a case with two fluorophores, Eqs. (4)–(6) are simplified to:

Eq. (7)

S_{j}^{A D} = A_{1} D_{1} α_{1 j} + A_{2} D_{2} (1 - α_{1 j}),

Eq. (8)

S_{j}^{A} = A_{1} α_{1 j} + A_{2} (1 - α_{1 j}),

Eq. (9)

S_{j}^{D} = D_{1} α_{1 j} + D_{2} (1 - α_{1 j}) .

We note that normalization requires that $α_{1} + α_{2} = 1$ . Therefore, $α_{2}$ can be replaced by $1 - α_{1}$ . For all pixels (different values of $j$ ), the phasors fall onto lines given by:

Eq. (10)

S_{j}^{A D} = u^{A D} + v^{A D} S_{j}^{A},

Eq. (11)

S_{j}^{D} = u^{D} + v^{D} S_{j}^{A},

where

Eq. (12)

u^{A D} = \frac{A_{1} D_{1} - A_{2} D_{2}}{A_{1} - A_{2}}, v^{A D} = \frac{A_{1} A_{2} (D_{2} - D_{1})}{A_{1} - A_{2}}

Eq. (13)

u^{D} = \frac{D_{1} - D_{2}}{A_{1} - A_{2}}, v^{D} = \frac{A_{1} D_{2} - A_{2} D_{1}}{A_{1} - A_{2}} .

These equations are valid for both the real and imaginary parts of the Fourier transformations.

Solving $A_{1}$ , $A_{2}$ , $D_{1}$ , $D_{2}$ results in:

Eq. (14)

A_{1, 2} = \frac{v^{A D} - u^{D} \pm \sqrt{Δ}}{2 v}, D_{1, 2} = \frac{v^{A D} + u^{D} \pm \sqrt{Δ}}{2},

where

Eq. (15)

Δ = \sqrt{{(u^{D})}^{2} - 2 u^{D} v^{D} + {(v^{D})}^{2} + 4 v^{D} u^{A D}} .

The experimental analysis is performed by linear fitting of $S_{j}^{A D}$ versus $S_{j}^{A}$ and $S_{j}^{D}$ versus $S_{j}^{A}$ to estimate $u^{A D}$ , $v^{A D}$ and $u^{D}$ , $v^{D}$ , respectively and solving Eqs. (14) and (15). After finding $A$ and $D$ , unmixing proceeds as follows:

Eq. (16)

α_{1} = \frac{S_{j}^{D} - D_{2}}{D_{2} - D_{1}} or α_{1} = \frac{S_{j}^{A} - A_{2}}{A_{2} - A_{1}} .

This corresponds to the fractional distance of the phasor of the mixture to the pure phasor, see Fig. 1. Interestingly, the spectra and decays of the pure components can also be calculated.

This is achieved by generating the spectrally resolved temporal phasor and time resolved spectral phasors. This can be performed by:

Eq. (17)

S_{m}^{A} = \sum_{n = 0}^{N - 1} s_{n m} e^{i \frac{2 π}{N} (n + \frac{1}{2})} / s_{m} = \sum_{k} A_{k} d_{m k} α_{k} / \sum_{k} d_{m k} α_{k},

Eq. (18)

S_{n}^{D} = \sum_{m = 0}^{M - 1} s_{n m} e^{i \frac{2 π}{M} (m + \frac{1}{2})} / s_{n} = \sum_{k} D_{k} a_{n k} α_{k} / \sum_{k} a_{n k} α_{k} .

$S_{m}^{A}$ represents the time resolved spectral phasor; this is the Fourier transformation of the spectrum at time bin $m$ . Similarly $S_{n}^{D}$ represents the spectrally resolved temporal phasor which is the Fourier transformation of the decay at spectral channel $n$ .

For a binary system, the decays and spectra of the pure components can now be written as:

Eq. (19)

d_{m 1, 2} = \frac{S_{m}^{A} - A_{2, 1}}{A_{2, 1} - A_{1, 2}} s_{m}

and

Eq. (20)

a_{n 1, 2} = \frac{S_{n}^{D} - D_{2, 1}}{D_{2, 1} - D_{1, 2}} s_{n} .

The pure decays ( $d_{m}$ ) and spectra ( $a_{n}$ ) can be found by simply multiplying the total decay $s_{m}$ or spectrum $s_{n}$ with the appropriate fractional intensity at each time bin [Eq. (19)] or spectral channel [Eq. (20)], respectively.

In the case of a three-component system, Eqs. (4)–(6) can be rewritten as:

Eq. (21)

S_{j}^{A D} = A_{1} D_{1} α_{1 j} + A_{2} D_{2} α_{2 j} + A_{3} D_{3} (1 - α_{1 j} - α_{2 j}),

Eq. (22)

S_{j}^{A} = A_{1} α_{1 j} + A_{2} α_{2 j} + A_{3} (1 - α_{1 j} - α_{2 j}),

Eq. (23)

S_{j}^{D} = D_{1} α_{1 j} + D_{2} α_{2 j} + D_{3} (1 - α_{1 j} - α_{2 j}) .

Normalization allows replacing $α_{3}$ by $1 - α_{1} - α_{2}$ . All pixels with different fractional intensities lie in four planes defined by:

Eq. (24)

S_{j x, y}^{D} = v_{x, y}^{D} S_{j x}^{A} + w_{x, y}^{D} S_{j y}^{A} + u_{x, y}^{D}

and

Eq. (25)

S_{j x, y}^{A D} = v_{x, y}^{A D} S_{j x}^{A} + w_{x, y}^{A D} S_{j y}^{A} + u_{x, y}^{A D} .

Here, the $x$ and $y$ indices correspond to the real and imaginary parts of the Fourier transformations. We define:

Eq. (26)

M (X, Y, Z) = X_{1} Y_{2} Z_{3} - X_{1} Y_{3} Z_{2} - X_{2} Y_{1} Z_{3} + X_{2} Y_{3} Z_{1} + X_{3} Y_{1} Z_{2} - X_{3} Y_{2} Z_{1} .

Then $v$ , $w$ , and $u$ can be expressed as:

Eq. (27)

v_{x, y}^{D} = - \frac{M (T_{x, y}, A_{y}, 1)}{M (A_{x}, A_{y}, 1)} w_{x, y}^{D} = - \frac{M (T_{x, y}, A_{x}, 1)}{M (A_{x}, A_{y}, 1)} u_{x, y}^{D} = \frac{M (T_{x, y}, A_{x}, A_{y})}{M (A_{x}, A_{y}, 1)}

and

Eq. (28)

v_{x, y}^{A D} = - \frac{M (T_{x, y} A_{x, y}, A_{y}, 1)}{M (A_{x}, A_{y}, 1)} w_{x, y}^{A D} = - \frac{M (T_{x, y} A_{x, y}, A_{x}, 1)}{M (A_{x}, A_{y}, 1)} u_{x, y}^{A D} = \frac{M (T_{x, y} A_{x, y}, A_{x}, A_{y})}{M (A_{x}, A_{y}, 1)} .

The unmixing of three components relies on minimization of $| v_{x y}^{D} S_{j x}^{D} + w_{x y}^{D} S_{j y}^{D} + u_{x y}^{D} - S_{j x, y}^{D} | / \sqrt{{(v_{x y}^{D})}^{2} + {(w_{x y}^{D})}^{2} + 1}$ and $| v_{x y}^{A D} S_{j x}^{D} + w_{x y}^{A D} S_{j y}^{A D} + u_{x y}^{A D} - S_{j x, y}^{A D} | / \sqrt{{(v_{x y}^{A D})}^{2} + {(w_{x y}^{A D})}^{2} + 1}$ and yields $v$ , $w$ , and $u$ . This procedure corresponds to fitting a plane through the experimental points. The reference points of the spectral and temporal phasors are related by:

Eq. (29)

D_{k x} = v_{x}^{D} A_{k x} + w_{x}^{D} A_{k y} + u_{x}^{D} D_{k y} = v_{y}^{D} A_{k x} + w_{y}^{D} A_{k y} + u_{y}^{D}

and

Eq. (30)

A_{k x} D_{k x} = v_{x}^{A D} A_{k x} + w_{x}^{A D} A_{k y} + u_{x}^{A D} A_{k y} D_{k y} = v_{y}^{A D} A_{k x} + w_{y}^{A D} A_{k y} + u_{y}^{A D} .

Substituting the $D_{k x}$ and $D_{k y}$ from Eq. (29) into Eq. (30) yields:

Eq. (31)

v_{x}^{D} A_{k x}^{2} + (u_{x}^{D} - v_{x}^{A D} + A_{k y} w_{x}^{D}) A_{k x} - (u_{x}^{A D} A_{k y} w_{x}^{A D}) = 0 w_{y}^{D} A_{k y}^{2} + (u_{y}^{D} - w_{y}^{A D} + A_{k x} v_{x}^{D}) A_{k y} - (u_{y}^{A D} A_{k x} v_{y}^{A D}) = 0,

which can be solved numerically once

v

,

w

, and

u

are known. The fractional intensities can be found from the spectral phasor or temporal phasor. These can be calculated from the areas in the phasor diagram, see Fig. 2. The area (

T

) of a triangle with vertices

C 1

,

C 2

, and

C 3

can be calculated by:

Eq. (32)

T (C_{1}, C_{2}, C_{3}) = \frac{1}{2} (C_{1 x} C_{2 y} - C_{1 y} C_{2 x} - C_{1 x} C_{3 y} + C_{1 y} C_{3 x} + C_{2 x} C_{3 y} - C_{2 y} C_{3 x}) .

Using this, the fractional intensities as derived from the spectral phasor plot are as follows:

Eq. (33)

α_{1 j} = \frac{T (A_{2}, A_{3}, S_{j}^{A})}{T (A_{1}, A_{2}, A_{3})} α_{2 j} = \frac{T (A_{1}, A_{3}, S_{j}^{A})}{T (A_{1}, A_{2}, A_{3})} .

And from the temporal phasor:

Eq. (34)

α_{1 j} = \frac{T (D_{2}, D_{3}, S_{j}^{D})}{T (D_{1}, D_{2}, D_{3})} α_{2 j} = \frac{T (D_{1}, D_{3}, S_{j}^{D})}{T (D_{1}, D_{2}, D_{3})} .

The pure decays and spectra are found from:

Eq. (35)

d_{m 1} = \frac{T (A_{2}, A_{3}, S_{m}^{A})}{T (A_{1}, A_{2}, A_{3})} s_{m} d_{m 2} = \frac{T (A_{1}, A_{3}, S_{m}^{A})}{T (A_{1}, A_{2}, A_{3})} s_{m} d_{m 3} = \frac{T (A_{1}, A_{2}, S_{m}^{A})}{T (A_{1}, A_{2}, A_{3})} s_{m}

Eq. (36)

a_{n 1} = \frac{T (D_{2}, D_{3}, S_{n}^{D})}{T (D_{1}, D_{2}, D_{3})} s_{n} a_{n 2} = \frac{T (D_{3}, D_{1}, S_{n}^{D})}{T (D_{1}, D_{2}, D_{3})} s_{n} a_{n 3} = \frac{T (D_{1}, D_{2}, S_{n}^{D})}{T (D_{1}, D_{2}, D_{3})} s_{n} .

Similar to the binary system, the pure decays and spectra are found by multiplying the total decay spectrum with the appropriate fractional intensity at each time bin (or spectral channel).

2.2.

Instrumental Details

Two excitation sources were used for our experiments. A 473-nm solid-state diode laser (Becker and Hickl, Berlin, Germany, BDL-473-SMC) with a pulse repetition rate of 80 MHz was used for single-photon excitation experiments. In particular, it was used for the simultaneous excitation of two fluorophores to validate the unmixing of binary systems.

A mode-locked Titanium: Sapphire (Ti:Sa) laser (Tsunami, Spectra-Physics, Sunnyvale, CA) tuned to 780 nm was used in the two-photon excitation experiments. It was employed to simultaneously excite three fluorophores in the test samples. The laser light was scanned in the $X Y$ direction using a galvanometer mirror scanner (040EF, LSK, Stallikon, Switzerland). The laser was focused into the sample by a microscope objective and the fluorescence emission was collected by the same objective lens. The results reported here were acquired using an infinity-corrected water-immersion objective (CFI Fluor $60 \times$ W, $numerical aperature (NA) = 1.3$ , Nikon, Japan). The emission passed through a dichroic mirror (680-nm short pass or 490-nm long pass) and was filtered by a multiphoton emission filter (FF01-680/SP-25, Semrock, Rochester, NY) when using the Ti:Sa laser and a 490-nm long pass filter when using the diode laser. The fluorescence was coupled into a 900-nm multimode fiber which was connected to a prism-based spectrograph equipped with a multichannel photo multiplier tube (PMT) with 32 anodes of which only 7 anodes are used. The output of the PMT was amplified by a high bandwidth multichannel amplifier (Philips scientific 775-s-50) and the amplified signals were connected to a multispectral lifetime imaging system (Lambda–Tau) (see Ref. 14). The spectral width of each detection channel was about 20 nm and the time gates had a width of 200 ps. The 80-MHz repetition rate corresponds to a total time measurement window of 12.5 ns. All images shown here have a field of view of $60 \times 60 {μ m}^{2}$ and the image size is $160 \times 160 pixels$ . All experiments were carried out at room temperature.

To account for the instrument response function, the temporal phasor is calibrated by recording the single exponential fluorescence decay from an aqueous solution of Fluorescein (0.1 mM at pH 7).

2.3.

Analysis

The spectral images shown here are visualized using a real color visualization method. Briefly, the spectral information is used to assign an RGB value to each pixel that corresponds to the color as perceived by eye. To this end, the spectra are multiplied by the wavelength-dependent real-color RGB value and integrated over the spectral channels.¹⁵ The temporal and spectral phasor plots were generated using homemade ImageJ plugins (plugins are available from www.Spechron.com). The Fortran IMSL library was used for data fitting and root finding. The lifetimes of both unmixed and original data were estimated by means of the Weber method.¹⁶

2.4.

Sample Details

Aqueous solutions (0.1 mM) with $pH = 7$ were prepared from Fluorescein (Sigma–Aldrich F6377), Alexa 532 (Invitrogen A-20001), Alexa 555 (Invitrogen A-20009), and Alexa 647 (Invitrogen A-20006). Imaging experiments on fixed cell were carried out using FluoCells prepared slide #2 (Invitrogen F-14781) and #6 (Invitrogen F-35925).

Slide #2 contained bovine pulmonary artery endothelial cells (BPAEC). Here, Texas Red-X phalloidin was used for labeling F-actin, anti-bovine $α$ -tubulin mouse monoclonal antibody in conjunction with BODIPY FL goat anti-mouse IgG antibody was used for labeling microtubules, and 4',6-diamidino-2-phenylindole (DAPI) was used for labeling the nucleus.

Slide #6 contained fixed muntjac skin fibroblast cells. Mitochondria were labeled with mouse anti-OxPhos Complex V inhibitor protein antibody and visualized using orange-fluorescent Alexa Fluor 555 goat anti-mouse IgG antibody. F-actin was labeled with green-fluorescent Alexa Fluor 488 phalloidin, and the nucleus was stained with TO-PRO-3 iodide.

The plant leaf sample was freshly picked near the laboratory and immersed in water during the experiments to match the water immersion objective.

3. Results and Discussion

3.1.

Unmixing of Two Components

In order to validate the unmixing method, separate images of Alexa 532 and Alexa 555 solutions were recorded, which exhibited a clear intensity gradient. Next, they were added digitally to yield a well-defined binary image.

This approach is preferred over the use of synthetic images with added Poisson noise¹³ It yields images closely resembling the images of real mixtures. The images had $400 counts / pixel$ on average. Figure 3 shows the phasor plots of the mixture as well as the intensity images, decay curves, and spectra after unmixing. The decay curve and the residual of the unmixed curve for Alexa 532 are shown in Fig. 3(d). Using the Weber method for lifetime, estimation resulted in average lifetimes of 2.83 and 2.69 ns for the original and unmixed Alexa 532, respectively. The same approach resulted in lifetimes of 0.33 and 0.33 ns for the original and unmixed Alexa 555, respectively. The spectra from the original and unmixed data are in excellent agreement; this is demonstrated by the small residuals of both components, see Fig. 3(e) and 3(h). In order to quantify the quality of the unmixed data including the curves and images, the average relative error was calculated:

Eq. (37)

δ = \sum_{n} \frac{| s_{unmixed} - s_{original} |}{s_{original}} / N .

Fig. 3

Temporal (a) and spectral (b) phasor plots of the Alexa 532-Alexa 555 mixture image. The images (c, f), decays (d, g) and spectra from unmixed and original data are shown (e, h).

The summation runs over all ( $N$ ) spectral, time channels or pixel numbers. The results are summarized in Table 1. The average unmixed fraction of Alexa 532 in the digitally added image was calculated from the two original images and found to be 66.66%. This is in good agreement with the value of 66.62% found using the blind unmixing procedure.

Table 1

Summary of the unmixed and original data for the Alexa 532–Alexa 555 image.

Component	$τ_{original}$ (ns)	$τ_{unmixed}$ (ns)	$δ_{τ}$ (%)	$δ_{λ}$ (%)	$δ_{image}$ (%)
Alexa 532 ( $400 counts / pixel$ )	2.83	2.69	5.7	8.3	1.67
Alexa 555( $400 counts / pixel$ )	0.33	0.33	7.1	1.7	1.58
Alexa 532 ( $150 counts / pixel$ )	2.87	2.79	6.2	11.8	6.09
Alexa 555 ( $150 counts / pixel$ )	0.33	0.34	11.4	1.9	5.88

The temporal and spectral phasor diagrams are shown in Fig. 3(a) and 3(b), respectively; the elongated shape of the phasors is characteristic for a two-component system. The reference points for the phasors are indicated in both diagrams and show that the two components can be unmixed when they show different lifetimes by employing the temporal phasor diagram. To this end, a line was fitted through the phasor points and the intersections of the line with the reference phasor curve yields the lifetimes of the pure components [Fig. 3(a)]. In addition, the fractional intensities were estimated from the distance between the phasor point and the intersections.¹³ By applying this method to the Alexa 532 and Alexa 555 images, lifetimes of 3.03 and 0.32 ns were obtained for Alexa 532 and Alexa 555. The average fractional intensity for Alexa 532 was estimated at 58% which is different from the 66.66% found from the combined analysis. The difference is due to the fact that the fractional intensities are estimated based on the assumption that the reference points lie on the reference semicircle. The reference points calculated from the combined analysis fall on the fitted line but do not lie on the semicircle. The discrepancy between the positions of the reference points is explained by the multiexponential decay of the dyes; only monoexponential decays fall on the semicircle.

The combined analysis does not require the decays to be single exponentials; it works also with multiexponential decays. We note that the emission maxima of Alexa 532 and Alexa 555 are separated by only 15 nm and the images used for Fig. 4(a) and 4(b) contain only $400 counts / pixel$ on an average. Nevertheless, good unmixing results are obtained.

Fig. 4

The comparison between original and unmixed image intensity profile binned in $y$ direction for (a) Alexa 532 and (b) Alexa 555 with $400 counts / pixel$ and with $150 counts / pixel$ for (c) Alexa 532 and (d) Alexa 555.

Figure 4(a) and 4(b) show comparisons between two intensity traces in the original and unmixed Alexa 532 (a) and Alexa 555 (b) images at an average intensity of $400 counts / pixel$ . The pixels are binned in the $y$ direction and the graph shows the total intensity in the $x$ direction for every component separately. Very good overlap is observed for all pixels in the image. Quantitative comparison of the Alexa 532 and Alexa 555 intensities in the original and unmixed images yields average relative errors of 1.67% and 1.58%, respectively (see Table 1). This indicates good overlap between the original and unmixed images.

The unmixing process was repeated for a lower number of average counts ( $150 counts / pixel$ ); the results are shown in Fig. 4(c) and 4(d). Quantitative comparison of the intensities yields average relative errors of 6.09% and 5.88% for Alexa 532 and Alexa 555, respectively. The summary of the results are shown in Table 1.

To test the method for the unmixing of three components, images of solutions of Fluorescein, Alexa 532, and Alexa 555 were recorded and added digitally. The images have gradients from right to left for Fluorescein, from left to right for Alexa 532, and from top to bottom for Alexa 555, see Fig. 5(c), 5(f) and 5(j). The gradients were generated by shifting the excitation beam slightly off center. The reference spectra and decays are shown in Fig. 5(d) and 5(e) for Fluorescein, Fig. 5(g) and 5(i) for Alexa 532 and Fig. 5(k) and 5(l) for Alexa 555. The unmixed images for Fluorescein, Alexa 532, and Alexa 647 are shown in Fig. 5(c), 5(f) and 5(g), respectively. The average fractional intensities of Fluorescein, Alexa 532, and Alexa 555 in the original images amount to 50.96%, 29.89%, and 19.13%, respectively. Blind unmixing of the digitally added images results in average fractional intensities of 51.99%, 28.25%, and 19.75%, respectively.

Fig. 5

Temporal (a) and spectral (b) phasors for the Fluorescein–Alexa 532–Alexa 647 mixture image. Unmixed images (c, f and i), corresponding decays (d, g, j) and spectra from unmixing and original data are shown (e, h, k).

The temporal and spectral phasor clouds in Fig. 5(a) and 5(b) show a triangle with clear vertices at the Fluorescein and Alexa 532 phasors. This is due to the fact that there are clear regions in the composed image that contain pure components from each dye. The average lifetimes of the original and unmixed dyes are found to be 3.79 and 3.62 ns for Fluorescein, 2.85 and 2.71 ns for Alexa 532, and 0.33 and 0.32 ns for Alexa 555. The results for the lifetime estimations and average relative error $δ$ are summarized in Table 2.

Table 2

Comparison between the unmixed and original data for the Fluorescein–Alexa 532–Alexa 647 mixture image.

Component	$τ_{original}$ (ns)	$τ_{unmixed}$ (ns)	$δ_{τ}$ (%)	$δ_{λ}$ (%)	$δ_{image}$ (%)
Fluorescein	3.79	3.62	2.93	5.97	2.61
Alexa 532	2.85	2.71	5.79	2.51	2.80
Alexa 555	0.33	0.32	7.14	5.49	1.7

3.2.

Unmixing of Two Components in Labeled Cells

Images of FluoCells (prepared slide #6) were recorded and blind unmixed, see Fig. 6. The unmixed fractional intensity images extracted decay curves and spectra are shown for both components. Visual inspection of the unmixed images suggests good separation of the signals of the two dyes. One of the components shows an emission maximum in the second spectral channel with a central wavelength of 522 nm and has an average fluorescence lifetime of 2.23 ns. This component corresponds to Alexa 488 used for staining the actin fibers in the specimen.

Fig. 6

The RGB image (a) and unmixed images (d, g) of FluoCells prepared slide #6. The temporal (b) and spectral phasors (c), estimated decays (e, h) and spectra (f, i) are shown.

The maximum emission wavelength of the second component is found in the fourth spectral channel with a central wavelength of 586 nm and it has an average lifetime of 1.65 ns. This component corresponds to Alexa 555 used to label the mitochondria. The fraction of this component amounts to less than 3% of the total intensity in the image. This demonstrates the sensitivity of the method for resolving small fractions. We note that the average lifetimes of the two components are close and they are separated by almost 700 ps, which explains the small size of the phasor cloud in the time domain [Fig. 6(b)]. In the spectral domain, however, the components are separated by almost 60 nm which results in the elongated shape of the phasor cloud in the spectral domain [Fig. 6(c)]. The Alexa 555 overlaps with Alexa 488 everywhere in the image and for that reason, there is almost no spectral phasor point close to the phasor of pure Alexa 555. On the other hand, there are some pixels which show pure Alexa 488 and this is reflected in both the spectral and temporal phasor plots where some phasor points reside very close to the pure Alexa 488 phasor point.

3.3.

Unmixing of Three Components in Labeled Cells

The image of the FluoCells (prepared slide #2) was obtained by two-photon excitation imaging at an excitation wavelength of 780 nm. This enabled exciting all three dyes in the specimen simultaneously. The spectral range of the spectrograph was shifted to include the blue emission from DAPI. The RGB image, the unmixed components, and the extracted spectra and decays are shown in Fig. 7. Clear segmentation is achieved and the maximum emission wavelengths of DAPI, BODIPY, and Texas red are found at 457, 505, and 636 nm, respectively, in reasonable agreement with the literature values of 461, 505 and 620 nm.¹⁷ We note that the average width of the spectral channels is 20 nm at the blue end of the spectrum and 60 nm at red end, which explains the comparatively large error in the estimation of the Texas red emission maximum.

Fig. 7

RGB image (a) and unmixed images (d, g, j) of FluoCells prepared slide #2). The temporal (b) and spectral phasors (c) and estimated decays (e, h, k) and spectra (f, i, l) are also shown.

The reference temporal and spectral phasors of the pure dyes as obtained from the analysis are depicted in Fig. 7(b) and 7(c). Inspection of the reference points in the temporal phasor shows a comparatively short lifetime for DAPI and close lifetimes for BODIPY and Texas Red. BODIPY shows clearly a multiexponential behavior, its temporal phasor falls inside the reference semicircle. Moreover, examination of the phasor plots shows a distinct cloud for DAPI which indicates little overlap with the other components, but the elongation of the phasor points between Texas red and BODIPY is an indication of overlap between these two components. This is also confirmed by the unmixed images, DAPI has no overlap with BODIPY but it does overlap with Texas red. The average lifetimes of 2.24, 2.98, and 3.54 ns are estimated for DAPI, BODIPY, and Texas red, respectively. This lifetime was estimated by the Weber method. The close lifetimes explain the condensed shape of the temporal phasor distribution, whereas the spectral phasor distribution is less packed together.

We note that the fast and straight forward segmentation of images can be realized using the phasor approach. The phasor points and the image pixels are correlated; a region of interest in the phasor plot can be projected back and highlighted in the original image. This provides a means to visualize a specific population of fluorescent molecules in the image. An example for the DAPI, BODIPY, and Texas red–labeled specimen is shown in Fig. 8. Here, the regions of interests are circular regions selected around the reference points of the pure components. The segmented (back projected) images and the estimated spectra are shown. This segmentation yields less accurate results than the above real unmixing. The segmented image and spectrum of DAPI show a reasonable agreement with the unmixed image. However, the Texas red and BODIPY results are of poor quality due to the large spectral overlap between Texas red and BODIPY.

Fig. 8

RGB image (a) and segmented images (c, e, g) of FluoCells prepared slide #2 based on back projecting regions in the phasor plot to the original image. The spectral phasor (b) and region of interests and estimated spectra of DAPI, BODIPY and Texas red (d, f, h) are also shown.

4. Conclusion

A novel phasor-based method is introduced for blind unmixing of spectral lifetime images. The combined analysis of temporal and spectral phasors yields fractional intensities, spectra and decay curves for up to three fluorescent dyes without any prior knowledge about the fluorophores. The temporal and spectral phasors for a mixture of dyes are correlated and for a binary system we demonstrated a linear relationship between the spectral and temporal phasors. This can be extended to a system with three fluorophores. In this case, the phasors with different fractional intensities behave linearly and lie in a plane.

Using a phasor approach considerably simplifies the unmixing process. By converting the spectra and decay curves into phasors and plotting the phasors of all pixels in 2-D diagrams, the unmixing process is reduced to finding the reference points of the pure dyes in the phasor plot and a linear interpolation. Both the real and imaginary parts of the Fourier transforms show the linear and planar relationships between the phasor points. For unmixing, a sample with only two components, only the real or the imaginary part of the phasors, are required.

For a sample with three components, all the phasor points fall in a plane. In this case, we use both the real and imaginary parts of the spectral phasors as $x$ and $y$ coordinates and the real and imaginary parts of time and spectral-time phasors as $z$ .

We note that we only use the first harmonic of the Fourier transform because it has the highest amplitude and is therefore least sensitive to noise. A very important feature of this method is that it does not need the decay curves to be single exponential, and multiexponential decays can be analyzed and used in the unmixing as well.

Generalization of the method to a higher number of components is feasible but requires the employment of higher harmonics to relate the different Fourier transformations. Visualizing the phasors of more than three components is not trivial and requires a three-dimensional (3-D) representation to visualize the cloud of phasor points. However, compared to SPA, the reference points are not estimated manually from the phasor cloud and visualization in 3-D is therefore not necessary.

The blind unmixing method presented here is a global method; lifetimes and spectra are assumed to be invariant in the image. This enhances the estimation accuracy considerably due to employment of all photons in the image. All the pixels take part in the unmixing process, which yields superior results even in the presence of components with low fractional intensities. In one of the examples presented here, a component with a contribution of less than 3% to the image is well resolved.

The assumption of invariant lifetimes and spectra is not valid in all cases. For instance, in Förster Resonance Energy Transfer (FRET) experiments¹⁸ the lifetime of the donor dye may vary from pixel to pixel which is not compatible with our unmixing method.

We assume the fluorescence lifetime to be invariant over the whole spectral range. This is a good approximation in many cases. The fluorescence lifetime usually shows picosecond variations over the spectrum due to the effects of internal conversion.¹⁹ Such small variations, however, hardly influence the unmixing accuracy. Even small lifetime changes of Alexa 647 over the spectrum hardly affect the unmixing process. In spite of this, the fractional intensities could be estimated with high accuracy.

The reference points for a two-component system can be found analytically. This simplifies the unmixing procedure and reduces the total time required for unmixing, extracting the fractional intensity maps, decay curves and spectra. This whole procedure takes less than 5 s for a $160 \times 160 pixel$ spectral lifetime image. When more components are present, the reference points cannot be found analytically any longer and numerical methods have to be used. Using the FORTRAN IMSL library for solving, the equations numerically slightly increase the total time for unmixing to 6 s. Optimization of the code and implementation of multithreading should significantly speed up the unmixing.

Finally, the phasor-based unmixing of spectral lifetime images is a fast and robust technique. It does not require any prior knowledge and it provides information about the spectra, decay curves, and fractional intensities of the components.

References

1.

T. Zimmermann, “Spectral imaging and linear unmixing in light microscopy,” Adv Biochem Engin/Biotechnol., 95 245 –265 (2005). MRTEEO 1059-910X Google Scholar

2.

F. FereidouniA. N. BaderH. C. Gerritsen, “Spectral phasor analysis allows rapid and reliable unmixing of fluorescence microscopy spectral images,” Opt. Express, 20 (12), 12729 –12741 (2012). http://dx.doi.org/10.1364/OE.20.012729 OPEXFF 1094-4087 Google Scholar

3.

R. A. Neheret al., “Blind source separation techniques for the decomposition of multiply labeled fluorescence images,” Biophys. J., 96 (9), 3791 –3800 (2009). http://dx.doi.org/10.1016/j.bpj.2008.10.068 BIOJAU 0006-3495 Google Scholar

4.

I. Jolliffe, Principal Component Analysis, Spinger-Verlag, Spinger-Verlag (2002). Google Scholar

5.

D. Chorvat Jr.et al., “Spectral unmixing of flavin autofluorescence components in cardiac myocytes,” Biophys. J., 89 (6), L55 (2005). http://dx.doi.org/10.1529/biophysj.105.073866 BIOJAU 0006-3495 Google Scholar

6.

P. PandeB. E. ApplegateJ. A. Jo, “Application of non-negative matrix factorization to multispectral FLIM data analysis,” Biomed. Opt. Express, 3 (9), 2244 –2262 (2012). http://dx.doi.org/10.1364/BOE.3.002244 BOEICL 2156-7085 Google Scholar

7.

H. ShirakawaS. Miyazaki, “Blind spectral decomposition of single-cell fluorescence by parallel factor analysis,” Biophys. J., 86 (3), 1739 –1752 (2004). http://dx.doi.org/10.1016/S0006-3495(04)74242-3 BIOJAU 0006-3495 Google Scholar

8.

S. Schlachteret al., “A method to unmix multiple fluorophores in microscopy images with minimal a priori information,” Opt. Express, 17 (25), 22747 –22760 (2009). http://dx.doi.org/10.1364/OE.17.022747 OPEXFF 1094-4087 Google Scholar

9.

H. D. Vishwasraoet al., “Conformational dependence of intracellular NADH on metabolic state revealed by associated fluorescence anisotropy,” J. Biol. Chem., 280 (26), 25119 –25126 (2005). http://dx.doi.org/10.1074/jbc.M502475200 JBCHA3 0021-9258 Google Scholar

10.

J. A. Paleroet al., “In vivo monitoring of protein-bound and free NADH during ischemia by nonlinear spectral imaging microscopy,” Biomed. Opt. Express, 2 (5), 1030 –1039 (2011). http://dx.doi.org/10.1364/BOE.2.001030 BOEICL 2156-7085 Google Scholar

11.

Q. HanleyA. Clayton, “AB-plot assisted determination of fluorophore mixtures in a fluorescence lifetime microscope using spectra or quenchers,” J. Microsc., 218 (1), 62 –67 (2005). http://dx.doi.org/10.1111/jmi.2005.218.issue-1 JMICAR 0022-2720 Google Scholar

12.

Y. C. ChenR. Clegg, “Spectral resolution in conjunction with polar plots improves the accuracy and reliability of FLIM measurements and estimates of FRET efficiency,” J. Microsc., 244 (1), 21 –37 (2011). http://dx.doi.org/10.1111/jmi.2011.244.issue-1 JMICAR 0022-2720 Google Scholar

13.

F. Fereidouniet al., “A modified phasor approach for analyzing time-gated fluorescence lifetime images,” J. Microsc., 244 (3), 248 –258 (2011). http://dx.doi.org/10.1111/jmi.2011.244.issue-3 JMICAR 0022-2720 Google Scholar

14.

F. FereidouniK. ReitsmaH. Gerritsen, “High speed multispectral fluorescence lifetime imaging,” Opt. Express, 21 (10), 11782 –11782 (2013). http://dx.doi.org/10.1364/OE.21.011769 Google Scholar

15.

J. A. PaleroH. S. De Bruijn, “Spectrally resolved multiphoton imaging of in vivo and excised mouse skin tissues,” Biophys. J., 93 (3), 992 (2007). http://dx.doi.org/10.1529/biophysj.106.099457 BIOJAU 0006-3495 Google Scholar

16.

S. Schlachteret al., “mhFLIM: resolution of heterogeneous fluorescence decays in widefield lifetime microscopy,” Opt. Express, 17 (3), 1557 –1570 (2009). http://dx.doi.org/10.1364/OE.17.001557 OPEXFF 1094-4087 Google Scholar

17.

F. Bestvateret al., “Two-photon fluorescence absorption and emission spectra of dyes relevant for cell imaging,” J. Microsc., 208 (2), 108 –115 (2002). http://dx.doi.org/10.1046/j.1365-2818.2002.01074.x JMICAR 0022-2720 Google Scholar

18.

T. W. Gadella, FRET and FLIM Techniques, Elsevier Science, Amsterdam, Netherlands (2008). Google Scholar

19.

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

Citation Download Citation

Farzad Fereidouni, Gerhard A. Blab, and Hans C. Gerritsen "Blind unmixing of spectrally resolved lifetime images," Journal of Biomedical Optics 18(8), 086006 (12 August 2013). https://doi.org/10.1117/1.JBO.18.8.086006

Published: 12 August 2013

Access the abstract

JOURNAL ARTICLE
13 PAGES

DOWNLOAD PAPER SAVE TO MY LIBRARY

GET CITATION

CITATIONS

Cited by 13 scholarly publications.

Explore citations on Lens.org

RIGHTS & PERMISSIONS

Get copyright permission Get copyright permission on Copyright Marketplace

KEYWORDS

Image segmentation

Luminescence

Clouds

Visualization

Binary data

RGB color model

Superposition

1.

Introduction

1.1.

Temporal and Spectral Phasor Analysis

Fig. 1

Fig. 2

2.

Methods

2.1.

Theory

Eq. (1)

Eq. (2)

Eq. (3)

Eq. (4)

Eq. (5)

Eq. (6)

Eq. (7)

Eq. (8)

Eq. (9)

Eq. (10)

Eq. (11)

Eq. (12)

Eq. (13)

Eq. (14)

Eq. (15)

Eq. (16)

Eq. (17)

Eq. (18)

Eq. (19)

Eq. (20)

Eq. (21)

Eq. (22)

Eq. (23)

Eq. (24)

Eq. (25)

Eq. (26)

Eq. (27)

Eq. (28)

Eq. (29)

Eq. (30)

Eq. (31)

Eq. (32)

Eq. (33)

Eq. (34)

Eq. (35)

Eq. (36)

2.2.

Instrumental Details

2.3.

Analysis

2.4.

Sample Details

3.

Results and Discussion

3.1.

Unmixing of Two Components

Eq. (37)

Fig. 3

Table 1

Fig. 4

Fig. 5

Table 2

3.2.

Unmixing of Two Components in Labeled Cells

Fig. 6

3.3.

Unmixing of Three Components in Labeled Cells

Fig. 7

Fig. 8

4.

Conclusion

References

Show All Keywords

Keywords/Phrases

Search In:

Publication Years