1.

Introduction

Raman spectroscopy is a powerful analytical method due to its chemical bond specificity. The molecular vibrations of each substance produce a characteristic fingerprint spectrum that can be used to determine the chemical and structural composition of the sample. Applications in biomedical research are wide-ranging, including the diagnosis of cancer,^{1, 2, 3} atherosclerosis,⁴ and blood analysis.⁵ The major advantages of Raman spectroscopy compared to fluorescence applications are that Raman signals do not bleach, and that Raman bands are spectrally narrow, which improves the separability of the signals.

A basic challenge for Raman spectroscopic measurements in turbid media is the signal quantification. The spectral line shape and overall intensity of the measured spectra depend on the absorption and scattering properties of the sample. In the case of biological tissue, absorption arises from chromophores such as hemoglobin and melanin. Light scattering in tissue occurs due to the refractive index inhomogeneities of the cell nuclei, cell organelles, or collagen fibers. Due to different biochemical tissue compositions and microstructures, the optical properties of tissues may vary strongly between tissues, tissue sites, and individuals. Since a contrast in optical properties can, on its own, indicate pathological changes such as tumor angiogenesis or dysplasia, extensive research has been carried out on the measurement of optical properties for diagnostic purposes.^{6, 7, 8, 9} Nevertheless, most Raman spectroscopic tissue measurements are not corrected for these influences, which hinders the quantitative comparison of the measured Raman spectra.

Since the measured Raman scattered light and purely elastically scattered light are both influenced by variations in the optical properties, correction methods have been suggested based on diffuse reflectance measurements. Some approaches have been developed for applications where either absorption or scattering is assumed to change^{10, 11} and are therefore not suitable for in vivo applications, where independent and simultaneous variations of the optical properties are common. The first theoretical and experimental studies on the correction of Raman signals in biological tissue were published recently^{12, 13, 14} for a parameter range of absorption coefficient ${\mu }_a$ and scattering coefficient ${\mu }_s$ , which are typical for measurements in the near infrared wavelength range (NIR). The correction functions utilize the elastic diffuse reflectance, which is measured in the same excitation-detection geometry as the Raman spectrum. However, knowledge of the total attenuation coefficient given by ${\mu }_t={\mu }_a+{\mu }_s$ is required,^{12, 13} or the reduced scattering coefficient ${\mu }_s^{\prime }$ is used.¹⁴ However, the determination of optical properties from a single elastic reflectance signal requires an elastic diffuse reflectance spectrum, together with sample specific assumptions regarding the wavelength dependency of the absorption and scattering coefficients,¹⁵ which reduces the method’s flexibility.

The aim of our investigation was to develop a method that could correct Raman signals for the influence of the optical properties in a large parameter range without a priori assumptions of the absorption and scattering properties. Therefore, we combined Raman spectroscopy with spatially resolved reflectance measurements.¹⁶ By performing measurements on various homogenous tissue phantoms as well as Monte Carlo simulations, we investigated how the Raman signal depends on the optical properties. To include parameters typical for measurements in the NIR as well as in the visible wavelength range (VIS), ${\mu }_a$ was varied from $0.01\text{to}10{\mathrm{mm}}^{-1}$ and ${\mu }_s^{\prime }$ from $0.1\text{to}10{\mathrm{mm}}^{-1}$ . These investigations exceed the parameter range of previous studies.^{12, 13, 14} To correct the Raman signals for the influence of the optical properties, we investigated possible correction functions that are dependent on the elastic diffuse reflectance. We also implemented a correction function based directly on the Monte Carlo simulation of the Raman signals.

2.

Experimental Setup and Monte Carlo Simulation

In this section, we explain the Monte Carlo code that was used to simulate the Raman signal and describe our combined Raman and spatially resolved reflectance setup as well as the production of the tissue phantoms.

2.1.

Monte Carlo Model for the Simulation of Raman and Reflectance Signals

In accordance with our experiments, we assumed that the absorption coefficient ${\mu }_a$ and the coefficient for Raman scattering ${\mu }_{\mathit{\text{Raman}}}$ were constant throughout the sample. Because ${\mu }_{\mathit{\text{Raman}}}$ was very small, multiple Raman scattering was neglected and the Raman signal was assumed to be proportional to ${\mu }_{\mathit{\text{Raman}}}$ . We further assumed a Henyey-Greenstein phase function¹² and an anisotropy factor $g=0.8$ for the elastic scattering processes.

We defined a flux density ${\Phi }_{\mathit{ex}}$ of photons that were not Raman scattered and a flux density ${\Phi }_{\mathit{\text{Raman}}}$ of Raman-scattered photons. Since multiple Raman scattering was neglected, both densities could be calculated with the conventional Monte Carlo method described in Ref. 17, which takes into account elastic scattering and absorption. Only the source term was different for both flux densities. For ${\Phi }_{\mathit{ex}}$ , the source term was a boundary condition at the sample surface that corresponded to the excitation light entering the sample. For ${\Phi }_{\mathit{\text{Raman}}}$ , the source term was the number density of Raman scattering events given by $n_{\mathit{\text{Raman}}}={\Phi }_{\mathit{ex}}\times {\mu }_{\mathit{\text{Raman}}}$ , which was proportional to the number density of the absorbed excitation photons $n_{\mathit{abs}}={\Phi }_{\mathit{ex}}{\mu }_a=({\mu }_a∕{\mu }_{\mathit{\text{Raman}}})\times n_{\mathit{\text{Raman}}}$ . Therefore, $n_{\mathit{abs}}$ was used as the source density to calculate ${\Phi }_{\mathit{\text{Raman}}}$ and scale the results with ${\mu }_{\mathit{\text{Raman}}}∕{\mu }_a$ . In the Monte Carlo calculation for ${\Phi }_{\mathit{ex}}$ , the density $n_{\mathit{abs}}$ was given by the ensemble of absorption process positions. Thus, the following scheme was used to calculate ${\Phi }_{\mathit{\text{Raman}}}$ : At each position where an excitation photon was absorbed, a Raman scattered photon was launched with an isotropic initial direction. The photon propagated with optical parameters that corresponded to the Stokes shifted wavelength. A fraction $N_{\mathit{\text{Raman}}}$ of these photons was remitted at the sample surface and reached a Raman detector, which was defined as a circular area on the sample surface. Normalization with the number of incident excitation photons $N_{\mathit{ex}}$ and scaling led to the detected Raman signal ${\mathit{Ram}}_T=(N_{\mathit{\text{Raman}}}∕N_{\mathit{ex}})\times ({\mu }_{\mathit{\text{Raman}}}∕{\mu }_a)$ . In this paper, ${\mathit{Ram}}_T$ corresponds to the calculated or measured relative Raman signal of a Raman scatterer embedded in a turbid matrix. The elastic diffuse reflectance of the sample was determined by the number of excitation photons $N_{\mathit{elastic}}$ detected at the sample surface. In the results shown below, the relative elastic diffuse reflectance signal is given by $N_{\mathit{elastic}}∕N_{\mathit{ex}}$ .

The size and position of the excitation and detector area corresponds to the images of the fiber endfaces on the sample surface in Fig. 1 . For simplicity, the fiber bundle for the Raman detection in Fig. 2 was assumed to correspond to one fiber with a $0.7\text{-}\mathrm{mm}$ radius, and the slant of excitation and receiving optics were neglected. The propagation direction of the remitted photons were analyzed by taking into account the finite numerical aperture (NA) of the imaging optics. For our application, we used $\mathrm{NA}=1$ and scaled the results by a constant factor accordingly. We used $10^5$ photons for ${\mu }_a>0.1{\mathrm{mm}}^{-1}$ , leading to a computation time of a few seconds for each set of parameters but $10^6$ photons for ${\mu }_a<0.1{\mathrm{mm}}^{-1}$ .

Fig. 1

Monte Carlo simulation for Raman measurements in turbid media illustrated by photon trajectories. Only the trajectories of Raman photons that reached the detector in our measurement geometry (excitation spot with radius $r_{\mathit{ex}}=1.5\mathrm{mm}$ and detection spot with radius $r_{\mathit{det}}=0.7\mathrm{mm}$ ) are plotted. An isotropic Raman scattering event occurred somewhere on each of the shown trajectories. (a) Schematic for three detected Raman photons. (b) Simulation for ${\mu }_a=1{\mathrm{mm}}^{-1}$ and ${\mu }_s^{\prime }=10{\mathrm{mm}}^{-1}$ . (c) Simulation for ${\mu }_a=0.1{\mathrm{mm}}^{-1}$ and ${\mu }_s^{\prime }=1{\mathrm{mm}}^{-1}$ . For simplicity, ${\mu }_a$ and ${\mu }_s^{\prime }$ correspond to the mean of the optical properties at the Raman and Stokes wavelength. The trajectories in (b) and (c) correspond to $10^3$ incident photons, in contrast to the $10^5$ to $10^6$ photons used to calculate our data.

Fig. 2

Experimental setup. S1 and S2: shutters; DM: dichroic mirror; BS: cover slide as beamsplitter; FL: focusing lens; PD1: photodiode for laser reference; PD2 through 4: photodiodes for each ring of reflectance detection fibers of OF3; OF1: optical fiber for Raman excitation; OF2: fiber bundle for Raman detection; and OF3: fiber bundle for spatially resolved reflectance excitation and detection. The central illumination fiber is surrounded by concentric rings of detection fibers at distances $r=0.4$ , 0.8, and $1.3\mathrm{mm}$ from the center. L1: lens system with magnification $M=5$ , including laser cleanup filter; L2: lens system with $M=1.5$ , including long-pass filter; L3: lens system with $M=1$ ; L4: lens system with $M=1.8$ for coupling to spectrometer (SP).

3.

Dependency of the Raman Signal on the Macroscopic Optical Properties

First we investigated how the Raman peak at $1524{\mathrm{cm}}^{-1}$ [ $\mathrm{C}\mathrm{C}$ bonds, see Fig. 3 ] depends on the optical properties of the silicone phantom, which correspond to typical values for skin in the VIS.⁷ The Raman signal dropped roughly proportional to $1∕{\mu }_a$ , while the influence of ${\mu }_s^{\prime }$ was fairly small, see Fig. 3. The optical properties of the phantoms are the mean optical coefficients at excitation and Stokes wavelength. These coefficients were measured with the integrating sphere spectrometer, since the whole parameter range was not accessible using our implementation of the spatially resolved reflectance setup. Since the silicone color paste contributed to the scattering coefficient of the phantom, ${\mu }_s^{\prime }$ of the phantoms was not constant within one series of increasing absorber concentration.

Fig. 3

(a) Raw spectra of silicone phantoms containing the same concentration of $\beta$ -carotene but with different ${\mu }_a$ of the silicone matrix ( ${\mu }_s^{\prime }=2+{\mathrm{mm}}^{-1}$ , where “+” indicates that ${\mu }_s^{\prime }$ was not constant within one series but increased slightly with the absorber concentration). The baselines are aligned, and ${\mu }_a$ and ${\mu }_s^{\prime }$ correspond to the mean of the optical properties at the Raman and Stokes wavelength in units of $1∕\mathrm{mm}$ . (b) Comparison of experimental ${\mathit{Ram}}_T$ ( $\mathrm{C}\mathrm{C}$ bonds at $1524{\mathrm{cm}}^{-1}$ , filled symbols) and calculated ${\mathit{Ram}}_T$ (lines or open symbols, fitted to experimental values) dependent on ${\mu }_a$ for various series of ${\mu }_s^{\prime }$ . The Raman signal error bars represent the standard deviation (SD) of five Raman measurements at different positions. The optical parameter error bars correspond to the uncertainty of the optical property measurements with the integrating sphere.

Next we compared our phantom measurements with simulations [Fig. 3]. The simulation results of the Raman signal were all scaled with the same factor to fit the experimental results. The factor includes the unknown in vivo Raman cross section as well as setup-specific constants. The trends of the experimental values corresponded well to the simulations. However, the experimental and simulated values did not match for all phantoms within the given errors, which could be explained by errors incurred by variations in the sample preparation.

We further performed Monte Carlo simulations to investigate the dependency of the Raman signal on ${\mu }_a$ and ${\mu }_s^{\prime }$ in a larger parameter range, as with the phantom measurements (Fig. 4 ). For our measurement geometry, ${\mathit{Ram}}_T$ was roughly proportional to $1∕{\mu }_a$ in the VIS ( ${\mu }_a>0.1{\mathrm{mm}}^{-1}$ and ${\mu }_s^{\prime }>1{\mathrm{mm}}^{-1}$ ). For smaller ${\mu }_a$ (typical for the NIR), the curves deviated significantly from this behavior. When the same simulation was performed for a pencil-beam excitation and a semi-infinite detection geometry, which is a typical imaging geometry, the $1∕{\mu }_a$ behavior continued into the NIR parameter range (data not shown). When the excitation and detection geometry corresponded to a $200\text{-}\mu \mathrm{m}$ fiber, the ${\mu }_a$ dependency was weaker, as for the measurement geometry used here (data not shown). This can be explained by the influence of the measurement geometry on the pathlength distribution of the detected Raman photons, which determines the influence of the optical properties on the Raman scattered light. For our measurement geometry, the influence of ${\mu }_s^{\prime }$ in the NIR is larger than the VIS.

Fig. 4

Dependency of the simulated Raman signal ${\mathit{Ram}}_T$ on the optical properties of the tissue matrix for our measurement geometry, where ${\mu }_a$ and ${\mu }_s^{\prime }$ correspond to the mean of the optical properties at Raman and Stokes wavelength in units of $1∕\mathrm{mm}$ . (a) ${\mu }_a$ dependency: ${\mu }_s^{\prime }$ is increasing from bottom to top. (b) ${\mu }_s^{\prime }$ dependency: ${\mu }_a$ is increasing from top to bottom. ${\mathit{Ram}}_T$ is defined as $N_{\mathit{\text{Raman}}}∕N_{\mathit{ex}}$ and is normalized by the unknown ${\mu }_{\mathit{\text{Raman}}}$ . In the VIS, ${\mathit{Ram}}_T$ is roughly proportional to $1∕{\mu }_a$ . The influence of ${\mu }_s^{\prime }$ is smaller but not neglible.

In addition to the Raman signal, we also analyzed the elastic diffuse reflectance of the laser light, which was extracted from the same spectra as the Raman signals (Fig. 5 ). The elastic signal was much more sensitive to variations of ${\mu }_s^{\prime }$ and less sensitive to ${\mu }_a$ . As a consequence, a simple ratio of Raman and elastic diffuse reflectance signal still depended on ${\mu }_a$ and ${\mu }_s^{\prime }$ and therefore was not sufficient for correction in this parameter range, which covered tissue optical properties in the VIS.

Fig. 5

(a) Comparison of experimental (filled symbols) and calculated results (lines, open symbols) for the dependency of the elastic diffuse reflectance signal $\left(R\right)$ on ${\mu }_a$ for various series on ${\mu }_s^{\prime }$ , where ${\mu }_a$ and ${\mu }_s^{\prime }$ correspond to the mean of the optical properties at the Raman and Stokes wavelength in units of $1∕\mathrm{mm}$ . The data for ${\mu }_s^{\prime }=4{\mathrm{mm}}^{-1}$ and ${\mu }_s^{\prime }=6{\mathrm{mm}}^{-1}$ are not included in the graph for clarity. (b) Comparison of experimental and simulated results (fitted to experiments) for the dependency of the Raman signal $\left({\mathit{Ram}}_T\right)$ and elastic diffuse reflectance signal $\left(R\right)$ on ${\mu }_s^{\prime }$ for ${\mu }_a=0.25{\mathrm{mm}}^{-1}$ . The elastic scattering signal was attenuated by a long-pass filter $(\approx \mathrm{OD}6)$ .

4.

Determination of ${\mu }_{\text{Raman}}$ by Correction Functions

To link Raman signals from turbid samples directly to Raman scattering coefficients, we had to correct the Raman signal for the influence of the macroscopic optical properties. The measured or calculated relative Raman signal, which depended on ${\mu }_a$ , ${\mu }_s^{\prime }$ , and ${\mu }_{\mathit{\text{Raman}}}$ but also on the sample geometry, measurement geometry, and sensitivity of the setup, is called ${\mathit{Ram}}_T$ . We assumed that the application of a dimensionless correction function $F({\mu }_a,{\mu }_s^{\prime })$ , which was sample and measurement geometry-dependent, led to

4.1.

Correction Utilizing the Elastic Diffuse Reflectance Signal

It would be very convenient if a simple ratio of ${\mathit{Ram}}_T$ and $R$ could completely remove the influence of the optical properties. However, as shown for our measurement geometry in Fig. 6 , the ratio only decreases the influence of ${\mu }_a$ . The influence of ${\mu }_s^{\prime }$ increases except for the parameter range ${\mu }_a=0.01\text{to}0.1{\mathrm{mm}}^{-1}$ and ${\mu }_s^{\prime }>0.4{\mathrm{mm}}^{-1}$ , where the influence of ${\mu }_s^{\prime }$ decreases as well. This explains why the ratio improved the results in Ref. 14, where a measurement spot of similar size was used $(r_{\mathit{ex}}\approx r_{\mathit{det}}\approx 1\mathrm{mm})$ .

Fig. 6

Dependency of the ratio of simulated ${\mathit{Ram}}_T$ to elastic diffuse reflectance $R$ on the optical properties (given in units of $1∕\mathrm{mm}$ ) for our measurement geometry. (a) ${\mu }_a$ dependency: ${\mu }_s^{\prime }$ is increasing from bottom to top. (b) ${\mu }_s^{\prime }$ dependency: ${\mu }_a$ is increasing from top to bottom.

In Fig. 7 , ${\mathit{Ram}}_T$ is plotted versus $R$ for various combinations of ${\mu }_a$ and ${\mu }_s^{\prime }$ and for three different measurement geometries. Obviously, ${\mathit{Ram}}_T$ cannot be described as a function of $R$ only without knowledge of the other variables. However, it was shown recently for measurements in the NIR that ${\mathit{Ram}}_T\times {\mu }_t$ , with ${\mu }_t={\mu }_a+{\mu }_s$ , can be described by a function of $R$ only; see Refs. 12, 13. This leads to the correction function $F_1={\mu }_t\times l_1∕f_1\left(R\right)$ , where $f_1\left(R\right)$ is a geometry-dependent function of the diffuse reflectance $\left(R\right)$ , and $l_1$ is a length to make $F$ dimensionless as defined in Eq. 1. Using our Monte Carlo model, we searched for similar empirical correction functions and evaluated their validity with respect to the range of optical properties and the measurement geometry. While we found some correction functions that were only valid in the semi-infinite measurement geometry or the VIS parameter range, only $F_1$ and $F_2={\mu }_s^{\prime }\times l_2∕f_2\left(R\right)$ could be used for our measurement geometry in the whole parameter range investigated here. Since $F_1$ and $F_2$ performed very similarly, we present only the application of $F_2$ in Fig. 7. It is obvious that $F_2$ works in the semi-infinite as well as in our measurement geometry, but not for a measurement spot with a diameter of $200\mu \mathrm{m}$ . In order to find a function $f_2\left(R\right)$ for the whole parameter range, we applied a 4th-order polynomial to the log of the data. We accounted for different optical properties at the excitation and Stokes wavelength by using the mean optical properties. A wavelength-dependent variation of $R$ can be accounted for by using the geometric mean $\sqrt{R_{\mathit{ex}}{R}_s}$ instead.

Fig. 7

(a) ${\mathit{Ram}}_T$ cannot be described as a function of $R$ only without knowledge of the optical properties and the geometry. (b) ${\mathit{Ram}}_T∕{\mu }_{\mathit{\text{Raman}}}\times {\mu }_s^{\prime }$ can be described as a geometry-dependent function of $R$ only, but not for all possible measurement geometries. The optical properties used to calculated the data points are ${\mu }_a=[0.01,0.02,0.05,0.1,0.2,0.5,1,5{\mathrm{mm}}^{-1}]$ from top to bottom, and ${\mu }_s^{\prime }=[1,2,5,10{\mathrm{mm}}^{-1}]$ from left to right. The different symbols correspond to three different measurement geometries in terms of the radius of excitation $r_{\mathit{ex}}$ and detection fiber $r_{\mathit{det}}$ : (1) $r_{\mathit{ex}}=0.1\mathrm{mm}$ , $r_{\mathit{det}}=100\mathrm{mm}$ ; (2) $r_{\mathit{ex}}=1.5\mathrm{mm}$ , $r_{\mathit{det}}=0.7\mathrm{mm}$ ; (3) $r_{\mathit{ex}}=0.1\mathrm{mm}$ , $r_{\mathit{det}}=0.1\mathrm{mm}$ .

When $F_2$ was applied to the simulated data in Fig. 4, the variation of the data decreased from 109% to 6%. This means that the accuracy of the ${\mu }_{\mathit{\text{Raman}}}$ to be determined is strongly improved. When $F_2$ was applied to our silicone phantom data (Sec. 2.3), the variation of the data decreased from 88% to 21%. The remaining variability of the signal could result from variations in the sample preparation and from the fact that the diffuse reflectance was only measured at the excitation wavelength. Therefore, an experimental investigation of the correction effect was limited with this setup.

5.

Determination of Molecular Raman Cross Section or Raman Scatterer Concentration

In Sec. 4 we described how ${\mu }_{\mathit{\text{Raman}}}$ could be determined once the calibration factors were measured. Another way to determine ${\mu }_{\mathit{\text{Raman}}}$ independent from a calibration measurement is by comparing the simulated ratio of the Raman to elastic signal (see Fig. 6) to the measured ratio given in Fig. 5. We estimated the magnitude of ${\mu }_{\mathit{\text{Raman}}}$ to be $\approx 3\times {10}^{-7}{\mathrm{mm}}^{-1}$ for our phantoms with $10\mu \mathrm{g}∕\mathrm{ml}$ $(\approx 2\times {10}^{-5}\mathrm{M})$ .

For known Raman scatterer concentration c and ${\mu }_{\mathit{\text{Raman}}}$ , the differential molecular Raman scattering cross-section $\sigma$ can be calculated by

Assuming the molecular Raman scattering cross-section $\sigma$ can be found in the literature, the concentration $c$ of the Raman scatterer can be calculated by Eq. 3 using the measured ${\mu }_{\mathit{\text{Raman}}}$ .

6.

Summary and Discussion

We introduced a method to combine Raman and spatially resolved reflectance spectroscopy that allows the measurement of the optical properties together with the Raman spectrum. To determine the absolute Raman scatterer coefficients ${\mu }_{\mathit{\text{Raman}}}$ independent from the optical properties, we compared to different methods to correct for the turbidity-induced signal variation. One correction method utilized empirical functions that are dependent on the elastic diffuse reflectance but also require knowledge of the optical properties. An alternative method used a Monte Carlo model to generate and propagate Raman signals in turbid media to directly calculate the correction factors, regardless of any functional relationship.

Both correction methods reduced the turbidity-induced signal variation by 95% or more for our simulated data, and by 76% or more for our phantom data. Both correction methods need to be adjusted once to the excitation and the detection using phantom measurements or Monte Carlo simulations. Since the empirical equations still require knowledge of the optical properties, we see no advantage to this method. The correction method based on the Monte Carlo simulation, however, does not require the measurement of elastic reflectance in the Raman measurement geometry and can even be used in geometries such as a $200\text{-}\mu \mathrm{m}$ measurement spot, where an empirical correction function has not been described yet. For the data presented here, we used correction functions obtained for samples that were homogenous within the sampling volume. However, this simplification may be critical if a sampling volume includes different layers, like the epidermis and the dermis in skin measurements. The influence of the sample structure on the correction function is the subject of a current investigation.

From the determined ${\mu }_{\mathit{\text{Raman}}}$ , we calculated the molecular Raman scattering cross-section $\sigma$ of a given $\beta$ -carotene concentration in our silicone phantoms. For a known $\sigma$ , the absolute Raman scatterer concentration can be determined from ${\mu }_{\mathit{\text{Raman}}}$ . For unknown molecular Raman cross-sections, which is a possible scenario especially in complex biological matrices, the methods can also correct Raman signal variations induced by sample-to-sample turbidity variation, which allows the quantitative comparison of results obtained from different individuals, tissues, or tissue sites.

Although we took advantage of the outstanding high-resonance Raman scattering cross-section of $\beta$ -carotene, the correction methods are transferable to any other Raman scatterer. For the correction of multiple Raman peaks at different wavelengths, the optical property measurement can be easily performed using white light instead of lasers. In addition to biomedical applications, the correction method might also be beneficial in the food industry or environmental research, where quantitative measurements of analytes independent of turbidity variations are necessary. Furthermore, we suggest that the correction procedures may be applied similarly to fluorescence measurements.