Quantitative Chemical Imaging and Unsupervised Analysis Using Hyperspectral Coherent Anti-Stokes Raman Scattering Microscopy

In this work, we report a method to acquire and analyze hyperspectral coherent anti-Stokes Raman scattering (CARS) microscopy images of organic materials and biological samples resulting in an unbiased quantitative chemical analysis. The method employs singular value decomposition on the square root of the CARS intensity, providing an automatic determination of the components above noise, which are retained. Complex CARS susceptibility spectra, which are linear in the chemical composition, are retrieved from the CARS intensity spectra using the causality of the susceptibility by two methods, and their performance is evaluated by comparison with Raman spectra. We use non-negative matrix factorization applied to the imaginary part and the nonresonant real part of the susceptibility with an additional concentration constraint to obtain absolute susceptibility spectra of independently varying chemical components and their absolute concentration. We demonstrate the ability of the method to provide quantitative chemical analysis on known lipid mixtures. We then show the relevance of the method by imaging lipid-rich stem-cell-derived mouse adipocytes as well as differentiated embryonic stem cells with a low density of lipids. We retrieve and visualize the most significant chemical components with spectra given by water, lipid, and proteins segmenting the image into the cell surrounding, lipid droplets, cytosol, and the nucleus, and we reveal the chemical structure of the cells, with details visualized by the projection of the chemical contrast into a few relevant channels.

epi direction, transmitted through the dichroic beamsplitter and a Semrock BLP01-532R-25 long-pass filter, spectrally dispersed by a Horiba Jobin-Yvon iHR550 imaging spectrometer with a 600 l/mm grating, and detected by a Andor Newton DU971N CCD camera, with 2 cm −1 resolution. To enable imaging, the excitation and the detection is confocal, with the detection spatial filter given by the spectrometer input slit (30 µm) and a 3 pixels (48 µm) vertical range on the CCD in a relayed intermediate image plane of the microscope with a 24× magnification. The spectral sensitivity was determined with a calibration source and corrected in the Raman spectra to represent detected photons per second and wavenumber.
The Raman spectra have furthermore been corrected by the factor ν −3 , where ν is the frequency of Raman scattered light, to take into account 1 the photon density of states factor in the spontaneous Raman scattering, which is not present in the CARS susceptibility.

Samples
Polystyrene Polystyrene (PS) grains (Sigma Aldrich, molecular weight ∼ 280000) were dissolved in toluene (5% w/w final concentration) and then drop cast on a microscope slide. A #1 coverslip was placed on the solution, to provide and optically flat structure, prior to solvent evaporation. We measured a resulting film thickness of 8 µm using DIC microscopy.

Glycerol trioleate
Glycerol trioleate (GTO) was sourced from Sigma Aldrich. It is a symmetrical triglyceride derived from glycerol and three units of the unsaturated fatty acid oleic acid. 15 µl of GTO was pipetted into a well created by the 9 mm diameter opening of a 0.12 mm thick adhesive imaging spacer (Grace BioLabs) attached to a microscope slide, and subsequently capped with a #1 coverslip.

Lipid mixtures
A series of mixtures of octanoic acid and α-linolenic acid from Sigma Aldrich were prepared at different volume concentration. The two fatty acids show distinctive vibrational frequencies as the former is a saturated lipid while the latter is polyunsaturated (see Fig. 2). The samples have been prepared in the same format as GTO.

Cells
Mouse 3T3L1 cells were maintained in DMEM/F12 supplemented with 10% (v/v) foetal bovine serum and 2 mM L-glutamine in a humidified incubator (5% CO 2 ). Cells were grown until confluent and differentiated for up to 10 days by addition of 85 nM insulin and 2 nM T3 to 3T3L1 maintenance medium with medium being replaced every 48 hours. IMT11 2 mouse embryonic stem (mES) cells were maintained in a pluripotent state in alpha MEM supplemented with 10% (v/v) foetal bovine serum (batch tested), 10% (v/v) newborn calf serum (batch tested), 10 −4 M beta mercaptoehtanol, 2 mM L glutamine and LIF (1000units/ml) (mES media) in a humidified incubator (5% CO 2 ). The mES cells were differentiated according to a protocol developed by Dani and coauthors . 3 Briefly, mouse ES cells were induced to differentiate by creating embryoid bodies (EBs) via the hanging drop method. This was achieved by spotting 1000 cells/10 µl in a bacteriological dish in mES media without LIF. To create hanging drops bacteriological dishes were inverted atop smaller bacteriological dishes filled with PBS for 2 days. Subsequently were inverted and EB's were treated with mES media supplemented with 10 −6 M all-trans retinoic acid for 3 days, with media being changed daily. EB's were cultivated for a further 2 days in maintenance media and approximately 2-4 EB's/cm 2 transferred to gelatin (0.1%) coated glass slides and allowed to settle for 12 hours before culturing in mES media supplemented with 85 nM insulin, 2 nM Triiodothyronine with media being replaced every 48 hours for up to 23 days.
The cells were subsequently washed in PBS and fixed using 4% (w/v) paraformaldehyde and stored at 4 degrees in PBS/50 units/ml penicillin/streptomycin prior to CARS analy-sis. The cells were routinely screened for mycoplasma using the VenorGeM Mycoplasma Detection kit (Cambio, Cambridge, UK).

Spectral extrapolation for field retrieval algorithms
The fast Fourier transform (FFT) assumes a periodic data set, but the CARS intensity ratio spectraĪ C are not periodic. We create an input vector for FFT which suppresses the resulting artifacts by extrapolating the CARS intensity ratio to zero frequency with a constant equal to the measuredĪ C at the lowest measured frequency. We use a frequency step given by the experimental spectral step and we interpolate, if needed, the measuredĪ C to an even spacing of the frequency including zero. We then mirror the frequency range to negative frequencies, resulting in N 1 frequency points. To avoid wrap-around effects, we then extend the length of the vector to N = 2 m elements symmetric around zero frequency, where m is chosen as the smallest integer so that N > 2N 1 . This extension uses the constant value of the maximum measured frequency. The resulting periodic spectrum is continuous.

Details of the iterative Kramers-Kronig (IKK) method
The IKK method uses the self-consistency iteratioñ with loop index i = 0, 1, .. and damping factor α. During the iterations the RMS σ {i} of the ) over the measured spectral region is used to evaluate the convergence, which is estimated by the ratio r = σ {i+1} /σ {i} . We found that using α = 1, i.e. without convergence control, the iteration can fail to converge for |χ v /χ e | ≤ 2. We therefore use α<1, and we initialize α with 0.

Use of CARS reference in the MKK method
In the MKK method, 4 the phase ϕ of the susceptibility is calculated by generating a temporal response using the inverse Fourier transform (IFT) of ln for positive times and the for negative times. This is equal to using the IFT of ln for positive times as in PCKK and thus provides a correction of ϕ for the effect of I ref , as we have analytically verified. However, in MKK a quantity proportional to the imaginary part of the susceptibility is retrieved by multiplying sin(ϕ) with √ I C . This quantity is dependent on the instrument response T (ω). In PCKK instead, the retrieved susceptibility is normalized to the susceptibility of the reference medium, and is independent of T (ω). To illustrate the differences we show in Fig. S1 the comparison of the retrieved imaginary part of the susceptibility in a lipid droplet close to the CH 2 band using PCKK without phase correction (black line) and MKK (red line) normalized to the maximum of the PCKK result. The results deviate according to the specific T (ω) of our instrument.

Phase retrieval of simulated spectra
To verify the validity of the phase retrieval methods IKK and PCKK and compare their performance with the methods MKK and MEM reported in the literature, we first use simulated spectra exactly adhering to the assumptions of the methods. We use the susceptibilitȳ The spectrum of the simulatedĪ C = |χ| 2 , is shown in Fig. S2a, together with the imaginary parts ℑ(χ) retrieved with MKK, 4 MEM, 5 IKK, and PCKK. The MEM retrieval was performed using the tool available . 6 For the MKK method, which uses we assume a frequency-independent T and |χ ref |. The retrieval methods are applied to the extended spectrum as described in the supplement, except for MEM where we use a singlesided frequency range which resulted in a smaller error. The deviation ∆ of the retrieved ℑ(χ) from the exact ℑ(χ) given by Eq. (2) is shown in Fig. S2c. Off resonance, a weakly frequency dependent ∆ is found for all methods, which increases with wavenumber. The re- trieval of the amplitude of the resonances become less accurate as the CARS ratio increases.
The MKK and PCKK methods coincide since T is constant and the phase offset correction is zero due to the inclusion of zero frequency. We have also studied the dependence of the phase retrieval on the spectral range considered, as detailed in the SM Fig. S3 to Fig. S5.
In general, a smaller spectral range introduces larger deviations. However, the errors are in the percent range which is typically below the noise in experimental spectra, and they are thus not a major limitation. In essence, all methods lead to an acceptable result for this Phase retrieval versus width of measured spectral range Phase retrieval is a spectrally nonlocal operation and thus depends on the spectral range over which the CARS ratioĪ C is known. The different retrieval methods for the susceptibility for a simulated spectrum were compared in the main manuscript for a frequency range of 0-4000 cm −1 . In experiments, the measured frequency range is limited due to excitation and detection constraints, as well as measurement duration limitations. It is therefore relevant to discuss the influence of the measured frequency range on the retrieval fidelity for the different methods, which is shown in Figs. S3-S5, using the procedures discussed in the main manuscript. In general, all used retrieval methods produce an estimate of ℑ(χ) with an error increasing with decreasing wavevector range, in particular for resonances with large |χ v |/|χ e |.
Notably, IKK is most sensitive to the reduction in spectral range. PCKK and MEM show a similar error, with the peculiar property that the retrieval is best for a somewhat limited Singular values which contain a signal contribution larger than the noise then have a value larger than √ 2 times the noise value estimated by this fit. We therefore define i max as the largest number for which the singular values up to i max are larger than √ 2 times the fit, i.e.
they contain a signal contribution larger than the estimated noise. Phase retrieval of a measured lipid droplet spectrum sum of the concentrations to be equal to 1, and convergence to 0.01% was requested. The dashed lines represent the spectra of the pure lipids, while the solid lines are the spectra obtained by MCR using random spectra and concentrations given by a uniform distribution between zero and one as initial guess. We found that the MCR results are dependent on the initial spectra, and we therefore represent the spectra as averages with error bars given by plus and minus their standard deviation over random initial spectra and concentrations.
The error is specifically evident in the α-linolenic acid spectrum, corresponding to the second MCR component, around the peak at 3020 cm −1 . The resulting concentrations of the components are shown in the inset versus the nominal concentration, again with the error bars representing plus and minus their standard deviation over random initial spectra and concentrations. The MCR concentrations deviate significantly from the nominal ones, with errors in the 20% range. This is about an order of magnitude larger than for FSC 3 (see Fig. 2), and the FSC 3 method also does not show a significant dependence on the initial spectra and concentrations, with a standard deviation of about 10 −4 in the concentration over random initial spectra and concentrations. Similar results are obtained using MCR on Raman data as shown in Fig. S8b. In terms of computational effort MCR is about two orders of magnitude slower than FSC 3 . A possible reason for the stronger fluctuations of the MCR results is that it uses the concentration constrain explicitly in the minimization algorithm. Any spatially dependent variation in the signal created for example by the refractive index structure of the signal is disturbing the retrieved spectra. In the FSC 3 method we only rescale the concentrations and spectra globally to minimize the variation of the sum of concentrations from 1 over the whole image, which is not having an explicit influence on the retrieval of the spectra. We can therefore conclude that in typical experimental datasets the spatial intensity variations are too strong to use the concentration constrain explicitly for the spectra retrieval.
Instead of using random initial spectra, they can be obtained by using the evolving factor analysis (EFA) algorithm available within the MCR toolbox. EFA provides estimates for the concentration matrix which can be used to initialize the MCR analysis. However the EFA algorithm is computationally even slower than MCR so that only a sub-ensemble of the data can be analyzed within a reasonable time. The obtained estimation of the sub-ensemble concentration matrix can be used to run the MCR analysis to find a set of spectra which can be used than as initial guess for the MCR on the full dataset. The choice of the subensemble is important, as it has to be statistical relevant and unsupervised. We used 10%  of this behaviour showed that the iterative nature of the IKK can produce spurious parts due to limited convergence, which can be above the noise of the data and thus visible in the higher SVD components. This is specifically relevant at sharp spatial structures where spectral artifacts can appear due to the sequential acquisition of the spectra.
In the SVD images of √ I C shown in Fig. S10, the lipid droplets are surrounded by thin layer of a different singular component. This is an result of the interference between χ e and In the corresponding images of ℑ(χ) (see Fig. S9 for IKK and Fig. 4 for PCKK

SVD basis rotation
The SVD spectra of ℑ(χ) do not correspond to individual chemical components, but rather to differences between chemical components fluctuating independently. The SVD analysis and subsequent visualization helps to determine regions of the samples which present the same color and that can be identified as a single substance. Once a substance is identified one can project the first vector of the matrix U onto the defined spectrum and generate a new set of vectors as new basis for the spectra. To do that we calculate an average spectrum in a manually defined spatial region representing one substance. We rotate U to U ′ = UR with the rotation matrix R such that the rotated basis U ′ has a first vector given by this This rotation process of U' can be repeated in the remaining subspace i.e. excluding the already projected dimensions, until all left vectors of U ′ are determined. We limit this operation to the first i max components of U, i.e. the ones containing signal above noise.
We show an example in Fig. S11, where panel a) gives the color reconstructed ΣV T image after SVD analysis of ℑ(χ) with the corresponding singular spectra in b). Since the region outside the cells is dominated by water, we project the first vector of U onto this spatial region (indicated with solid lines in Fig. S11a), resulting in the D ′ displayed in Fig. S11b.
In the resulting first vector of U ′ (see Fig. S11d) the lipid band at 2850 cm −1 is suppressed.
We then project the second vector of U' onto the spectrum measured on a lipid droplet (see    Similarly to what was observed in the high frequency region, the maximum error of the total concentration E C of about ±0.2 is localized at the lipid droplet edges and in the largest lipid droplets. This is consistent with its origin being the optical aberration introduced by the variation of refractive index between the droplets and the cellular matrix, reducing the signal intensity and thus artificially increasing the concentration in the analysis. It is important to note that in this spectral region the component of highest average concentration, water, does not have any vibrational resonance, and can only bee seen due to the inclusion ofχ e p in the analysis. Indeed, using the NMF on ℑ(χ) alone returns a large concentration error E C ranging from -0.4 to 2.1, and a lacking contrast of the nucleoli.