|Home | About | Journals | Submit | Contact Us | Français|
We present a novel technique, intrinsic Raman spectroscopy (IRS), to correct turbidity-induced Raman spectral distortions, resulting in the intrinsic Raman spectrum that would be observed in the absence of scattering and absorption. We develop an expression relating the observed and intrinsic Raman spectra through diffuse reflectance using the photon migration depiction of light transport. Numerical simulations are employed to validate the theoretical results and study the dependence of this expression on sample size and elastic scattering anisotropy.
Sample-to-sample turbidity variations are a limiting factor in non-invasive optical techniques. Macroscopically, light propagation in turbid media such as biological tissue is governed by elastic scattering and absorption of the media. For example, prominent absorption features of hemoglobin can lead to spectral shape distortions in tissue fluorescence spectroscopy [1–5], which subsequently confound interpretation of underlying fluorophores. Similarly, in diffuse reflectance and Raman spectroscopy, turbidity variations can cause sampling volume variations across samples [6–11]. Such turbidity-induced sampling volume variations introduce additional non analyte-specific variance into subsequent data analysis. This additional variance results in increased prediction error of analyte concentrations and a calibration algorithm that is not robust.
Many researchers have developed methods to correct for spectral distortions in biological fluorescence spectroscopy in which the shape of the observed spectrum is significantly altered by the presence of absorbers such as hemoglobin [3, 12–14]. Our laboratory has utilized diffuse reflectance spectroscopy (DRS) in the development of intrinsic fluorescence spectroscopy (IFS) to correct for turbidity distortions, particularly, absorption-induced spectral distortions of the fluorescence line shape [1, 2, 5]. Diffusely reflected light is the backward emission which undergoes numerous elastic scattering events before re-emerging from the tissue, and thereby provides a metric for the amount of tissue absorption and scattering present. The optical properties of a given sample at a particular wavelength can therefore be measured in situ by monitoring the diffuse reflectance at that wavelength. Similarly, DRS can monitor optical properties at multiple wavelengths. By measuring fluorescence and diffuse reflectance at the excitation wavelength and over the fluorescence emission wavelengths using the same excitation-collection geometry, the method of IFS can be used to remove these distortions. The underlying principle is that in a turbid medium fluorescence excitation-emission undergoes similar scattering and absorption as diffuse reflectance. For IFS, the main goal is to remove spectral shape distortions, exemplified by the hemoglobin absorption peak near 420 nm.
We are developing quantitative biological Raman spectroscopy to measure analyte concentrations in the near infrared (NIR) wavelength range (~ 830–1000 nm) . In these studies, shape distortion is less of an issue because of the lack of absorbers with prominent spectral features in the NIR wavelength range . However, for quantitative analysis, turbidity-induced sampling volume variations become very significant. Consider, for example, two identical biological tissue samples containing a Raman analyte, except that the second sample has a larger scattering coefficient than the first. The increased scattering causes light to be localized in a smaller volume, with a corresponding higher efficiency for the collection lens. As a result, the size of the Raman signal in the second sample will be larger than that of the first, and the measured concentration of the Raman analyte will be different. Furthermore, since both samples have some degree of scattering, the measured Raman concentration in both samples will differ from the actual concentration.
In the Raman spectroscopy literature, some researchers have applied corrections based on direct absorption spectroscopy [17–18]. Waters extended the formalism developed by Kubelka and Munk to relate the Raman spectrum to the measured diffuse reflectance as a function of either the Kubelka-Munk absorption or scattering coefficient . This method assumes that only one optical property is changing at a time. Thus, for powdered samples, where the size of the particles and therefore their scattering characteristic change little over time, the effect of absorption from a progressively darkening sample on the Raman spectrum can be sufficiently removed [20–21]. However, the Kubelka-Munk formalism is not generally applicable to biological tissue because it assumes isotropic scattering, and biological tissue scattering is known to be anisotropic . Recently, a method of retrieving Raman spectra from subsurface layers in diffusely scattering media was implemented, in which Raman scattered light from surface regions laterally offset from the excitation laser spot was collected and analyzed via multivariate techniques . However, correction for sampling volume was not considered.
Under the photon migration framework developed by Wu et al. [1, 24], the same general principle that applies for IFS should hold true for Raman spectroscopy as well. The goal of this paper is to present a method which corrects measured analyte concentrations for turbidity-induced sampling volume variations. We refer to this method as intrinsic Raman spectroscopy (IRS). Starting from the photon migration theory, we review the analytical model for extracting intrinsic fluorescence (the fluorescence as it would be observed in the absence of scattering and absorption) and derive a parallel expression for the intrinsic Raman spectrum. Monte Carlo simulations are then employed to demonstrate its validity and elucidate the relationship between observed Raman and diffuse reflectance for semi-infinite samples and samples of finite dimension. An analysis of the methodology with respect to sample size and scattering anisotropy is presented. An experimental study is presented in the companion paper to validate this method and demonstrate that it can be used experimentally to correct for turbidity-induced spectral distortions and sampling volume variations.
Most analytical and numerical models of light propagation in tissue employ macroscopic optical properties, including the absorption coefficient, μa (cm−1), the scattering coefficient, μs (cm−1), and the elastic scattering anisotropy, g = <cosθ>, the average cosine of the single scattering angle, θ. The absorption and scattering coefficients are the probability of a photon being absorbed or scattered, respectively, per unit path length. The sum of μa and μs is called the total attenuation coefficient, μt, with its inverse defined as the mean free path. As mentioned in the Introduction, turbidity, caused by variations in μa and μs, can introduce sampling volume variations which prevent accurate concentration measurements.
Light propagation in turbid media can be described by the radiative transfer equation . However, the analytical solution to this integro-differential equation can be found only for very special conditions, with approximations. Diffusion theory is one of the most extensively studied approximations. The diffusion equation, along with appropriate boundary conditions dictated by the sample geometry, may be solved to provide the fluence distribution inside the sample and the reflected flux. Diffusion theory is often used to model photons that experience multiple scattering events and thus propagate “diffusively”, usually with a certain amount of source-detector separation .
Photon migration theory, developed by Wu et al.  employs probabilistic concepts to describe the scattering of light and to set up a framework that allows the calculation of the diffuse reflectance from semi-infinite turbid media. In this framework the diffuse reflectance from a semi-infinite medium can be written as:
with fn(g) the photon escape probability distribution, n the number of scattering events before escaping, g the scattering anisotropy, and a the albedo, μs/(μs+μ a). Two fundamental assumptions are made in photon migration theory: (1) The photon escape probability distribution of a semi-infinite medium depends only on the number of scattering events and the anisotropy; and (2) the escape probability distribution can be approximated by an exponential function of the form fn(g)= k(g)e−k(g)n. Using this, Eq. (1) can be further approximated by an integral form:
with k(g) an anisotropy and geometry dependent parameter that can be expressed as the product S(1-g) where S depends on the excitation-collection geometry. This quantity must be calibrated for each experiment. As demonstrated by Zonios et al. , this expression models the experimental results well.
Using photon migration theory, Wu et al.  derived an analytical expression relating the observed fluorescence (FOBS) and the diffuse reflectance (Rd) to the intrinsic fluorescence (FINT), defined as the fluorescence measured from an optically-thin slice of tissue, free of distortions due to scattering and absorption. This expression and its variants have been employed to recover turbidity-free fluorescence spectra from various types of tissue. The resulting intrinsic fluorescence facilitates interpretation of the underlying fluorophores and consequently improves the accuracy of disease diagnosis .
Following Wu’s derivation, the observed Raman spectrum (RamOBS) at various turbidities can be written as:
with μsR the Raman scattering coefficient, μt the total attenuation coefficient at the excitation wavelength, and RamINT the turbidity-corrected “intrinsic” Raman spectrum. RamINT is a dimensionless quantity equal to μsR/l, where l is a characteristic length. Subscripts ‘x’ and ‘R’ denote quantities at the excitation and Raman wavelengths, respectively, and the other symbols are defined above. The intrinsic Raman spectrum under semi-infinite conditions is therefore:
From this expression, the intrinsic Raman spectrum must be related to the observed Raman spectrum through a calibration procedure. In intrinsic fluorescence spectroscopy the calibration is applied to each term, k(g), Rd,x, and Rd,R, individually, through a combination of modeling and experimental data . This procedure is not feasible for intrinsic Raman spectroscopy in the case of interest here, where both Raman and diffuse reflectance are measured in the region of NIR, because the absorption features are exceedingly weak and indistinct and the individual parameters are difficult to specify. Therefore, an alternative procedure must be found. One approach is to approximate Eq. (5) by an expression proportional to a function of the ratio of RamOBS to Rd. This was explored in reference , and will be the subject of a forthcoming publication. An alternative, hybrid approach is taken here. Because k(g), Rd,x, and Rd,R, are all related to μs, μa and g, we can combine k(g)/(Rd,xRd,R) into a function of the easily measurable diffuse reflectance at the Raman wavelength (Rd,R, or simply Rd hereafter),
We call f(Rd) the calibration function. It includes the factor l and has dimensions of inverse length, and can be obtained from a set of tissue phantoms with varying μs and μa. Note that RamOBS and Rd must be measured in the same excitation-collection geometry. An additional advantage of this approach is that instrument-dependent parameters such as the excitation-collection geometry can be calibrated in one step. Such a calibration function can be applied to future Raman spectra provided μt and Rd have been determined.
In the following, Monte Carlo simulations are employed to elucidate the relationship between the observed Raman and the intrinsic Raman spectra. We also study the dependence of this calibration function on sample size and scattering anisotropy.
The Monte Carlo code employed here is based on an existing open source code developed by Jacques  for diffuse reflectance and fluorescence in a single layer medium. The code has two steps: In the first step, photons are injected from the top surface of the sample and propagate in the presence of elastic scattering and absorption. The absorbed quantity at each bin is recorded and used as the initial photon weight for the second step, in which fluorescence photons are launched from each bin isotropically and then propagated. The outcome includes diffuse reflectance, from the first step, and fluorescence, from the second step. A few modifications were made in our code: (1) A fixed-weight scheme was employed for photon weight bookkeeping, i.e., the weight (intensity) of a photon stays the same as long as it is not absorbed or Raman scattered; when absorption or Raman scattering occurs, the photon weight is reduced to zero. This approach is similar to that of Welch et al.  and gives comparable results to Jacques’ albedo-based approach, but is perhaps more representative of a natural process. (2) Secondary Raman scattering is neglected, i.e., a Raman photon cannot generate another Raman photon. This is a good assumption, since Raman scattering is a very rare event. (3) Angular resolution for the exiting flux was added and we provide for finite sample size. This allows us to compare results from different numerical apertures as well as various sample sizes. Cylindrical coordinates have been employed in the Monte Carlo code, thus the sample geometry reported in the format of radius by depth.
The simulation begins with injection of a photon into the medium with a calculated step size sampled from the probability distribution px=μsexp(-μsx) (Note that the average step size is 1/μ s). A probability check (ξ<Pa; Pa=μ a/μt) determines whether the photon is absorbed or scattered. If the photon is absorbed, it is terminated, and the location of absorption recorded. If it is not absorbed, it is scattered (with probability Ps=1-Pa). If the scattered photon passes a Raman probability check (ξ<PR; PR=μsR/(μs+μsR)), a new photon is launched isotropically at the Raman wavelength and propagates until it exits or is absorbed; otherwise the scattered photon continues to propagate without wavelength shift. Subsequent scattering angles are determined using a Heney-Greenstein phase function, a good approximation for biological tissue. Diffuse reflectance consists of the collected photons at its original wavelength, and the observed Raman spectrum consists of the collected photons with wavelength shift.
A series of Monte Carlo simulations were performed. The goal was to study the relationship between Raman scattering and diffuse reflectance under various turbidity conditions. We simulated 49 samples, following a 7×7 matrix of scattering and absorption properties with ranges similar to those found in biological tissue , in which scattering usually dominates absorption. The scattering coefficient, μs, was varied from 18.4 to 99.4 cm−1 (18.4, 36.8, 50.6, 62.6, 73.6, 87.4, 99.4) and the absorption coefficient, μa, was varied from 0.1 to 1.4 cm−1 (0.1, 0.15, 0.2, 0.36, 0.5, 0.95, 1.4). The scattering anisotropy, g, was kept at 0.8 unless mentioned otherwise. A Raman scatterer of constant strength was present in each sample to serve as an indicator of the Raman spectrum, i.e., RamINT in Eq. (6) was kept constant throughout the data set. The excitation beam was collimated with 0.1 cm radius for all simulations. The collection angle was between ±45 degrees; the collection spot was concentric with the excitation spot but with various radii, specified for each simulation in the Results and discussion section.
Using the Monte Carlo code described above, we studied the relationship between the observed and the intrinsic Raman spectra in two scenarios: semi-infinite and non-semi-infinite sample geometries, at a single Raman wavelength. From the simulation results, qualitatively we observe that Raman intensity changes with the diffuse reflectance, supporting the basic principle that both of them experience similar turbidity-induced distortions.
The one-to-one relationship between the observed fluorescence intensity and the diffuse reflectance intensity has been demonstrated  with either μs or μa fixed while the other varies. Fig. 1 shows that the one-to-one relationship does not hold when both optical properties are allowed to vary simultaneously, a situation where extra dependence on μt has to be considered as can be seen in Eq. (6). When μt is factored in, a new one-to-one relationship is revealed (Fig. 2). A fit to the curve provides the calibration function, f(Rd), which then provides the IRS correction when both optical properties are allowed to vary. Note that (RamOBSμt), RamOBS, and Rd have been normalized to their respective maximum values in the 7×7 sample set. The functional form used for fitting is not important; the best fit to the data may be used. We find that an exponential fit may best approximate the observed function and is employed in the companion experimental paper. However, for illustrative purposes in this paper, we find it useful to fit with a power law as the power law exponent is a convenient metric for curvature changes. It is interesting to note that RamOBS appears to be only weakly dependent on Rd, particularly for larger μa values, as shown in Fig. 1. This implies that RamOBS is proportional to 1/μa. Thus, if one substitutes 1/μa for RamOBS in Fig. 2, the plot is approximately an inverse Rd vs. μt/μa. This interpretation provides a heuristic explanation to why the data points in Fig. 2 collapse onto a smooth curve.
Under the semi-infinite condition with all return light collected, the diffuse reflectance observed with various turbidities is well described by a function of the ratio μs/μa according to scale invariance [26, 31–32]. We next consider non semi-infinite sample geometries, and vary the turbidity and study the diffuse reflectance with respect to μs/μa. Three sample geometries were simulated with various collection spot radii and fixed delivery spot radius (0.1 cm) and g (0.8). We observe a general trend in the results shown in Fig. 3–Fig. 5 that the diffuse reflectance is not simply a function of the ratio μs/μa for arbitrary size samples and collection spot radii. However, it approaches a function of only the ratio μs/μa when the sample or the collection spot radii becomes sufficiently large. This demonstrates that the above-mentioned scale invariance, i.e., the ideal case in which Rd=f(μs/μa), is only preserved when both the sample and the collection radius are semi-infinite. We find that the sample has to be larger than 8δ (the penetration depth, δ = [3μa(μa+μs′)]−−0.5) in any dimension to be considered truly semi-infinite without any boundary effects, in agreement with reference . The average penetration depth in our simulated sample design is 0.33 cm, and thus the 8δ criterion is often violated, however, we observe that the diffuse reflectance can be well approximated by a one-parameter (μs/μa) function using a 2 cm (r) × 2cm (z) sample with a 2-cm collection spot radius, but not with a small collection spot radius such as 0.4 cm. Given that the collection spot radius is ~0.3 cm in our instrument, we do not expect scale invariance to hold for the diffuse reflectance.
It is important to note that although diffuse reflectance versus μs/μa deviates significantly from the semi-infinite behavior, Eq. (6) is still applicable. In other words, the variability introduced by sample geometry and collection spot size can be captured by both the diffuse reflectance and the Raman spectrum, and thus the relationship still holds. Fig. 6 shows the simulation results of calibration function for 4 different sample sizes. Comparing with Fig. 2, we observe that the curvature, as indicated by the exponent in the power law fit, increases as the sample size increases. The larger data spread for smaller sample size is likely due to higher photon loss through the boundaries, which reduces the signal-to-noise ratio. It may also be due to a decrease of variability captured by the calibration function with smaller sample size, as the semi-infinite condition was assumed in deriving Eq. (6). However, the curve still appears to be a good, unbiased fit to the data.
From Monte Carlo simulations we have learned that the curvature of the calibration function increases with the sample size. An analogous phenomenon can be observed in Fig. 7 when the anisotropy (g) is varied from 0.99 to 0.7. In the high anisotropy (g=0.99) case, the scattered light is highly forward directed, resulting in an effective path length much longer than in the low anisotropy case; the exponent appears smaller, as in the smaller (i.e., non-semi-infinite) samples. Further, the fact that the relationship between RamOBSμt and Rd becomes close to linear suggests that the Raman scattering and diffuse reflectance become more linearly correlated in effective path length. This is certainly the case when light scattering is more forward directed.
The effect of sample size and scattering anisotropy on the curvature of the calibration function can be studied collectively using the Monte Carlo results shown in Fig. 8. Consideration of these two parameters and the interplay between them is important for experimental design. For example, it is known that whole blood is highly forward scattering with μs > 300 cm−1 and g ~0.99 . This implies that the exponent will be close to 1 and therefore the calibration curve will be more linear.
To apply the intrinsic Raman spectroscopy correction, one needs to know the total attenuation coefficient μt. In Eq. (6), when μs >> μa and g is relatively constant, μt can be conveniently replaced by μs′, the predominantly obtained quantity in the literature [32, 35–38] based on diffusion theory or variants of it. Our laboratory routinely measures optical properties from biological tissue in other wavelength ranges , and a similar method can be employed for this purpose.
Turbidity variation is one of the major limitations in non-invasive quantitative biological Raman spectroscopy. To overcome this limitation, we have developed the theory of intrinsic Raman spectroscopy and provided validation using numerical simulations. We employed the photon migration theory for light propagation in turbid media, and related the intrinsic and the observed Raman spectra through diffuse reflectance. We employed Monte Carlo simulation to validate the modeling results and study the influence of sample geometry and scattering anisotropy on the curvature between the intrinsic and the observed Raman spectra. These results provide a systematic way to correct turbidity-induced sampling volume variations in NIR quantitative biological Raman spectroscopy. An immediate benefit of intrinsic Raman spectroscopy is the improvement in non-invasive determination of analyte concentrations in turbid media. The companion paper presents an experimental study to validate this method and demonstrate that it can be used experimentally to correct for turbidity-induced sampling volume variations.
This work was performed at the MIT Laser Biomedical Research Center and supported by the NIH National Center for Research Resources, grant P41-RR02594, and a grant from Bayer Health Care, LLC. Wei-Chuan Shih was a MIT Martin Fellow for Sustainability. We thank Ishan Barman, Gajendra Singh, and Ramachandra Dasari for valuable discussions.
OCIS codes: (170.1470) Blood or tissue constituent monitoring; (170.3660) Light propagation in tissue; (170.5280) Photon migration; (170.5660) Raman spectroscopy; (170.7050) Turbid media.