2.1.

Experimental Studies on Human Subjects

The concentration datasets as well as the tissue fluorescence and photobleaching profile used in our numerical simulations (Sec. 3) was originally described in one of our laboratory's previous publications.⁴ Briefly, NIR Raman spectra were acquired from the forearms of 20 human volunteers undergoing oral glucose tolerance tests (OGTT). An 830-nm diode laser (Process Instruments) was used for Raman excitation, and the back-scattered light from the tissue was collected on a liquid-nitrogen cooled charge-coupled device (CCD) (Roper Scientific, Ottobrunn, Germany) coupled to a f/1.8i spectrograph (Kaiser Optical Systems, Ann Arbor, Michigan) For each volunteer, spectral measurements were initiated following the ingestion of a 220-ml SUN-DEX solution (Fisher Hamilton LLC, Two Rivers, Wisconsin) containing 75 g of glucose. Spectra were acquired for a total of 3 min per data point at 5-min intervals over approximately a two-hour time period, forming a measurement time series for each volunteer. A representative spectral time series acquired from a volunteer is shown in Fig. 1a. The second panel, Fig. 1b, gives the photobleaching profile as measured on this volunteer. It is a curve of I _λ/I _exc as a function of time, where I _λ is the area under the spectral curve over a 20-cm^–1 band surrounding any arbitrarily selected central value (600 cm^–1 used for data shown here), and I _exc is the measured intensity of the excitation source. The measured values of I _λ/I _exc were fitted with a double exponential curve,²¹ as shown in the figure.

Fig. 1

(a) Raman spectra acquired from a human volunteer during an OGTT. (b) Normalized tissue autofluorescence decay (I _λ/I _exc) obtained as a function of time during the OGTT performance (the measured data points are in circles and the dotted line is the best-fit double exponential curve).

Concomitant with the spectral measurements, blood samples were also collected using finger sticks at regular 10-min intervals for estimation of glucose levels by a clinical glucose analyzer (HemoCue, Lake Forest, California). To correlate with the spectral measurements, the reference values were approximated at intermediate time points using spline interpolation. For comparison, two such profiles taken from two different volunteers are shown in Fig. 2. Profile 1 shows a rise in blood glucose concentration (to almost double its normal value), followed by a return to original levels, indicative of proper functionality of the blood glucose regulation system in the volunteer. Profile 2, in sharp contrast, exhibits a nearly monotonic rise in glucose levels over the test period, characteristic of a weak insulin-mediated response in this subject.

Fig. 2

Representative glucose concentration profiles taken from two human volunteers: profile 1 (top) and profile 2 (bottom).

2.2.

Experimental Studies on Tissue Phantoms

Physical tissue models (tissue phantoms) in which the fluorescence and Raman scattering could be precisely varied were used to investigate the prediction performance of calibration models based on the correlated and uncorrelated spectral datasets, respectively. Additionally, in conjunction with the numerical simulations performed in Sec. 3, this study was used to understand the effect of fluorescence removal strategies on the prospective prediction accuracy of the developed calibration models.

A total of 50 tissue phantoms were prepared with aqueous solutions of glucose (analyte of interest), urea (spectral interferent), and indocyanine green (ICG). ICG was employed to produce a fluorescence background similar to the tissue autofluorescence observed with near-infrared (NIR) excitation. To mimic the (photobleached) fluorescence correlation observed in glucose tolerance testing of human subjects with impaired glucose tolerance (see Profile 2 in Fig. 2), 20 tissue phantoms were constructed, where the correlation between the glucose and ICG concentrations was high (R ² ≈ 0.95). The other 30 tissue phantoms had negligible correlation between the glucose and ICG levels (R ² ≈ 0.1). For both correlated and uncorrelated tissue phantoms, glucose and ICG concentrations varied from 0.9 to 1.2 M and 0.18 to 0.36 μM, respectively. The relatively high concentration of glucose in the tissue phantoms enabled higher fluorescence values to be studied while retaining a satisfactory signal-to-noise ratio. The urea concentration in all the tissue phantoms was randomized in the range of 250 to 400 mM.

For Raman excitation, an 830-nm diode laser was focused into the delivery fiber (200-μm core diameter) of an optical fiber and filter probe assembly²² to deliver an average power of ∼100 mW on a spot size of ∼1 mm² at the sample. The sample solution was held in a cylindrical glass vial (one inch diameter) with a top opening, through which the probe tip was immersed directly into the solution without contacting the container (Fig. 3). The back-scattered Raman light was collected by ten collection fibers (200-μm core diameter) before being dispersed via a tunable f/2.0 spectrograph (LS-785, 5 cm^–1 resolution, Princeton Instruments, Trenton, New Jersey) onto a TE-cooled CCD detector (PIXIS 256, Princeton Instruments). The acquisition time for each Raman spectrum was 20 sec. Spectra of the phantoms were acquired randomly with respect to the constituents' concentration to eliminate any potential temporal correlation. Since the optical fiber probe itself generated prominent Raman spectral features at wavelengths below 800 cm^–1, the spectral range of 950 to 1800 cm^–1 was used for data analysis (PLS calibration and prediction) after vertical binning and cosmic ray removal. Additionally, to investigate strategies for removing the effects of the fluorescence background, SSRS was performed by acquiring the tissue phantom spectra after tuning the spectrograph grating position by ∼25 cm^–1 (approximately the full width half maximum of the characteristic glucose peak at 1100 cm^–1). The Raman shift axis was calibrated using the standard spectra measured from a 4-acetaminophenol powder sample. The tuning of the grating was monitored by the position of the laser line on the CCD detector such that an 18-pixel shift was maintained between the two positions. Further details of the calibration procedure followed for SSRS implementation are described elsewhere.²³ The difference between the pair of Raman spectra (acquired at two different grating positions) was used for PLS calibration and prediction, and the performance was compared with that of the PLS calibration model constructed using the raw Raman spectra. We have also employed first derivative spectra for comparison with the aforementioned calibration models. For the first derivative spectra, discrete differencing of the acquired spectra was performed. Additionally, we undertook a derivative application after smoothing of the acquired spectra (Savitzky and Golay,²⁴ nine-point window). All mathematical processing and numerical simulations were performed using in-house software codes written in Matlab (The Mathworks, Natick, Massachusetts).

Fig. 3

A schematic diagram of the experimental setup. Raman spectra were obtained from tissue phantom solutions using an optical fiber probe, which included a laser light delivery fiber and ten collection fibers. The spectrograph was equipped with a micrometer, which was used to precisely tune the grating for implementing SSRS. F1: laser line filter. S: shutter. L1: focusing lens for optical fiber coupling. F2: Rayleigh rejection edge filter.

3.1.

Effect of Photobleaching Study

Simulation calibration spectra were first generated by forming weighted linear combinations of constituent analyte spectra of glucose, creatinine, and urea, as measured by our Raman system. For our application, glucose represents the analyte of interest, with the other two constituents playing the role of spectral interferents in the mixture model. The assigned weights in the linear combinations correspond to the concentrations of the analytes in the sample at that particular time instance. Two sets of simulation spectra were constructed: a correlated spectral dataset and an uncorrelated spectral dataset. Specifically, the correlation coefficient ρ _A,B between two random variables A and B is given by:

To generate the correlated spectral dataset, the concentrations shown in profile 2 (Fig. 2) were assigned to glucose to formulate the mixture Raman spectra. The glucose concentration values of profile 2 show a strong negative correlation with the photobleached fluorescence levels plotted in Fig. 1b (ρ _A,B ≈ 0.95, R ² ≈ 0.9). For the uncorrelated dataset, the concentration values from profile 1 (Fig. 2) were allocated for glucose, as the correlation between this profile and the aforementioned photobleaching curve is insignificant (R ² ≈ 0.1). It should be noted that the concentrations of creatinine and urea in both sets of simulation spectra were kept nearly constant (i.e., within 5% variation) to replicate the negligible changes observed in the concentration levels of the other analytes during a typical OGTT.

To model the fluorescence, we performed hypersmoothing of a typical tissue spectrum to remove the sharper Raman features and retain only the autofluorescence background. For both correlated and uncorrelated dataset simulations, this hypersmoothed spectrum was scaled by a factor specified by the photobleaching profile [Fig. 1b] at that particular time before addition to the mixture of Raman spectra. Moreover, to fully evaluate the impact of fluorescence-to-Raman amplitude ratio on the performance of the calibration models for both the correlated and uncorrelated cases, the relative intensity of the fluorescence and Raman signals was also varied in the range of 5 to 20. Finally, noise was generated in accordance with the shot-noise-limited condition (which has previously been observed to hold true for spectra acquired using our Raman system²⁶) for addition to the fluorescence-Raman combination spectra. We revisit the shot-noise limited assumption and its implications in Sec. 4.2. The signal-to-noise ratio (SNR), defined here as the ratio of the maximum intensity of the fluorescence-Raman combination signal to the mean noise magnitude, was varied from 50 to 200. The range of fluorescence-to-Raman amplitudes and SNR values mentioned before are consistent with that observed in Raman spectra acquired for these types of analytical measurements (typical values of observed SNR ∼ 100 and fluorescence-to-Raman ratio ∼ 10). At each setting of the tunable parameters, i.e., fluorescence-to-Raman ratio and SNR, 100 calibration spectra were generated for the correlated as well as the uncorrelated datasets.

The prediction spectra were generated following a scheme similar to that stated before, except that the concentrations of all the analytes were completely randomized in this dataset. The fluorescence levels in the prediction spectra were also completely randomized to simulate those observed in a set of unknown samples (such as a test set of human volunteers), where the fluorescence levels and the concentrations of the analytes would be expected to possess minimal correlation with one another. This situation typically occurs when a calibration model developed during a tolerance or loading test is used prospectively on different subjects to determine their glucose levels. For a given setting of the fluorescence-to-Raman ratio and SNR level, 50 such spectra were included in the prediction dataset.

Partial least squares (PLS) analysis was used to develop the calibration models based on leave-one-out cross-validation on the calibration (training) set. The reported final errors [root-mean-square error of prediction (RMSEP), Eq. 2] were obtained by employing these calibration models on the independent prediction sets.

3.2.

Fluorescence Removal Study

In this simulation study, we followed an algorithm similar to that outlined in Sec. 3.1, except that all the spectra (i.e., calibration and prediction, both for correlated and uncorrelated cases) were processed for fluorescence removal. Here, we employed SSRS to quantify its benefits, relative to that of numerical postprocessing schemes, for both the photobleaching correlated and uncorrelated cases.

For our simulations, we generated a second set of calibration and prediction spectra for both the correlated and uncorrelated datasets by shifting the composite mixture spectra (Raman plus fluorescence) by 25 cm^–1. These shifted spectra were then subtracted from the original set of spectra to obtain the Raman difference spectra (plus shot noise and residual fluorescence background). While most of the previous research has primarily been focused on reconstruction of the intrinsic Raman spectrum from the difference spectra,^{27, 28, 29} we have employed the difference spectra directly in analysis to extract concentration information through PLS. The direct incorporation of the difference spectra into calibration analysis is advantageous, as it avoids the artificial spectral features and noise (ringing behavior) that may appear due to the application of a reconstruction method (such as Fourier deconvolution).

For fair comparison with numerical strategies, we implement three variations of least- squares polynomial subtraction as alternative fluorescence removal tools. First, a lower order polynomial is fit to the broad fluorescence background and subtracted. Specifically, we employ a second or third order polynomial for the PLS analysis, depending on which order provides the least error in cross-validation. This approach is widely pursued for signals with broad fluorescence backgrounds, because higher order polynomials tend to overfit the data, as also noted by Beier and Berger.³⁰ Second, an algorithm developed by Lieber and Mahadevan-Jansen,¹⁷ which automatically fits a modified polynomial just below the original spectrum, is used. Third, a recently developed adaptive “minmax” scheme³¹ is used to subtract the background fluorescence. PLS analysis is undertaken in conjunction with each of the three numerical processing schemes.

4.1.

Numerical Simulations

4.1.1.

Effect of photobleaching study

As detailed in Sec. 3.1, PLS was first used in a leave-one-out cross-validation routine on the calibration spectral set. The number of loading factors used to build the final model for prediction was in accordance with the number of constituents that provided the minimum error in cross-validation. The results of the simulations at varying levels of Raman-to-noise ratio for the photobleaching correlated and uncorrelated cases are plotted in Fig. 4, with error bars indicating the standard deviation of the RMSEP values at each setting. It is evident that the calibration model developed on the correlated calibration spectra always performs worse than the model developed on the uncorrelated set, irrespective of the SNR value or the fluorescence-to-Raman amplitude ratio. For example, at a setting of fluorescence-to-Raman of 10 and SNR of 100 (values observed typically for tissue Raman spectra), the mean prediction errors were measured to be 16.3 and 5.3 mM for the correlated and uncorrelated cases, respectively, thereby demonstrating an approximately three-fold increase in the prediction error. Expectedly, for lower levels of Raman-to-noise ratio, the prediction performance of both models drops considerably. (It is noted that changing the fluorescence-to-Raman ratio or the SNR has an identical effect, as they both modulate the Raman-to-noise ratio.) Importantly, the performance of the calibration model built on the correlated dataset drops off much more significantly than that developed on the uncorrelated dataset—as evidenced by the widening gap between the bars at higher fluorescence-to-Raman (i.e., lower Raman-to-noise) settings. This can be attributed to the fact that it becomes more difficult to distinguish the Raman spectrum of glucose from the confounding fluorescence background in a noisier spectral set. The simplest experimental solution to reducing the disparity between the correlated and uncorrelated datasets is to increase the overall SNR, which would then decrease the RMSEP; however, this would come at the cost of either increased exposure time or increased laser power.

Fig. 4

Bar plot of RMSEP obtained for calibration models developed on photobleaching correlated and uncorrelated datasets, respectively, as a function of increasing fluorescence-to-Raman ratio (i.e., decreasing Raman-to-noise ratio). Here, the SNR is held constant (100). Identical results are obtained by changing the SNR while holding the fluorescence-to-Raman ratio fixed.

The comparison of PLS prediction performance between correlated and uncorrelated datasets (Fig. 4) can be further generalized to show that the amount of deterioration in the prediction performance is dependent on the magnitude of correlation (R ²) between glucose concentration and tissue fluorescence. In Fig. 5, RMSEP increases with increasing R ² for all values of fluorescence-to-Raman ratios. At low values of R ², the RMSEP appears to be independent of R ², but displays a steep rise after a threshold value of R ². This suggests that the prediction error at low values of R ² is determined primarily by the spectroscopic SNR, while at higher values of R ², the prediction error depends on both SNR and R ². This can also be visualized by observing the widening gap between any pair of RMSEP curves with an increase in R ² in Fig. 5.

Fig. 5

RMSEP of simulations as a function of the correlation (R ²) between glucose concentration and fluorescence intensity. The fluorescence-to-Raman ratio (F/R) was varied from 10, 15 to 20, and 20 simulations were performed for each case.

The results of Figs. 4 and 5 also have extensive clinical implications, particularly for the correlated case. Such a situation is representative of experimental observations from OGTTs performed on diabetic and insulin-resistant patients. Since it is difficult to perform glucose/insulin clamping studies in normal, and particularly, diabetic populations, tolerance tests have remained the most viable protocol for the development of a calibration model using spectroscopic techniques.^{32, 33} Moreover, while test-strip-based measurements and glycosylated hemoglobin (HbA1c) examination are more extensively used at the present time, tolerance tests remain the gold standard of clinical tests of insulin/glucose control systems, and have been used on a regular basis to screen for gestational diabetes^{34, 35} and to determine the level of insulin resistance.³⁶ However, our results indicate that photobleaching introduces a systematic error in Raman spectroscopic calibration model development using tolerance tests, particularly for people with impaired glucose regulation. Put another way, any model developed on a diabetic population is susceptible to substantial errors in prospective prediction, even when used for diabetic patients only. This is clearly undesirable. The potential remedies to break the unwanted correlation include a change in experimental study design (e.g., point-of-care measurements across a large population) and removal of the fluorescence background. The inclusion of a large number of human volunteers in the calibration set may, however, introduce additional complexities such as large variations in tissue optical properties and skin heterogeneity.

4.1.2.

Fluorescence removal study

For this study, the fluorescence removed spectra—SSRS-processed and least squares polynomial subtracted—were used in conjunction with the concentration datasets to build the calibration models (as detailed in Sec. 3.2). The results of prospective prediction of the five models (i.e., based on raw, lower order polynomial subtracted, modified polynomial subtracted, minmax fit subtracted, and SSRS processed spectra) are summarized in Fig. 6. For this set of simulations, experimentally observed values of SNR (of 100) and fluorescence-to-Raman ratio (of 10) were used.

Fig. 6

Bar plot of RMSEP values obtained for glucose concentrations from simulations on photobleaching correlated (red) and uncorrelated (blue) datasets. The groups represent calibration and prediction performed using the following types of spectra: (from left to right) unprocessed, lower order polynomial subtracted, modified polynomial subtracted, minmax fit subtracted, and SSRS processed spectra, respectively.

For the correlated case, SSRS provides substantial reduction of prediction error (from 16.7 to 3.8 mg/dl), even lower than the levels obtained using the raw uncorrelated calibration model (5.5 mg/dl). Conventional lower-order polynomial subtraction, on the other hand, provides only about half the amount of improvement (from 16.7 to 10.3 mg/dl) as SSRS processing for the correlated case. This disparity in results could be attributed primarily to the introduction of nonphysical negative regions in the spectra and associated artifacts due to the process of fluorescence removal using the polynomial subtraction technique. Newer variants of the polynomial subtraction technique, namely the modified polynomial (8.3 mg/dl) and adaptive minmax fitting (6.7 mg/dl), provide lower average RMSEP values than conventional polynomial subtraction. Clearly, the inclusion of meaningful constraints in these techniques enhances the prospective prediction accuracy, though not to the levels obtained using SSRS.

Our results for the uncorrelated case suggest that application of SSRS or numerical processing do not cause similar changes in prediction performance from the raw spectra results. This is not unexpected, because implicit calibration schemes (PLS) do a reasonable job of separating the analyte of interest (e.g., glucose) from the confounding analytes (including tissue fluorescence), especially when reasonable levels of SNR are achieved.

Finally, it is important to understand that our numerical simulations have been based on the assumption that the main factor limiting SNR is the shot noise in the system. Although this assumption holds good for experiments on tissue samples where exposure times are limited by other constraints, this may not always be the case. As has been pointed out by Bell, Bourguignon, and Dennis,^{23, 37} and O’Grady,³⁸ there are many situations where it is observed that the random but fixed variations in pixel-to-pixel response on the CCD might actually be more dominant than shot noise. This is especially true when high signal levels can be achieved and exposure times can be increased without any significant downside (such as when chemical mixtures and powder samples are tested). SSRS will clearly provide a bigger improvement in such applications, because the major cause of noise arising from the irregularity in detector response can be canceled when the shifted spectra are subtracted.

4.2.

Experimental Studies on Tissue Phantoms

To validate our simulation results, we analyzed the spectral datasets obtained in our tissue phantom studies. Figure 7a shows a raw Raman spectrum obtained from a representative tissue phantom, Fig. 7b shows SSRS data obtained by subtracting two raw spectra obtained at slightly different spectrograph grating positions 25 cm^–1apart, and Fig. 7c shows the numerical first derivative of the raw spectrum shown in Fig. 7a. It is clear that the fluorescence background can be suppressed using either shifted subtraction at two grating positions (SSRS) or taking the first derivative. However, the SNR of the SSRS data appears to be significantly better than that of the first derivative spectra. This can be attributed to the cancellation of the fixed pattern noise on the CCD detector, as discussed before. It should be noted that having the correct amount of shift between the two acquired spectra is critical to the application of SSRS, as too small a shift would result in poor SNR and too large a shift would cause improper background suppression. For the SSRS spectrum, we also observe the presence of a small residual background because of the corresponding shift in fluorescence in the two sets of spectra.

Fig. 7

(a) Representative Raman spectrum acquired from a tissue phantom. (b) SSRS spectrum obtained by subtracting two spectra, obtained at spectrograph grating positions 25 cm^–1 apart. (c) First derivative spectrum.

We first created the prediction dataset by randomly extracting ten uncorrelated tissue phantom data. PLS calibration models were then constructed based on the raw spectra, SSRS spectra, and first derivative spectra for both the fluorescence-correlated and uncorrelated datasets (which consisted of 20 tissue phantoms each). Model performance during optimization was judged by comparing the standard error for leave-one-out cross-validation in the calibration dataset. The prediction set was withheld from all optimizations and was solely employed to evaluate the performance of the final optimized models. The whole procedure was iterated 50 times to obtain an average and standard deviation of the RMSEP values. Table 1 summarizes the results of the prospective prediction.

Table 1

Summary of mean and standard deviation (in parentheses) of the RMSEP values obtained from the tissue phantom experiments for both the correlated and uncorrelated datasets.

We observe that the prediction performance of the calibration model built on the fluorescence-correlated dataset is substantially worse than that built on uncorrelated phantoms. The mean prospective prediction errors are 53.8 and 39.1 mM for the correlated and uncorrelated calibration models, respectively. Importantly, we find that the application of SSRS reduces the error in the correlated calibration model from 53.8 to 42.7 mM to the levels achieved with the uncorrelated dataset. We also observe that the application of SSRS marginally increases the prediction error of the uncorrelated dataset from 39.1 to 42.4 mM. This may be attributed to the imperfect rejection of the fluorescence background by SSRS application. Nevertheless, the two primary results observed in our previous numerical simulations are reinforced by this tissue phantom study: 1. calibrations developed from the fluorescence-correlated datasets have a greater prediction error than calibrations developed with fluorescence-uncorrelated datasets, and 2. the application of frequency-shifted spectra provides more accurate concentration predictions than those obtained from the unprocessed spectra, particularly for the fluorescence-correlated dataset.

In sharp contrast to SSRS, the first derivative application appears to adversely impact the prediction performance for both the fluorescence-correlated and uncorrelated datasets. One would expect that this is an outcome of the reduction in SNR of the first derivative spectra shown in Fig. 7c. This is consistent with previous observations, especially when an additional filtering procedure is not employed.³⁹ It has also been previously reported that the derivative method makes the lower intensity Raman bands difficult to identify in the processed spectrum.¹⁷ Smoothing of the acquired spectra before the differentiation procedure, however, reduces the average RMSEP value, probably due to the resultant increase in SNR¹⁶ compared to simple derivative application.

Taken in conjunction with the least squares polynomial subtraction results from our simulations, these results indicate that while the numerical fluorescence removal strategies are able to suppress the fluorescence background, the instrumentation-based SSRS technique provides better quantitative predictions. Furthermore, the tissue phantom studies verify the robustness of the SSRS methodology in a realistic setting.