1. Introduction

Fluorescence hyperspectral imaging is a powerful tool for biological studies as it provides an additional dimension of information compared to filter based multispectral fluorescence microscopes [1, 2]. The expansion in the spectral domain enables accurate classification and quantitative analysis of multiple spectrally overlapping components in biological structures for studying complex biological processes [3, 4] and clinical diagnoses [5, 6]. In order to achieve hyperspectral imaging, a 3 dimensional (3D) data structure (x, y, λ) is typically populated. Traditional hyperspectral imaging techniques typically rely on spatial scanning [7, 8] or spectral scanning [9, 10]. Spatial scanning interrogates the emission profile (λ) of a single pixel or line scanned across the field of view (FoV) to construct the 3D data cube. Conversely, spectral scanning methods image the entire FoV (x, y) on a 2D area detector and probe a single wavelength (λ) at a time. The speed of acquisition is typically limited by the mechanical movement of the dispersive optics and integration time for the detector array, which also introduces 1/f noise into the measurement. Capturing inherent dynamics of living systems such as Ca²⁺ oscillations associated with neuron firing typically requires a frame rate at around 20 Hz [11, 12]. Furthermore, significant blur from cardiovascular motion degrades resolution for in vivo imaging [13]. The bandwidth of current computer I/O interfaces also bottlenecks the achievable frame rate for uncompressed hyperspectral imaging. For a system with 2000 fluorescence channels, 16 bits per channel, 512 × 512 images at 20 frames per second (fps) speed, the overall data throughput is about 20.97 GB/s compared to the theoretical limit of a modern PCI-E 3.0 × 16 interface at 15.76 GB/s. The write speed of a commercially available storage device is typically less than a few gigabytes per second, representing a practical limit on data throughput.

Recently, several snapshot hyperspectral imaging systems were developed to improve the image acquisition speed while still maintaining reasonable data throughput rates [14–17]. Coded aperture snapshot spectral imager (CASSI) directly acquires the projection of a rotated data cube and reconstructs the hyperspectral images using a sparse sampling algorithm [16]. Alternatively, an image mapping spectrometer (IMS) spatially offsets the data cube layer by layer and projects each layer onto an isolated region of a large area detector [17, 18]. Both of these methods provide innovative strategies to overcome the scanning mechanism and bandwidth limitation. However, most of the existing snapshot methods are not directly compatible with multi-photon excited fluorescence microscopy, which significantly increases the depth of penetration for in vivo imaging. Furthermore, images captured using focal plane array detectors generally introduce trade-offs between spatial and spectral resolution, as both sets of information must be encoded on the same fixed number of imaging channels.

In this work, we demonstrated an experimentally simple approach to achieve video-rate hyperspectral imaging for two-photon excited fluorescence (TPEF) through non-descanned spatial/spectral multiplexing (NDSSM) without degraded the spatial resolution. A 16-channel photomultiplier tube (PMT) array (Hamamatsu, H12311-40) was used as spatial-spectral multiplexing detector with fast single photon response time (ns), and paired with a 16-channel digital oscilloscope card (AlazarTech, ATS9416) for fast data digitization. Spatial-spectral multiplexing was achieved by performing non-descanned detection, such that the projections of the spectra illuminated onto the PMT array varied as a function of beam’s scan position. Demultiplexing allowed recovery of spectra with over 2000 fluorescence channels at video rate (17 fps for frame size of 512 × 512 pixels). The data throughput for the 16-channel digital oscilloscope was compressed at 2.24 GB/s with the multiplexing, as opposed to 336 GB/s without compression.

It is interesting to consider the potential change in spatial resolution resulting from the spatial/spectral multiplexing. In CASSI measurements, multiplexing to recover spectral information necessarily comes at the expense of spatial resolution, since the same fixed number of pixels now report on combined spatial and spectral information. In contrast, all the spectral information in NDSSM is contained in the orthogonal spectral axis accessed by the excitation positions in the PMT array, such that hyperspectral images can potentially be acquired without sacrifices in spatial resolution or image quality.

2. Methods

2.1 Experimental Instrumentation

The hyperspectral imaging system is depicted in Fig. 1(A). In brief, a tunable 80 MHz Ti:sapphire femtosecond laser (Spectra-Physics, Mai Tai) was used for excitation. The beam was scanned across the sample using a resonant scanning mirror at 8.8 kHz (EOPC) for the fast-scan axis and a galvanometer mirror (Cambridge-Tech) for the slow-scan. A square aperture was placed in the rear conjugate plane to cut off the turning point of the resonant mirror to avoid sample damage. A 10 × 0.3 NA objective (Nikon) was used to focus the beam down to the sample with the fluorescence signal collected in the epi direction back through the same objective. The laser was tuned in a range of 800 - 1000 nm for different excitations with powers at the sample around 40 - 100 mW. A dichroic mirror and a band pass filter (Chroma, 350 – 700 nm) were used to isolate the fluorescence signals. A 4f lens pair directed the collimated fluorescence onto the center of a transmission diffraction grating (Wasatch Photonics, WP-600/600-25.4). The spectrally separated beam was then focused on to a 16 channel PMT array (Hamamatsu, H12311-40). The 16 channel PMT array electronics were custom built to improve the single photon response time by two orders of magnitude over the original design. The circuit for each channel consisted of a two-stage high-speed operational amplifier circuit connected directly to the output of each element on the PMT array. The output was AC-coupled into a 50 Ω impedance microstrip through a custom printed circuit board (PCB) trace in order to provide impedance matching to the output transmission line. Responses of the PMT array were digitized synchronously with the laser pulses by using a 16 channel PCI-E digital oscilloscope card and remapped into sixteen 512 × 512 images via custom software (MATLAB).

 

Fig. 1 Overview of spatial-spectral multiplexing. (A) Instrument schematic. (B) Illustration of spectra projected on the PMT array from different pixels in one image plane. (C) Spectral reconstruction from a single channel image for two components with different spectra.

Download Full Size | PDF

2.2 Sample preparation

Fluorescein, coumarin 6, sodium dodecyl sulfate, and chloroform were purchased from Sigma-Aldrich, with eFluor 450 obtained from Thermo Fisher. All aqueous solutions were prepared with Milli-Q water (18.2 MΩ cm⁻¹). For emulsion mixtures with two components, stock solutions of 5 μM fluorescein in water and 50 μM coumarin 6 in chloroform were prepared separately. The two stock solutions were mixed equally in amount with ~1mg of sodium dodecyl sulfate. Then the mixture was shaken by hand to produce a stable emulsion mixture for imaging. For emulsion mixtures with three components, stock solutions of fluorescein in water, coumarin 6 in chloroform, and eFluorTM 450 in chloroform were made separately with a concentration of 50 μM for each. Then two two-component emulsion mixtures were prepared separately for eFluor 450 / fluorescein and coumarin 6 / fluorescein system based on the procedure above, and the two two-component emulsion mixtures were mixed in equal parts to produce a three-component emulsion mixture.

For high speed live animal imaging, genetically modified nematode (C. elegans) strains were cultured on agar plates seeded with Escherichia coli at 20 °C. While imaging, C. elegans strains were sealed between two glass coverslips immersed in 0.9 % NaCl physiological saline.

2.3 Spectral reconstruction and component classification

A spectral retrieval algorithm was developed to analyze the raw data. Since the incident angle of the beam onto the diffraction grating varied during the beam scanning, each set of the 16 images from different channels was correlated spatially and spectrally for the same FoV. For each frame, pixels with the same galvanometer mirror position shared the same spectral projection pattern on the PMT array; while pixels along the galvanometer mirror axis produced unique projections as a function of position [Fig. 1(B)]. Therefore, spectrum can be generated by integrating all signals from the same fluorophore along the resonant mirror axis while the galvanometer axis contained wavelength information [Fig. 1(C)].

For an emission spectrum $s_{}^0(\lambda )$ from one fluorophore (i.e., all locations of identical classification), the detected spectrum from a single PMT channel $s_{}^{meas}(\lambda )$ can be described in Eq. (1), in which L is the spectral bandwidth of the PMT channel for a given galvanometer position. The measured spectral intensity $s_j^{meas}(\lambda )$ for a given measurement centered about the wavelength j is limited by L, which is equivalent to the entrance slit size of a spectrometer.

In Eq. (1), δλ corresponds the bandwidth of a single pixel in the final recovered spectrum, which is generally significantly smaller than the bandwidth dictated by the physical dimensions of the PMT L, and s⁰_j corresponds to the ground truth spectrum at the position j. Extending this analysis to each pixel in the image and the corresponding center position of the spectrum upon beam-scanning along the x-axis (galvo), the mathematical relationship described by summation in Eq. (1) corresponds to convolution of the ground truth spectrum with a rectangle function of width L.

Significant signal to noise enhancement in the measured spectrum arises through integration of all the similarly classified locations along the orthogonal y-axis (resonant mirror). For a δλ corresponding to each x-axis position in the image, a single detector can produce a spectrum with up to 512 elements in a 512×512 image.

While the single-detector approach described mathematically in Eq. (2) can in principle recover a complete spectrum at every position, the use of a PMT array has three critical practical benefits: i) measurements with an array produce detected intensity at each location to aid in spectral classification, ii) multiple array elements allow extension of the spectral dynamic range without requiring large displacements of the galvo-axis mirror, and iii) the higher collection efficiency of light and spectral overlap between channels results in signal to noise improvements. The spectral overlap between channels depended on the beam scanning path shown in Fig. 2(A) and the design of the 16-channel PMT array shown in Fig. 2(B). Each single PMT channel had a width of 0.8 mm with 0.2 mm gap between channels. For a FoV of 500 μm × 500 μm, 512 × 512 pixels with a 10 × objective, a narrow band emission was scanned across three adjacent PMT channels with the size of 2.5 mm × 2.5 mm, inducing ~2/3 overlap between the spectral windows for each PMT channel. For a single PMT channel, the detected narrow band emission had a width around 170 pixels along the galvo axis. In the case of no spectral overlap between adjacent detectors, the combined set of 16 detectors could in principle be used to generate a spectrum with 8,192 elements (16 × 512). To minimize edge effects, spectra with 460 spectral channels were used, correspond to a maximum of 7,360 channels. In practice, significant spectral overlap illustrated in Fig. 2(C) enabled spectral stitching and calibration, in which minimization of the squared deviations in the overlapping regions served to calibrate the relative sensitivities of adjacent detectors. As a result, ~2200 independent spectral channels were recovered after the stitching.

 

Fig. 2 (A) Beam scanning path (not the beam path). (B) The projection of spectrum of a narrow band emission on the PMT array. (C) Spectral window calibration with doubling crystal.

Download Full Size | PDF

Second harmonic generation (SHG) from a doubling crystal (BBO) in the object plane was used as the reference as it generated uniformly distributed light across the FoV with a well-defined narrow band spectrum. By tuning the laser to 800 nm, 900 nm and 1000 nm, three references SHG spectra (centered at 400 nm, 450 nm, and 500 nm) were generated for calibration purposes (Fig. 2(C)). As a planar detector was used in the setup, signal projected on the blue and red end on the PMT array had slightly larger spectral separations compared to the center due to the longer projection length. It was consistent with the observation of 553 data points between 400 / 450 nm peaks and 549 data points between 450 / 500 nm peaks. The average channel density was used for the calibration. With a 100 mm lens to focus spectral pattern down to the PMT array, the width of the whole spectral window was about 200 nm. In this case, the rectangular function in Eq. (2) had a width corresponding to 15 nm (170 pixels), which limited the operational spectral resolution. The full width half maximum of the reconstructed SHG peak was 10 nm, compared to 2.5 nm for the excitation source measured by a spectrometer (Thorlabs, CCS175).

The preceding analysis describes the approach taken to recover spectra from pooled measurements of locations with identical composition. Obviously, this analysis first requires prior classification of n composition, which is described herein. A flow chart of the algorithm is shown in Fig. 3. In brief, a threshold was first set to identify signal pixels from the background. Next, a simple Euclidean angle (Eq. (3)) method was used for fluorescence spectra classification. This method is consistent with previous classification approaches [19], selected primarily to minimize the analysis time.

 

Fig. 3 Flowchart for the iterative classification algorithm.

Download Full Size | PDF

Pixels (${\stackrel{\rightharpoonup }{s}}_{x,y}$) with large Euclidean angles relative to the n known reference spectra (${\stackrel{\rightharpoonup }{s}}_{ref}$) (i.e., $\cos \left(\theta \right)$< 0.99) were collectively classified as an “other” component. As a result, discrete masks for each of the n + 1 components were generated. Spectra of all the components were then extracted based on the new discrete classification masks. This process was iterated until 99% of the pixels in the FoV were classified as the same component with the previous iteration. Following classification, spectral recovery for each class was performed as described in the preceding paragraph, with the process iteratively performed until convergence was achieved in both the spectra recovered and the spatial classifications. For data acquired at video-rate, the number of components and their spectra recovered from the first frame was used as the reference for the rest of the frames to reduce the spectral data analysis time. The computation time was dependent on the initial guess and the assumed number of components. For a two-component system, processing a single frame took about 20 to 50 seconds on a 4^th generation hexa-core Intel Core i7 processor running in MATLAB. In order to initiate the iterative process, an initial guess spectrum was produced from the ensemble average of all pixels containing signal in the entire FoV, which served as the initial reference spectrum for a preliminary binary classification.

3. Results and discussion

3.1 Heterogeneous mixtures of fluorescent dye droplets

A fluorescein/coumarin 6 fluorescent dye mixture was used to assess the spectral imaging of the system and the accuracy of the classification algorithm. Figure 4(A) showed a TPEF image from the integrated intensity of all 16 channels for this suspension under 800 nm excitation. A two-component classification image produced the images in Fig. 4(B), with one component labeled in green (component a) and another in blue (component b). It is worth noting that the recovered spectra of the two components were highly overlapping (Fig. 4(B) insert), making them difficult to confidently discriminate using a filter-based fluorescence microscope. The recovered spectra of this hyperspectral imaging system were overlaid with independent measurements from a benchtop fluorimeter (FLS1000, Edinburgh Instruments) with 400 nm excitation (Figs. 4(D) and 4(E)). All spectra were normalized to the maximum intensity. The good agreement between the recovered spectra and the fluorimeter measurements was consistent with component a composed primarily of fluorescein in water and component b primarily coumarin 6 in chloroform, which also agreed with their anticipated respective solubility in water and chloroform. While the fluorescence emission spectra acquired by statistical decomposition were in good qualitative agreement with those obtained from independent measurements with a fluorimeter using pristine solutions, statistically significant deviations were present for both coumarin and fluorescein. In particular, the measurements obtained in the water/chloroform suspensions exhibited broader peaks than those interrogated in pure solutions. The increase in width for the suspensions is attributed to a combination of two effects; i) convolution of the emission spectrum with a rectangle function (vide supra) results in spectral broadening, and ii) the inhomogeneity of the environment in the water/chloroform suspension results in additional inhomogeneous broadening. From the measured SHG response, the contributions from convolution are insufficient to account for the observed differences in peak shapes. Therefore, the primary origin of the deviations observed between the emission spectra in Figs. 4(D) and 4(E) are attributed to spectral broadening associated with inhomogeneities from the greater diversity of environments available for partitioning (e.g., coumarin partitioning into the aqueous phase, fluorescence into chloroform, and both to the surfactant-stabilized interfaces). To study the fluorophores’ distributions at the water-chloroform interface, analysis was repeated assuming a maximum number of three components. Under this case, the boundary of the chloroform droplets was identified as a third component (red in Fig. 4(C)). The recovered spectrum for the third components (component c) was shown in Fig. 4(F). The spectrum of component c (red dot line) was fitted as a linear combination of component a and component b as shown in solid line in Fig. 4(F) with R² = 0.99. Based on this fitting, component a (fluorescein in water) was weighted at 60.6% in the boundary between water and chloroform and component b (coumarin 6 in chloroform) was weighted at 39.1%. This method could provide new strategies for the study of surface interactions between different solvents.

 

Fig. 4 (A) Integrated total fluorescence of a dye mixture (5μM fluorescein in water and 50 μM coumarin 6 in chloroform). (B) Classification under the assumption of two components (a in green and b in blue). Insert: recovered spectra of component a and b. (C) Classification under the assumption of three components (a in green, b in blue and c in red). Insert: recovered spectra for component a, b and c. (D) Overlay of recovered spectrum of component a with the fluorescein emission spectrum from a benchtop fluorimeter, and the difference between them. (E) Overlay of recovered spectrum of component b with the coumarin 6 emission spectrum, and the difference between them. Both plotted spectra were averaged over all the spatial positions of the same classified component. (F) Spectrum of the third component (red dotted line) fitted as a linear combination (black solid line) of the first two components.

Download Full Size | PDF

In Fig. 5, hyperspectral TPEF is demonstrated for measurements of fluorescence photobleaching, highlighting the advantages of high speed hyperspectral microscopy in the analysis of dynamic multicomponent mixtures. The three-component emulsions system of eFluor 450 (chloroform)/coumarin 6 (chloroform)/fluorescein (water) was used for this purpose. A video with both classification images and spectra can be found in the Supplementary Materials (Visualization 1). Figure 5(A) shows a single frame (frame number 55) from the video (green for fluorescein, blue for coumarin 6 and violet for eFluor 450). A target coumarin 6 droplet was circled in red for tracking the dynamic behavior. This frame was set as the start of the dynamic analysis (0 s), which corresponded to 3.2 s in the video. After 21.8 s, most of the coumarin 6 fluorescence was lost from photobleaching. Following the removal of the coumarin 6 signal, the classification of the target droplet was reclassified to ‘eFluor 450’ by the classification algorithm (circled in Fig. 5(B)). The total number of excitations driving photobleaching was plotted Fig. 5(D). In brief, the single-photon peak height distribution of the PMT was measured through a histogram analysis of a time-trace. The mean voltage per photon provided a point of calibration to convert the voltages acquired during imaging to absolute mean numbers of photons. The total number of emitted fluorescence photons for an isotropic emission was calculated based on the solid angle of collection through the objective. Previously measured values for the quantum efficiency of fluorescence for coumarin 6 [20] allowed conversion from the total number of emitted photons to the total number of excitations per unit time. The plot of total number of excitations was then fit to an exponential decay with a time constant of 5.16 ± 0.08 s, consistent with an exponential process routinely observed in photobleaching [21].

 

Fig. 5 Dye mixture with three fluorophores: fluorescein (green) in water, coumarin 6 (blue) and eFluor 450 (violet) in chloroform separately (also see Visualization 1). Target droplet with coumarin 6 in chloroform (circled in red) was tracked from its first appearance in frame 55, set as 0 s (A) and final misclassification in frame 386 with during time 21.8 s (B). (C) Spectra recovered from the target droplet at 0 s (green solid line) and 21.8 s (purple dash line). (D) Photobleaching curve of coumarin 6 in the target droplet (blue dot) with the exponential fit (red line).

Download Full Size | PDF

3.2 In vivo imaging of genetically modified C. elegans

Genetically modified nematodes (C. elegans) were used as a living biological sample for the demonstration of fast hyperspectral imaging. C. elegans strains were obtained from the Caenorhabditis Genetics Center. The main phenotypes of this C. elegans variant are: green fluorescence in the body wall muscle nuclei [ccIs4251, green fluorescence protein (GFP)], green fluorescence in the pharyngeal muscle (mIs12, GFP), dumpy (dpy-17), red fluorescence in the epidermis (frIs7, DsRed), and red fluorescence in pharyngeal muscle (uIs69, mCherry). A video was acquired for the C. elegans with 17 fps, 512 × 512 pixels (Visualization 2, Supplementary Materials) and 800 nm excitation. Since only GFP and DsRed could be accessed by TPEF with 800 nm excitation, a two-component classification was shown in the video and Figs. 6(A) and 6(B). The two components were labeled in red and green separately with the recovered spectra shown in Fig. 6(C). In order to assure the accuracy of this microscopy and classification algorithm, a fluorescence image from a conventional fluorescence microscope (Olympus BX-51) was acquired for the same type of C. elegans with an excitation filter at 460 - 490 nm, as shown in Fig. 6(D). A similar distribution of the two structures with different fluorescence properties was observed, in which the outside epidermis of the C. elegans was red and the pharyngeal muscle in the head was green. In addition, the green fluorescence from the body wall muscle nuclei (ccIs4251) could also be seen in images acquired from our instrument. Compared to one photon excited fluorescence in Fig. 6(D), significantly less out of plane fluorescence was observed with the TPEF. Detailed biological structures could be seen along with the spectra for each component. Component 1 and component 2 were identified to be predominantly GFP and DsRed, respectively. In Visualization 2, an injury was induced in the mitotic region of the C. elegans. Consistent with surface infection, injury, or osmotic stress [22, 23], green fluorescence due to frIs7 mutation was preferentially enhanced at the site of injury.

 

Fig. 6 Two-photon fluorescence image of gene coded C. elegans analyzed without knowing the emission spectra of the fluorophores as a priori at (A) frame 33 and (B) frame 104. Pseudo-color based on different components recovered from custom classification algorithm with green for component 1 and red for component 2 (also see Visualization 2). (C) Recovered spectra for component 1 (green) and 2 (red) for both frame 33 (solid lines) and frame 104 (dash lines). Both plotted spectra were averaged over all the spatial positions of the same classified component. (D) Fluorescence image of gene coded C. elegans with conventional fluorescence microscope under 460 - 490 nm excitation.

Download Full Size | PDF

4. Conclusion

A spatial-spectral multiplexing hyperspectral two-photon fluorescence microscope was developed with over 2000 effective spectral channels in a 200 nm wavelength window and an imaging acquisition rate up to 17 fps. Iterative demultiplexing enabled classification and full recovery of emission spectra with no prior knowledge of the sample. Good agreement was observed between the recovered spectra and those obtained independently using a commercial fluorimeter. The large number of spectral channels of this system made it possible to distinguish fluorophores with similar emission spectra. The high imaging speed made it possible to study highly dynamic multicomponent systems and live biological samples. This instrument is directly compatible with numerous conventional beam-scanning microscope platforms upon addition of a diffraction grating and a multi-channel detector array.