Method of transmission filters to measure emission spectra in strongly scattering media

: We describe a method based on a pair of transmission filters placed in the emission path of a microscope to resolve the emission wavelength of every point in an image. The method can be applied to any type of imaging device that provides the light in the wavelength transmission range of the filters. Unique characteristics of the filter approach are that the light does not need to be collimated and the wavelength response does not depend on the scattering of the sample or tissue. The pair of filters are used to produce the spectral phasor of the transmitted light, which is sufficient to perform spectral deconvolution over a broad wavelength range. The method is sensitive enough to distinguish free and protein-bound NADH and can be used in metabolic studies.


Introduction
Multi-photon excitation is the method of choice for deep imaging in microscopy because once a molecule is excited, the emission can be collected with high efficiency in the presence of scattering using large area detectors and index of refraction matching at each interface in the optical system [1][2][3]. One microscope design that achieves very high efficiency is the DIVER [4]. Once a molecule has been excited; the emitted photon must travel through the highly scattering media which can increase the chances of absorption of the fluorescence due to the effective long path in the scattering medium. While the transport of photons in a uniform scattering media is relatively well understood, the measurement of the emission spectrum is problematic in a microscope because the methods that we use to separate the emission wavelengths have high efficiency only if the light is well collimated [1,4]. One possible solution is use of absorption emission filters not based on interference, which can resolve the wavelength of emission independent of the collimation of the light beam, based on some type of ratiometric measurement [4,5]. One of these filters pair is the sin-cos wavelength dependent combination that can measure the ratio of intensities passing through the filters with high precision [5]. In this work we use this method of emission spectra filtering in combination with the phasor approach that greatly facilitates data analysis. We also analyzed other filters that have a linear dependence on wavelength, for example a pair with increasing-decreasing transmission as a function of wavelength. These filters provide information about the wavelength of emission sufficient to resolve close (in wavelength) spectral emissions.
The resolution of emission wavelengths in biological samples can have important applications, for example, to study metabolic activity within cells using reduced nicotinamide adenine dinucleotide (NADH) fluorescence. NADH is an enzyme cofactor involved in redox metabolic reactions and exists in two spectroscopically distinctive forms: free and protein-bound. The ratio of free-to-protein-bound NADH is related to changes in the cellular energy metabolism [6][7][8]. These two forms of NADH differ largely in lifetimes [9,10], but their emission spectra are very similar [11,12]. Consequently, the free to protein-bound NADH ratio has been extensively studied using Fluorescence Lifetime Imaging Microscopy (FLIM) [6,7,[13][14][15][16], but not as much

Filter method to obtain the spectral phasor
Assuming that we have an ideal sin-cos filter i.e., the transmission of the filter can be described mathematically by the sine and cosine functions or a linear filter combination, we can establish with simulations the best achievable spectral resolution as a function of the total photon counts for each kind of filter (sine-cosine or linear). If the filter response does not follow the mathematical form of the sin-cos functions, a correction can be done to use the filter to provide the correct mathematical shape as described previously [5]. We emphasize here that we are not measuring a spectrum, which is a complex quantity, but rather we measure the spectral phasor, which is a single point on the phasor plot. Of course, this reduces the information, but it allows application of some simple mathematics to obtain a quantity which corresponds to the spectral maximum and depending on the filter used we can also obtain the spectral width. More importantly we can utilize the property of the phasors to resolve linear combination of spectral component in the same pixel of an image [20][21][22][23].
First, we briefly describe the algorithm used to calculate the spectral phasor and the equations used in this work.
G and S are the coordinate of a pixel in the phasor plot. I(λ) is the emission spectrum at a given pixel. The emission spectrum is captured at the bandwidth L (the range of transmission of the filters) and n is the harmonic of the Fourier coefficients in Eq. (1) and (2). We consider only the first harmonic in this work. The transmission filter method for hyperspectral imaging consists of filtering the emission by passing it through a filter with a given wavelength response and measuring the total transmission through the filter. This procedure provides the numerator of Eq. (1) and (2). Note that the spectral phasor calculation is done by the transmission of the filter itself so that Eq. (1) and (2) are the result of single measurement rather than the calculation implicit in Eq. (1) and (2) that would require measurement of the entire spectrum. The filter transmission is normalized to the total emission intensity given by the denominator of Eq. (1) and (2). The resulting G and S coordinates are plotted in the phasor plot with G being the x coordinate and S the y coordinates. We plot the coordinate G and S in the spectral phasor plot which is a polar plot of the coordinates G and S. In an image we typically collect 10 5 -10 6 points, depending on the image size. Note that the method can be used in a raster-scan microscope, confocal microscope or in a camera-based microscope.
The phase of the point in the phasor plot is related to the phasor emission wavelength at that pixel. In the case of the sine-cosine filter combination, the modulus of the phasor point is related to the spectral width of the emission band at the pixel. Some important questions follow: 1) What happens if one pixel of an image has more than one emitting species? 2) How far apart in wavelength must the emissions of two or more species be in the same pixel to be distinguishable? 3) Another type of question is related to the best shape of the filter that is invariant for scattering and that also best separates two spectral bands. In this paper we perform simulations to answer these questions and then we show the performance of real filters in our microscope in the presence of very large scattering. Phasors follow the law of linear combination which means that if in one pixel the emission of 2 or more species is present, the position of the phasor point is the linear combination of the spectral components weighted by the spectral intensity of each component at that pixel.
For an acquisition of a spectral phasor using the sine-cosine filters, we only need to measure the transmission of the sine, the cosine filters, and the total intensity over the same spectral band of the filter as needed to calculate Eq. (1) and Eq. (2). We will perform simulations to determine how much deviations of the ideal shape of the filters can be tolerated, namely resulting in a different point in the phasor plot, (considering photon statistics) and how the total transmission of the filters affects the spectral position and the S/N ratio. We also examine other filters that could have a different response, for example that have a linear wavelength response but can still be used to perform a ratiometric measurement as a function of the light passing through the filters.

Simulations for sin-cos filters
In the simulation we define the G and S coordinates of the spectral phasor that determines the coordinate of one point in the spectral phasor plot using Eq. (1) and Eq. (2). The total count affects the S/N of the wavelength measurement. It affects the spread of the position of the coordinates G and S due to uncorrelated noise in the individual coordinates. However, if there is any kind of correlated intensity fluctuations in both filters, the ratiometric method will compensate for the common mode fluctuations. We assume that in each filter we can measure N photons distributed in wavelength according to the transmission of the filter.
In the sine-cosine filters method we use filters that have the response of the numerator of the expression for G (cosine) and S (sine) over given a wavelength range L, which means we need 2 filters for the data collection. We also need the denominator in the equations for G and S that corresponds to the total signal. The output of the filters can be measured in succession using only one detector or in parallel if at least 3 detectors are available. If we use simultaneous detectors, then the laser fluctuations or any other fluctuation can be compensated, if fluctuations occur in all detectors simultaneously.
If in the simulation we increase the level of noise to 0.1 (Gaussian noise is added to the G and S coordinates), but if we maintain the same parameters for the simulation, we get the phasor plot of Fig. 1

(C).
We also simulated sine-cosine filters with a smaller overall transmission along the spectral band but still the transmission is modulated by the sine and cosine function. We found that the total amount of modulation of the transmission has no effect on the calculation of the spectral phasor position, but that more averaging is needed to obtain a high S/N ratio. The noise in the G and S coordinates only depends on the square root of the total amount of counts in that phasor point in the image. If we change the modulation of the filter from 90% to 10% in the simulations, we still obtain the same spectral phasor position, but the S/N is affected. Note that the modulation of the filter transmission is not affecting per se the resolution. But if a change in modulation due to scattering is also changing the filter transmission at different wavelengths it will change the position of the spectrum in the phasor plot.
While the modulation of the filter transmission does not change the position of the phasor in the plot, it affects the noise because the total filter transmission is less and, as a consequence, the total photons collected are reduced. This consideration is important in scattering media since the filter will receive photons at all angles, but this unavoidable effect will not change the position of the pixel in the phasor plot.
Note that given the G and S equations, a common factor that affects both the numerator and the denominator of the G and S coordinates will not change the phasor position in the phasor plot. So, the laser intensity or the concentration of a substance in the emission pixel will not affect the phasor positions.

Linear transmission filters
Linear filters are commercially available for photography and they are inexpensive. These filters can be fabricated by adding dyes of a given absorption to a solution or a plastic matrix. These filters have the advantage that while using a single detector maximized for collecting scattered light, we can measure transmission of a spectrum through these filters as an integrated quantity (all the wavelength simultaneously) so that the noise is very low.
Examples of these filters with a linear transmission as a function of wavelength are shown in Fig. 2. Let us model these filters with an ideal linear response ( Fig. 2(B)), with a positive and negative slope, respectively. The wavelength range is set from 400 to 600 nm in the simulations to match the wavelength of our linear filter, but any other spectral range will give similar results. The shape of the filter suggests that a ratio of the intensities through the two filters should be sufficient to determine the wavelength of emission as shown on Fig. 2(C). However, the ratio is a non-linear quantity, i.e., it is not bound in a range and the ratio is not additive. Instead, we will use the GP function defined by the following equation: Where g f (orange) and s f (blue in Fig. 2) are the transmission of the two linear filters in the emission bandwidth as measured with a spectrophotometer. The GP (since filter transmissions are positive) is a mathematical quantity that is limited to the range -1 to +1. If the denominator is the total intensity, the GP function is additive. The plot of the GP for the linear filter GP transmission is shown in Fig. 2(C). It is a straight line that starts at -1 at 400 nm and reaches +1 at 600 nm.
The GP function is perfectly linear as a function of wavelength (for an ideal linear filter). The value of the GP directly provides the value of the wavelength of the emission due to the fluorescence in a pixel. Since the GP formula is independent of the intensity (any common factor in the numerator and denominator in the GP equation cancels out), the total intensity measurement is not needed for the linear filters, the GP is sufficient to determine the emission wavelength. The modulation of the G and S coordinates is instead dependent on the total filter transmission.
In terms of the original equations, Eq. (1) and (2) can be rewritten in terms of the total filter transmission (s f +g f ) as These equations show that the measurement of the filter transmissions, g f and s f , respectively, are sufficient to obtain the spectral phasor coordinates weighed by the total intensity. These equations are valid for any type of linear filters when the values of g f and s f can be directly measured. Any filter that is monotonic up, in pair with a filter that is monotonic down, will satisfy the requirement that the GP value is sufficient to determine the emission wavelength. However, in the case of an ideal linear filter the up and down filters are complementary and only one filter needs to be used in addition to the total intensity.

Simulation of the response of a linear filter to Gaussian spectral bands
A Gaussian function (the spectrum, sp1(λ)) passing through the "down" filter (ft1(λ)) will give the following output If the spectrum is very narrow and at the lower wavelength of the filter, G fdown transmission is 1 and at the other extreme of the wavelength range the G fdown is zero. The value of G fdown is between 0 and 1 and it can be used to map the phase between 0 and 360 degrees by Eq. (6) where the angle θ in the phasor plot is in radians.
The modulation of the classical (sine and cosine) spectral phasor is sensitive to the wavelength bandwidth of the original spectrum, where high modulation (close to 1) is reflective of a narrow band, and low modulation (towards zero) is due to a very wide band [20][21][22][23]. This situation is not achieved with a single linear filter, because the integral of the linear filter is the same, independent of the position of the Gaussian band in the filter range, so we simply define the modulation to be unity, or the expression in Eq. (7) which is always unity.

Using one down (or up) filter and the total intensity
If the total intensity is measured, then Eq. (1) can be written as in Eq. (8), where Sum is the total intensity, which means there are no filters in the optical path.
This approach is the easiest implementation of the filter based spectral phasor method. If the normalized transmission is obtained as shown by the normalizations of the down filter (orange line in Fig. 2) and up filter (blue line in Fig. 2) then the output of the down (or up) filter is directly proportional to the wavelength.

Exponential transmission filters
This type of filter is shown in Fig. 3 for the case where the transmission is described by an exponential dependence on wavelength, the value of 1 corresponding to the lowest wavelength and the value of zero corresponding to the maximum wavelength decaying exponential paired with another filter described by an increasing exponential (blue and orange lines in Fig. 3(A)). In this case the GP function is not linear (black line in Fig. 3(B)) but will still map uniquely to the emission wavelength. For this filter, the total intensity is 1 and is constant at all wavelengths (light blue line in Fig. 3(B)). A separate measurement of the total intensity is not needed for this filter, since it is constant and independent of wavelength, but it will depend on the laser power and sample concentration.

Simulation of 3 Gaussian spectra with different centers in the same pixel
We examine 3 different models of up-down filter when transmitting the same spectra. The output is identical for the 3 filters examined as shown in Fig. 4.

Fig. 4. A-B)
Sin-cos filter (using transmission in a restricted spectral region to simulate a "linear transmission" of Sin, Cos and the total intensity. This result is obtained using only part of the sin-cos filter. It could help in some arrangements when the total wavelength range can be limited. C-D) Transmission of the same spectrum through ideal linear-up and linear-down filters. E-F) Light is passing though exponential-up and exponential-down filters. In the last 2 cases, only the light transmitted by two filters is needed to calculate the spectral phasor.

Exponential filter (or any monotonic or down filter)
In this case we use the G and S coordinates to calculate the spectral position. The same spectrum position is recovered perfectly (compared to linear filters) but the intensity measurement could not be obtained for any monotonic curved filter. The same reasoning applies to the commercial Kodak and Lee photographic filters.
As summary of this simulation on the recovery of Gaussian spectra, all 3 methods: "sincos-intensity", "linear up down ramp" and "arbitrary monotonic up and down" gave the same positions of the spectral phasors, but the monotonic curved filters only use one filter transmission and the total intensity, which can be an advantage in the use of the monotonic filter method in a confocal system with only 2 detection channels. Of course, an effective filter-factor (the so-called g-factor in spectroscopy) must be determined to calibrate the different sensitivities of the two channels used. This procedure is done using the monochromator calibration.

Simulation of images
Until now we have not examined the response of the 3 types of filters to noise in each pixel. For this purpose, we perform simulations of images with different levels of noise because the statistics in the recovery of images is easier to appreciate and quantify given the large amount of data in an image. Using the dual detector method, the noise that is common to both filters cancels out. Only photon noise that can be different in the two detectors will affect the position of a spectral component in the phasor plot. We simulated data with a noise of 0.05 which corresponds to about 400 photons in total collected in a filter. Then we simulated images with different wavelengths of emission. We simulated spectral data for an image with or without noise and we examine the image using the 3 simulated filter sets.
If noise is added equally to the intensity in both channels it has no effect on the position of the spectral phasor. There is, however, an effect on the wavelength recovery if the photons are few and if the photon noise in the two filters is not the same due to random fluctuations. In this case the position of the spectral phasor moves in the phasor plot in the same manner as already demonstrated with simulation of a single spectral phasor. In principle, the filter methods based on measuring only two quantities (filter transmission and intensity at each pixel) are very robust, and they are not sensitive to the signal levels that affect both channels.

Measurements using linear-up filter and total intensity
We made a colored linear-up transmission filter using dyes (for details see Materials). The spectral transmission of the filter is shown in Fig. 5(A). The spectral range used in our experiments was from 400 to 600 nm, where the filter is mostly linear. Due to the non-ideal linearity of the filter, when plotting the data in the phasor plot, we corrected the phase using a polynomial interpolation of the data as shown in Fig. 5(B).
To test our filter, images were obtained using the linear-up filter and the total intensity (no transmission filter, acquired sequentially) with collimated monochromatic light at 20 different bands in the range 400 to 590 nm. Figure 5(C) shows a calibrated phasor plot with the monochromatic light data, where we can see the correct distribution of light of different wavelengths in the phasor plot. To evaluate the performance of the method with scattered light, we repeated the monochromatic light measurements adding a ground glass diffusor in the optical path ( Fig. 5(D)), and the results obtained were the same as using collimated light. It is important to note that colored filters, unlike interference filters, are not sensitive to incident light direction.
To further evaluate the performance of our method with scattered light, we have imaged a mixture of blue and yellow-green fluorescent microspheres with and without scattering ground glass placed in the optical path between the sample and the detector. Figure 6 shows that corresponding phasor plots for images taken with and without scatter are the same. In both cases, using the reciprocity property of phasor plots and images, we are able to separate blue and yellow-green fluorescent microspheres in images (Fig. 6(B),(E)) using green and red cursors in the phasor plot ( Fig. 6(C),(F)).
Then, we tested the method with images of solutions of Rhodamine 110 and Rhodamine B, which have a maxima in the emission spectrum of ∼530 nm and ∼575 nm, respectively [24]. It is shown in Fig. 7(A) the distribution of pixels in the phasor plot of Rhodamine 110 (selected with blue cursor) and of Rhodamine B (selected with red cursor). Rhodamines appeared at the expected positions considering that we are plotting the center of mass of the emission spectrum. Also, we acquired images of a 1:1 mixture of both Rhodamines (selected with green cursor Fig. 7(A)) which appears in the phasor plot in the middle of the individual Rhodamines angular positions, as expected. Note that in this case the phasor does not follow the properties of linear  addition. The spectral phasor has the average phasor position and a spread or distribution around that position, which is plotted in terms of the normalized distribution histogram (occurrences, Fig. 7(B),(C)). The extent of the phase distribution in the phasor plot is user selectable. Finally, we wanted to demonstrate the ability of our method to discriminate between the two states of NADH (free or protein-bound) since their emission spectra are very similar. Emission maximum for free NADH is ∼470 nm in water, while for NADH bound to Lactate Dehydrogenase (LDH) it is ∼440 nm, hence only a few nm apart [12]. These data may vary in literature depending on how measurements were done [18,25,26]. We acquired images of free and LDH-bound NADH in solution, using a 2-photon laser set to 740 nm for excitation. The distribution of pixels of the two emitting species appeared in different positions in spectral phasor plot (Fig. 8(A) top), the distribution of pixels corresponding to free NADH were highlighted with a red cursor and the ones for LDH-bound NADH with a green cursor. Free NADH appears at longer wavelengths than bound NADH as previously reported [12]. Using FLIM we corroborated the free (0.4 ns) and bound states of NADH (3.4 ns, 80% of bound protein) [13,19].
Then, we applied our method to the autofluorescence of living NIH3T3 cells as shown in Fig. 8(B). In order to collect the NADH signal while avoiding FAD emission (λmax emission 540 nm) we added a bandpass filter (400-500 nm) in the emission optical path. It should be noted that NADH autofluorescence and its phosphorylated version NADPH cannot be spectrally separated, so we are measuring a mixture of both, and their combined fluorescence is referred as NAD(P)H [25,27]. Using the positions in the phasor plot of free and protein-bound NADH established with the solutions, we color-coded the pixel distribution in the phasor plot ( Fig. 8(G)), and each pixel of the image was painted accordingly (Fig. 8(D)). In Fig. 8(D) we can observe that the pixels with the higher contributions of free NAD(P)H (red-pink) are in the nucleus and the ones with higher contributions of bound-NAD(P)H (cyan-white) are in the mitochondria (which appear as bright structures in the cytoplasm), as has been previously reported using FLIM [28][29][30]. We corroborated our results with FLIM measurements, since FLIM has been widely used for this type of analysis, and we had the capability to do FLIM measurements simultaneously with the transmission filters, we acquired FLIM images to compare the results in the same cell ( Fig. 8(E)-(H)). Also, we used a commercial spectral detector with the Zeiss LSM 880 microscope to measure NADH spectra in solutions (Fig. 8(A) bottom) and in live cells (Fig. 8(C),(F),(I)). It also shows that there is a small difference in the emission spectrum of free and protein-bound NADH, and the positions where they appear in the phasor plot are close to each other. The information that we were able to extract from these measurements is comparable with the one obtained by our method. Both methods show higher contribution of free NAD(P)H in nuclei while bound-NAD(P)H is mostly present in mitochondria. Mitochondria in living cells are in constant movement, which may affect how they appear in the image, especially in our linear filter method that uses two sequential images. However, both methods give similar results.

Discussion and conclusion
Using the transmission filter method, the spectral information, which can be a complex curve, is transformed into a point in the phasor plot. This transformation results in a reduction of the information but at the same time we increase the spectral separation of the method at each pixel of the image, with very high resolution (less than 1 nm in wavelength). This resolution cannot be reached using conventional band pass filters. We have shown that our method separates scattered light of different wavelengths perfectly well, which can be especially useful for applications in highly scattered tissue samples or in vivo experiments.
The free to protein-bound NADH ratio has been correlated with metabolic changes in living cells and tissues and can be an important application for our method. The results we have obtained for this biological problem, using our method, are shown in Fig. 8. We were able to separate free from protein-bound NADH in solution and in living cells using spectral information obtained with our linear filter method. These results were similar to those obtained using a commercial spectral detector or FLIM.
Although the filter-up can provide very high spectral separation, since the modulation for this filter is treated as a constant, the resulting phasor cannot be used to obtain a linear combination of species in the same pixel. This issue is clearly a limitation for these up-filters. However, if the spectral components are spatially separated then the resolution of components can still be achieved.
It has been suggested that using the harmonic content very complex spectra could be analyzed [31,32].
A specific application of the sin-cos filter is with cameras. The filters can be inserted in a rotating wheel or using beam splitters to separate the light bands. Of course, if there is strong scattering, an image on a camera cannot be formed.
The most interesting application of the filter method to separate spectral components is for strongly scattering biological materials as we mentioned. In this case the image is obtained by raster scanning using multiphoton excitation and the phasor method can be used for a thick sample. This is a straightforward application and is the motivation for this work. In the case of thin samples, the method should be compared with other available methods based on prisms, gratings, and tunable interference filter.
The main advantage of the up-or sin-cos filter method is the very high spectral separation that can be achieved. There are also advantages from the point of view of S/N since the filter method measures all spectra at the same time using a hardware approach so that a large number of photons can be collected in a short time.
The filter method is also inexpensive. The pair of custom-made sin-cos filters costs about $2000 which is much less expensive than any hyperspectral camera or multiple anode detectors. The algorithms are available in our software (freely available from our website at www.lfd/.uci.edu/software) or the algorithm can be easily implemented in Matlab or Python.