Complex wavelet filter improves FLIM phasors for photon starved imaging experiments

: Fluorescence lifetime imaging microscopy (FLIM) with phasor analysis provides easy visualization and analysis of fluorophores’ lifetimes which is valuable for multiple applications including metabolic imaging, STED imaging, FRET imaging and functional imaging. However, FLIM imaging typically suffers from low photon budgets, leading to unfavorable signal to noise ratios which in many cases prevent extraction of information from the data. Traditionally, median filters are applied in phasor analysis to tackle this problem. This unfortunately degrades high spatial frequency FLIM information in the phasor analysis. These high spatial frequency components are typically edges of features and puncta, which applies to membranes, mitochondria, granules and small organelles in a biological sample. To tackle this problem, we propose a filtering strategy with complex wavelet filtering and Anscombe transform for FLIM phasor analysis. This filtering strategy preserves fine structures and reports accurate lifetimes in photon starved FLIM imaging. Moreover, this filter outperforms median filters and makes FLIM imaging with lower laser power and faster imaging possible.


Introduction
Fluorescent Lifetime Imaging Microscopy (FLIM) is an imaging modality that measures the time a fluorescent molecule remains in an excited state before radiatively relaxing back to ground state [1,2]. In time domain FLIM, this is achieved by plotting the fluorescent decay of fluorophores in each pixel with a Time Correlated Single Photon Counting (TCSPC) module [3,4]. However, TCSPC collection generates huge datasets which results in difficulties to visualize and analyze the FLIM data. To simplify analysis and visualization of FLIM, the lifetime data are decomposed to real (G) and imaginary (S) components by Fourier transform. By plotting the real component against the imaginary component and producing a 2D histogram based on their coordinates, we generate a 2-dimensional representation of the lifetime data called the phasor plot [5]. FLIM phasor plots provides an intuitive and easy analysis of the FLIM signal, especially when multiple components are represented in the data [6,7].
One direct application of FLIM is as an additional contrast for separating different fluorescence. Even when different species of fluorophores might have similar values on the photon intensity map, one can group the fluorophores depending on their respective phasor clusters on the phasor plot. Another useful application of FLIM with phasor analysis is metabolic imaging by probing nicotinamide adenine dinucleotide (NADH). The autofluorescence of NADH in living cells and tissues can be used to determine metabolic states. When employed in time lapse imaging, it is possible to monitor shifts in glycolytic versus oxidative states [8,9]. FLIM also provides easy observation of Förster Resonance Energy Transfer (FRET). FRET leads to substantial lifetime changes on the chromophores, which makes observation and analysis of FLIM much more robust compared to the traditional ratio-metric calculations [10,11]. Recently, a combination of Phasor-FLIM with Stimulated emission depletion (STED) microscopy has led to new imaging modalities such as the Leica τ-STED, which enabled an increased STED resolution and elimination of uncorrelated background noise, even at low excitation and STED powers [12,13].
Like with other photon-counting imaging modalities, collecting enough photons for image formation is a big challenge for FLIM. Profiling fluorescence lifetime requires more photon counts than does simple localization of molecules, which in many cases is already a challenge due to the few photons a fluorophore can emit [14]. The tight photon budget is made even more limiting when photobleaching and photo-perturbating cells are of major concern, namely time-lapse experiments [15]. Other groups have attempted to analyze photon-starved FLIM images with different approaches such as global analysis fitting methods [16,17], and Bayesian fittings methods [18][19][20] When applying phasor analysis to FLIM, one of the most common methods to extract meaningful data is to median filter the images [21]. Subjective adjustments are made within the median filter to achieve an optimal processing result while minimizing the number of pixels included in the median value calculation. Although effective in eliminating noise to a certain extent, median filters do not perform well when there are low photon counts. Moreover, median filtering generally diminishes high spatial frequency components such as the edges of features and puncta regardless of photon budget. These features often correspond to membranes, mitochondria, granules, and other fine structures, which are of biological interest [22,23].
To address the need of an image processing technique which preserves the fine details on photon starved phasor FLIM images, we propose a filtering strategy that combines a dual tree complex wavelet filter with an Anscombe transform [24,25]. This new filtering method needs no subjective adjustments by the user to give an optimal result. Our Complex Wavelet Filter (CWF) minimizes the effects of Poisson noise and adopts the strength of complex wavelet filtering in high dynamic range images. We show that CWF drastically outperforms median filters in low photon images. Furthermore, CWF will be especially useful in most situations where available photons are limited, for instance in metabolic and STED imaging.

Materials and methods
Anti-Tom20 AF488 in fixed Cos-7 cells ("Cell 4 color samples", produced by "Center for Microscopy and Image Analysis, University of Zürich, Switzerland") was imaged on a STELLARIS 8 FALCON & SP8 DIVE FALCON (Leica Microsystems, Germany) systems equipped with HyD X and HyD SMD detectors and an HC PL APO CS2 63x/1.40 OIL objective (Leica Microsystems, Germany). Excitation was set to 499 nm via the systems' tunable White Light Laser (WLL) and we acquired the emission from 504 nm to 548 nm with tunable SP detection. A lower intensity threshold was set to 2 photons for all image analyses.
Live-cell metabolic imaging of HEK 293 cell line was collected on an SP8 DIVE FALCON spectral multi-photon FLIM microscope (Leica Microsystems, Germany) and excited with a Spectra Physics Insight 3X ultrafast IR laser at 740 nm, 10 frame repetitions, with 390 nm -510 nm detection and processed using the LAS-X FLIM/FCS module.
τ-STED images were acquired on an SP8 STED FALCON (Leica Microsystems, Germany) using a HC PL APO CS2 93x/1.30 GLYC objective. Mitochondrial membranes of HeLa cells were labeled with Nile Red. Nile Red was excited with 561 nm, the emission range was from 573 nm to and 674 nm and a 775 nm pulsed laser was used for STED.

Equations
A fluorescence decay signature, I(t) with respect to time t of a pixel I, is measured by single photon counting in FLIM. The decay is decomposed into a real part G and an imaginary part S by Fourier transform: (1) where n denotes the chosen harmonic and ω = 2πf in which f is the laser repetition rate. The phasor coordinates G and S are further plotted as the horizontal and vertical axis, and a 2-D histogram is generated to form the phasor plot. Notice every pixel in collected image would have a corresponding phasor coordinate in the phasor plot, which enables selections of pixels in the image based on their phasor coordinates, as seen in Fig. 1. Complex Wavelet Filter (CWF) enables low photon count FLIM datasets to generate similar phasor plots and information extraction as ideal (high photon count) datasets. The grey-scale, intensity images of Tom20 stained mitochondria (a, b) are overlaid with FLIM phasor-selected pixels (red overlay). The unfiltered phasor was calculated from the lifetimes integrated from 50-frames (c) with the selected lifetimes in the phasor (c inset) that map back to the images in a and b and serves as our "ground truth". A photon starved image was created by extracting a single frame from our 50 frame ground truth dataset. That single frame with the unfiltered phasor (d-f) demonstrate how noisy the phasor becomes in the photon starved condition and the corresponding phasor selection (f inset) misses most of the Tom20 pixels (highlighted red in d and e). CWF application to the single photon starved image (g-i) shows that the mitochondria can be clearly visualized through the phasor selection (i and i inset). The phasor insets (c, f, and i) are slightly reduced scale.
The phasor G and S values in photon-starved images usually suffer from noise previously defined as scatter error. Poisson noise is the major contributor to this scatter error [26]. Traditionally, once the G and S values for every pixel are calculated, a median filter or a wavelet filter is applied to the G and S images, respectively. However, median filters degrade high spatial frequency components, even with favorable photon statistics. Directly applying wavelet filters to phasor components does not reconstruct meaningfully in low photon-count regions where structures are heavily obscured by scatter error. We therefore applied CWF. We will refer to these filtering methods as CWF and Real wavelet filtering, to distinguish between the two wavelet filtering methods.
To address this issue, first, we calculate the Fourier coefficients and intensity for every pixel i: which generate three separate images. Second, the Anscombe transform: is applied to each image. Anscombe transform makes the standard deviation of the data approximately constant; and therefore, transforms the G and S coordinate maps to approximate standard Gaussian distributions. After Anscombe transformation, we perform a dual tree complex wavelet transform on the transformed images. The Dual-Tree Complex Wavelet Transform (DT-CWT) generates 6 bands with directional sensitivity for −75°, −45°, −15°, +15°, +45°and +75°for each resolution level. We used the DT-CWT with LeGall 5, 3 tap filters for the first level and Kingsbury quarter sample shift orthogonal (Q-Shift) 10, 10 tap filters for higher levels [27]. The calculation of the wavelet coefficient and a demonstration of the procedure is shown in Figure S1 in Supplement 1.
Finally, the images are transformed back with the inverse Anscombe transform and we calculate the G and the S with:

Simulation results
The CWF wavelet filter was first tested on a simulated dataset along with median filter and the real wavelet filter as shown in Fig. S2 in Supplement 1. The simulated data was created by two components with partial overlap, the two-component having 1.5 ns and 2.5 ns lifetime respectively. The simulated has lower spatial frequency alternations at the periphery of the circle, with increasing spatial frequency moving toward the center. The simulated data were created assuming 0.1 photons per laser pulse with 80 MHz repetition rate and 3 µs pixel dwell time. A ground truth dataset was created by having 500 image accumulation, while different filtering was tested on a dataset with single image accumulation. The phasor plots in Fig. S2 (a-e) in Supplement 1 demonstrated the CWF reconstructed a phasor plot most similar to the ground truth, when the median filter and real wavelet filter still maintained substantial scatter error. If analyzing the G coordinate maps, we can see that the median filter and real wavelet filter fail to reconstruct the structural patterns at high spatial frequency regions, while the CWF was able to maintain the structural pattern to a much higher spatial frequency (Fig. S2f-o in Supplement 1). This was verified in the patterning within the error maps when comparing each filtered image to the 500 images accumulated ground truth (Fig. S2p-y in Supplement 1).
Further, we quantified the Mean Square Error (MSE) of G values with respect to spatial frequencies calculated based on the single frame accumulated image, taking the 500 accumulated image data as ground truth. The results shown in Fig. S3 in Supplement 1 show CWF outperforms the median filter and the real wavelet filter at all spatial frequencies. MSE of median filtering and CWF both increase at higher spatial frequency. However, the median filter increases in error more than CWF as spatial frequency increases. There seems to be a decrease in MSE for the median filter above 0.3 cycles/pixel spatial frequency, which is an artifact due to the breakdown of spatial frequencies in the ground truth. This led to an averaging effect of the G values at locations with the highest spatial frequencies. The median filter happens to reconstruct values similar to flawed regions in ground truth image which lead to the observed decrease in MSE values for the median filter.

Mitochondria imaging of fixed Cos-7 cells
To analyze the performance of CWF, we start with stained cell lines which are the simplest scenario for FLIM imaging. We imaged fixed Cos-7 cells with antibody labeled mitochondria, and created two separate datasets, one ideal dataset with vast number of photons as the ground truth and one simulating the real-life scenario with limited photon budget. The ground truth dataset was created by 50 integrated frames over the same field of view, and a single frame image was segmented out of the 50-frame dataset as our real-life sample. We further processed the 1-frame dataset with the complex wavelet filter along with other filtering methods in hopes to extract similar information as from the ideal dataset.
The phasor plot for the 50-frame dataset, the unfiltered 1 frame dataset and the complex wavelet filtered dataset are shown in Fig. 1(c), 1(f), 1(i) respectively. Comparing Fig. 1(c) and 1(f), we see that under low photon counts, the phasor plot is broadly scattered. In Fig. 1(i), we see that CWF was able to correct for the noise and produce a phasor distribution similar in size, shape, and most importantly in G and S positions compared to the ground truth.
To further demonstrate the validity of the CWF, we performed a pixel grouping of the structural image based on their phasor plot coordinates. As explained in the introduction, the phasor plot is a 2D histogram in nature. We highlighted all pixels above 1/3 of the maximum histogram count in of each phasor plot as an unbiased method of phasor selection. The resulting selections are shown in the insets of Fig. 1(c), 1(f), 1(i) in red, and the selected pixels are shown in Fig. 1(a), 1(d), and 1 g in red as well. For the unfiltered 1-frame dataset, although a large area of the phasor plot was selected, the selection region in the Fig. 1(d) is scattered and very limited, with a lot of structural information missing. However, when the complex wavelet filter was applied, the selected region became much more similar to that of the 50-frame dataset, with more structure covered.
To visualize the filtering effect of CWF along with other filtering strategy on structural details, we showed the G-coordinate maps of the different filtering method in Fig. 2. The unfiltered 1 frame data set contained substantial amount of noise and the G values in some regions deviated quite a bit compared with the 50-frame ground truth dataset. Median filters and real wavelet filters were able to reconstruct higher intensity and lower spatial frequency to limited extents and generated artifacts. In contrast, CWF performed well, especially with higher spatial frequencies, given the limited information in the photon starved sample.
We displayed the error maps of each filtering method by calculating the absolute difference of each pixel between the filtered image and the ground truth in Fig. 2 to display the artifacts of different filtering in detail. Red-shifted colors in the error maps indicate higher error values, and an image with no error would be completely blue, not shown. Median filtering and the real wavelet filter have much higher errors compared to CWF, indicated by the general red-shifted error map colors in Real (Fig. 2(j) and (n)) and Median (Fig. 2(k) and (o)) filtered images. Higher Fig. 2. G-coordinate-and error maps show that CWF achieves better structural reconstruction of FLIM phasor data than other filtering methods. FLIM G-coordinate maps (a-h) are arranged in left-to-right order of increasing error (i-p) highlighting the detail of those errors in the zoomed in views (i-l) of e through h, respectively and pixels with photon counts below 2 were thresholded out. This thresholding was important because pixels with such low counts are not typically considered in analyses. Error maps (i-l) are generated by plotting the absolute difference of the filtered G coordinate maps with respect to the accumulated 50-frame ground truth. The color-scale bars show the colorimetric values of the G coordinate and Mean Square Error (MSE). Single-frame CWF was able to recover the most information, without introducing artifacts (a) and has the least error (i and m) of all the applied filtering methods. The real wavelet filter seems to have diminished detail in the G-coordinate map (b) but surprisingly has a similar error calculation (j and n) to the median filter (c, k, and o). The unfiltered single-frame G-coordinate map is extremely noisy, with heavy deviation from the actual G values in some regions (d) and, not surprisingly, has the most error (l and p). Where the CWF really shines when filtering the higher spatial frequency components, finer features, highlighted by the arrowheads in panels m through p. spatial frequency components as indicated by the arrowheads in Fig. 2(m)-(p) demonstrate CWF can reconstruct FLIM signals in finer features with less error than other methods.
Quantitative analysis of complex wavelet filtering versus median and real wavelet filtering was accomplished through calculating and comparing MSEs of the real component (G) using the 50-frame accumulated image as the ground truth. The MSE value of each filtering is displayed inside the error map in the upper right corner of Fig. 3(i)-(l). Percent denoising was calculated in MSE reduction compared to our ground truth dataset. CWF has more than six-fold decrease in MSE and about 22% better noise reduction than Median Filtering. Fig. 3. CWF provides superior filtering for autofluorescence NADH metabolic live-imaging of HEK-293 Cells. A single-frame of unfiltered raw data (a), median filtered (b), & CWF (c) and within each condition, from left-to-right: phasor plots, phasor color coded selections, the corresponding images, and zoomed regions of interest (from dashed boxes). The twodimensional gradient plot shows the degree of bound NADH (Oxidative Phosphorylation) and free NADH (Glycolysis) in the image overlaid on the intensity values (gray-scale). CWF can extract valid and meaningful information from the noisy data while the median filter generates artifacts along high spatial frequency components such as along the nuclear membrane and creating what looks like structure in the nucleus (b). CWF has smooth reconstruction with valid alignment with the cell structures (c).

Autofluorescence metabolic FLIM imaging of HEK 293
FLIM phasors has been heavily used for NADH metabolic imaging. NADH is an important cofactor heavily involved in metabolism. NADH is also autofluorescent, requiring no addition staining to be fluorescent. Imaging of NADH allows the determination of metabolic states in a sample, as NADH has different lifetimes when bound or unbound to enzymes [28]. When bound, NADH has a longer lifetime, which results in a FLIM phasor signal on the upper left corner of the universal circle on the phasor plot, while free NADH results in a shorter lifetime and signal clustered to the bottom right. Higher percentage of bound NADH usually indicates more oxidative phosphorylation in the cell, while higher free NADH percentage indicates more glycolysis [8,9].
Live cell imaging has its own uncertainly principle that especially applies to metabolic FLIM imaging. Metabolic and physiological perturbations are caused by observation of NADH due to phototoxic, photobleaching, and thermal effects of the 2-photon laser. "Ground truth" FLIM images of autofluorescence are therefore not valid.
However, to analyze the filtering performance of CWF on autofluorescent NADH FLIM, we imaged HEK-293 cells (Fig. 3). NADH in the cytoplasm is more concentrated and has longer lifetimes since they have higher binding ratios to enzymes. In the nucleus, NADH is present with lower concentration and with higher ratio of soluble, unbound NADH. Therefore, two clusters arise in the phasor from the two population of NADH; one that come from cytoplasm and one from the nucleus. These two populations are hardly seen without filtering (Fig. 3(a)) or when median filters are applied (Fig. 3(b)). CWF, on the other hand, reveals both populations and their respective phasor clusters are clearly separated.
We color-coded the intensity images based on phasor distributions after filtering. When no filters are applied (Fig. 3(a)), we see that speckled color variance in the uniform structure like in the nucleus, caused by the scatter error in the phasor plot. The median filter (Fig. 3(b)) partially reduced the scatter error, but we can see structure artifacts forming, especially at edges on the nuclear membrane. Again, this is due to degradation of high spatial frequency components induced by median filters. With CWF (Fig. 3(c)), the nuclear membrane and mitochondrial puncta are well reconstructed in a smooth and realistic manner. We demonstrate that the extracted FLIM information processed by CWF aligns well with the expected lifetime and diffuse distribution of nuclear versus cytoplasmic NADH.

t-STED imaging of Nile-red stained mitochondria
Although STED provides sub-diffraction limit resolution for fluorescence imaging, it suffers from low photon budgets in most scenarios. Moreover, selective deactivation of fluorophores, which is the basis of STED, require high depletion power, which can lead to many problems like photon toxicity [29]. To tackle these problems, researchers have combined phasor FLIM with STED, realizing that valid information of STED are concentrated in certain parts of the phasor cluster. This has led to technologies like t-STED, which further improved the imaging resolution while also lowering the laser power required. However, we still face the problem of having low photon counts, and the degrading effect of high spatial frequency components in traditionally used median filters are fatal.
In Fig. 4, we show t-STED images of mitochondria in live Hela cells stained with Nile red and processed them with different filtering. Figure 4(a) shows the unfiltered image with substantial noise due to the low photon counts and the respective phasor plot also contains a lot of scatter error. Median filters slightly reduce the scatter error of the phasor plot; however, median filters generate erosions on the high frequency components that lead to unnatural edges along the mitochondria membrane ( Fig. 4(b)). CWF further reduces the scatter error in the phasor plot, and the processed membrane are more natural and smoother. This can be further confirmed by a line profiling along the arrows indicated in Fig. 4 (a-c). Structures were totally buried under noise for the raw image. The median filter reconstruction was noisier with more bumps at the edges, while the CWF reconstruction was more uniform and realistic.

Discussion
In this paper, the proposed CWF for FLIM phasors analysis outperforms the traditional median filters for FLIM phasor analysis. Median filtering of images is a relatively simple way to clean up a noisy image and it delivers useful results for a range of applications. However, as the matrices of pixels to perform the median filter grow larger, the sharpness of the image is reduced, and the high spatial frequency components are washed out. This problem is worse in the case of photon starved specimens.
Compared to median filters, we demonstrated that CWF reduces the scatter error of the phasor plot by at least 20%. Moreover, CWF preserved high spatial frequency components of the images well, including mitochondria and membrane, while the median filters generated filtering artifacts. CWF also has more outstanding performance when the photon budget is low, and therefore enables imaging with lower laser power and fast imaging speed.
There are two major changes in our filtering strategy compared to using traditional filters. First, we added an Anscombe transform pair to address the complications of low photon counts. The Anscombe transform mediates the high variance of photon limited images, making standard deviation approximately constant before doing the wavelet transform [25]. It needs to be noted that the Anscombe transform is performed on Fourier coefficients calculated by Eqs. (3) and 4, instead of on G and S maps as median filtering does. In median filtering, G and S coordinate maps are normalized to values ranging from 0 to 1 based on photon intensity counts. Anscombe transformation is therefore non-applicable to median filters. Second, oriented wavelets are useful to address directionality of features through application of the dual tree complex wavelet transform, preserving edges and fine structures [24]. Combining these two changes, we are pushing the preservation of fine structures when only limited photons can be collected for every pixel.
CWF has advantages over traditional wavelet filtering due to CWF anti-aliasing and is shift invariant [30]. Compared to median filters, CWF is also easier to use as it does not require manual adjustment of hyperparameters such as filter size. For a median filter, the optimal filter size will depend on the pixel size of the image and the size of the target imaging structure, which in many cases takes multiple attempts to adjust. This process is eliminated with CWF and therefore is not only more user friendly but also more objective. While phasors have also been used in other imaging modalities like hyperspectral imaging, this approach could also be easily transferable to those analysis as well [26].