Constrained regularization for noninvasive glucose sensing using Raman spectroscopy

Multivariate calibration is an important tool for spectroscopic measurement of analyte concentrations. We present a detailed study of a hybrid multivariate calibration technique, constrained regularization (CR), and demonstrate its utility in noninvasive glucose sensing using Raman spectroscopy. Similar to partial least squares (PLS) and principal component regression (PCR), CR builds an implicit model and requires knowledge only of the concentrations of the analyte of interest. Calibration is treated as an inverse problem in which an optimal balance between model complexity and noise rejection is achieved. Prior information is included in the form of a spectroscopic constraint that can be obtained conveniently. When used with an appropriate constraint, CR provides a better calibration model compared to PLS in both numerical and experimental studies.


Introduction
Multivariate calibration is a powerful analytical technique for extracting analyte concentrations in complex chemical systems with many hidden factors. 1 Numerous applications of multivariate calibration can be found in analyzing spectroscopic data, where information about all the analytes can be collected simultaneously at many wavelengths.Explicit calibration methods such as ordinary least squares (OLS) and classical least squares (CLS) are often used when all of the constituent spectra can be individually measured or pre-calibrated using pure constituents.In contrast, implicit multivariate calibration techniques such as principal component regression (PCR) and partial least squares (PLS) allow the construction of a numerical model entirely based on spectroscopic and concentration data. 2 They are powerful because virtually no prior information is needed.Explicit or implicit, the goal of multivariate calibration is to obtain a vector consisting of regression coe±cients, b, such that an analyte's concentration, c, can be accurately predicted by taking the scalar product of b with an experimental spectrum, s: where lowercase boldface type denotes a column vector, uppercase boldface type a matrix; and the superscript T denotes transpose.The regression vector, b, has the same length as the measured spectra, and thus can provide spectroscopic information regarding the calibration.b is unique in an ideal noisefree linear system without constituent correlations, and the goal of both implicit and explicit schemes is to ¯nd an accurate approximation to b for the experimental system of interest.
Recently, we proposed a novel hybrid technique, called constrained regularization (CR), where implicit multivariate calibration is viewed as an inverse problem. 3 It is hybrid owning to salient features of both explicit and implicit techniques.CR requires only the spectroscopic and concentration data, yet allows the incorporation of the target analyte spectrum as prior information.Starting with the inverse mixture model as the forward problem, we de¯ne the inverse problem with solution b.Instabilities associated with the inversion process are removed by regularization, and prior information is included by means of a spectroscopic constraint.
5][6][7] Raman spectroscopy o®ers several advantages such as no need for reagents and separation, nondestructiveness to the specimen, high sensitivity with plasmonic enhancement, capability of qualitative and quantitative measurements, providing molecular structure information with high spatial resolution, etc. [8][9][10][11][12] We demonstrated that CR provides improved performance using numerical models as well as experimental Raman spectra acquired from human subjects.We show that with CR both the standard error obtained using leave-one-out cross validation (SEV) and the standard error of prediction (SEP) improve compared to results obtained using PLS.To the best of our knowledge, this is the ¯rst time that CR has been applied to a Raman skin model and in vivo data for noninvasive glucose sensing.

Theory
The forward problem for our calibration model is de¯ned by the linear inverse mixture model for a single analyte: where S is a set of calibration spectra, with each spectrum occupying a column of S, associated with several known concentrations of the analyte of interest that are expressed as a column vector, c, the jth element of which corresponds to the jth column of S. The goal of the calibration procedure is to use the set of data [S, c] to obtain an accurate b by inverting Eq. ( 2).The resulting b can then be used in Eq. ( 1) to predict the analyte concentration, c, of an independent sample by measuring its spectrum, s.The \accuracy" of b is usually judged by its ability to correctly predict concentrations prospectively via Eq.( 1).CR enabled a convenient way to accomplish inversion with the °exibility of incorporation of a spectral constraint, b 0 , as described as the minimization of a quadratic cost function, È 13 : with jjajj the Euclidean norm (i.e., magnitude) of a, and a spectral constraint that introduces prior information about b.The ¯rst term of È is the model approximation error, and the second term the norm of the di®erence between the solution and the constraint, which controls the smoothness of the solution and its deviation from the constraint.If b 0 is zero, the inversion is identical to standard Tikhonov regularization.As mentioned above, for ¼ 0 only least-norm (LM) solution is obtained.In the other limit, in which goes to in¯nity, the solution is simply b ¼ b 0 .In a previous volunteer study using PLS calibration, we retrospectively exploited the similarity between b and g, the spectrum of the analyte of interest, by calculating the correlation coe±cient of b with g. 14 The result indicated that signi¯cant correlation between b and g o®ers con¯dence that the calibration is based on the analyte of interest rather than spurious e®ects such as co-variations among constituents.Therefore, the spectrum of the analyte of interest is a good choice for the spectral constraint.

Experimental Methods
We have employed multivariate calibration techniques such as PLS to noninvasive glucose sensing using Raman spectroscopy. 2,7,146][17][18][19][20][21] Glucose is a convenient analyte to study because its concentration can be conveniently altered and monitored in human volunteers.As shown below, numerical simulations and in vivo data collected noninvasively from the forearm of volunteers demonstrate that CR, when used with appropriate constraints, provides a calibration with signi¯cantly improved prediction accuracy compared to PLS.
Numerical Study I.The numerical data were chosen to closely simulate the human data.The predominant Raman spectral features sampled from the forearm are indicative of skin [Fig.1(a)].To simulate the forearm spectrum, we employed a model composed of nine representative constituents of the skin-blood-tissue matrix and the spectrum of glucose dissolved in water (Fig. 2).The choice of these constituents was based on the known composition of skin, and the relative amplitudes were chosen to approximate those of the skin-blood-tissue matrix.Glucose was included at physiological concentrations of 70 mg/dL to 210 mg/dL.The resulting simulations exhibited Raman signal ratios of glucose to the total matrix varying from 0.2% to 0.6%, which is the typical range measured in skin. 22,23Such relative amounts of glucose were con¯rmed by studies in our laboratory employing minced samples of porcine skin, a good spectral model of human skin, with elevated levels of glucose.In simulating sample-tosample variations, we varied all of the background constituent concentrations in a random fashion (SD 5% of the design spectral weights in Table 1), ensuring that there is no signi¯cant correlation between pairwise model constituents (r 2 $ 0:02).An appropriate amount of Gaussian random noise (SD 130 counts), estimated from the volunteer data, was added to each noiseless sample.We de¯ne the signal as the norm of the spectrum of interest.Total Raman signal-to-noise-ratio (SNR) (12,000) and glucose Raman SNR (24-72) can then be calculated by dividing the norm of the total Raman signal (1:5 Â 10 6 counts) and the glucose signal (3120-9360 counts) by the noise magnitude (130 counts), respectively.Finally, to simulate reference concentration measurement error, Gaussian random error (SD 6 mg/dL) was added to the glucose concentrations, as well.The  model constituents and parameters are summarized in Table 1.Since these parameters are similar to their experimentally observed counterparts, we expect the numerical data to closely simulate the in vivo Raman spectra.
A typical Raman spectrum of human skin is shown in Fig. 1(a).The broad slowly varying background was ¯t to a ¯fth-order polynomial and subtracted o® from the spectrum.A total of 25 simulated Raman spectra are shown in Fig. 1(b), with glucose and other constituents varied within the above design constraints.As can be seen, they approximate the observed features of skin Raman spectra very well.Figure 1(c) shows the di®erence between the ¯rst two simulated spectra, magni¯ed 10Â, and Fig. 1(d) shows the model glucose spectrum at 90 mg/dL, magni¯ed 100Â.The Raman signature of glucose is not apparent in either the sample spectra nor their di®erence, thus necessitating the use of multivariate calibration techniques.
Numerical Study II.In Numerical Study I, all model constituents were varied in an approximately random fashion.To study the e®ectiveness of CR with spurious correlations present, a second numerical study was performed, with signi¯cant constituent co-variations present in the calibration sample set: Strong correlation (r 2 $ 0:72) between hemoglobin and glucose concentrations, and exponential decays in the total Raman signal level from the ¯rst sample to the last.(The volunteer spectra manifested behavior of this type. 14) All other model parameters were identical to those of the ¯rst numerical study.
Data Analysis.In both numerical studies, PLS, OLS and CR calibration methods were applied to these simulated spectra using code written in MATLAB (The MathWorks, Natick, MA).In each case, a set of 25 calibration sample spectra were generated and a speci¯c method was applied, with leave-one-out cross validation, to calculate the SEV and b vector.This b vector was then used to predict the concentrations of a set of 25 independent ¯eld samples, generated using the same model parameters, and the SEP calculated.SEV/P were calculated using the formula: with a j and c j the calculated and reference concentrations, respectively, and n the number of samples.This procedure was repeated 50 times and the results were averaged to obtain the estimates of mean error and variance.The selection of the regularization parameter was implemented as an automated search within the cross-validation framework, and the exact value has been selected when the SEV for a particular dataset is minimized.Therefore, it varied across di®erent datasets.
Volunteer Study.We applied PLS and CR calibrations to in vivo transcutaneous Raman spectra collected from 17 nondiabetic human volunteers. 14The full details of the protocol are given in Ref. 14. Brie°y, volunteers fasted for 12 h prior to the experiment.Shortly after the start of the experiment, the subjects drank a highly concentrated glucose solution (Sun-dex, 75 g glucose).Raman spectra were collected from a $ 1 mm 2 region of the forearm of a subject, excited by 830 nm light.The sampling volume was $ 1 mm 3 .Each spectrum was the average of 90 consecutive two-second acquisitions.Spectra were acquired every 5 min over a period of 2 to 3 h.During this period, the blood glucose concentration typically doubled and then returned to its initial level.Reference capillary blood samples were collected every 10 min and analyzed by means of a Hemocue glucose analyzer a precision of SD 6 mg/dL.Spline interpolation was used to provide reference values at the 5-min intervals.Raw spectra were smoothed and then the wavelength was calibrated using an indene reference.

Experimental Results
Numerical Study I. Table 2 summarizes the results from the application of each calibration method to the numerical data, using the results of PLS with 9 factors as a baseline.We observed the following: First, without regularization, the SEV is reduced only slightly, with solution led to simply the LN solution.The SEV is reduced only slightly when using CR with collagen I, the largest constituent of the simulation model, as the constraint, indicating that collagen I is a poor constraint.Finally, using CR with either the glucose spectrum or b OLS as the constraint, there is signi¯cant improvement in accuracy.Results for OLS are further discussed in Sec. 5.
Numerical Study II.Table 3 summarizes the results for the second numerical study in which the simulated data exhibit spurious correlations among model constituents.CR with appropriate constraints works much better than PLS or CR with poor constraints.This demonstrates that proper prior spectral information indeed improves the calibration results even when confounding e®ects such as spurious correlations exist.Since the b vector of OLS is not derived from the calibration sample set, spurious correlations have no e®ect on it.Thus, the OLS results are identical to those of Table 2.
Volunteer Study.Results for the application of PLS and CR without constraint or with glucose as the constraint to 17 volunteers of the transcutaneous study are summarized in Table 4 which shows the SEV's of data from two representative volunteers and the average over all 17 volunteers.However, due to nonlinear e®ects such as di®erences in tissue composition, auto°uorescence background, sampling volume, and physiological conditions, di®erent regularization parameters were selected to minimize the SEV for each individual.This also re°ects that better results can be obtained if the prior information is weighed more in some volunteers, while less in others.A potential future direction is to assess whether CR can better handle individuals with different conditions such as pre-diabetic or with skin issues.Averaging reduces variations among di®erent individuals due to physiological and skin diversities, unmodeled in the numerical studies.CR with glucose as the constraint shows the most signi¯cant error reduction compared to the result of PLS. Figure 3 shows the glucose spectrum, the PLS b vector, and the CR b vector with the glucose spectrum as the constraint.Both PLS and CR b vectors have been smoothed using a Savitzky-Golay 13-point ¯lter to remove high-frequency °uctuations that are likely caused by noise.

Discussion
PLS and PCR are valuable tools in multivariate calibration.For spectral analysis of linear chemical systems, they allow °exible experimental design and are powerful in extracting correlations between spectra and analyte concentrations.To utilize prior information, hybrid methods such as CR allows the °exiblity of incorporating spectroscopic information about the analyte of interest.In the broader context, regularization methods tend to perform slightly better than either PLS or PCR. 24A heuristic explanation is that regularization provides a continuous \knob", and therefore can be used to ¯nd a better balance between model complexity and noise rejection.Our results show that in addition to this intrinsic advantage, an even more signi¯cant improvement can be obtained by incorporating a meaningful solution constraint.
In Numerical Study I, we showed that CR with appropriate constraints outperforms PLS under \normal" conditions, i.e., no spurious correlations among constituents.In this case, spectroscopic noise and concentration error are the major factors determining prediction accuracy.In Numerical Study II, CR in the presence of spurious correlations was investigated.The results showed that CR with appropriate constraints gives rise to b vectors that are signi¯cantly less a®ected by such spurious correlations, and thus provide much better predictions than PLS.Further, compared to the scenario simulated in Numerical Study I in which no signi¯cant constituent co-variation was present, these results suggest that CR e®ectively reduces the detrimental in°uence due to spurious correlations.
Table 2 also lists the OLS results for constituents of di®erent spectral quality.For OLS (noiseless), the b vector was calculated from noiseless constituent spectra, whereas for OLS (noise-added), b was derived from constituent spectra with SNR similar to that of the calibration sample set used by the implicit methods.Therefore, the results of Table 2 show that the performance of CR with a well-chosen constraint approaches that of OLS (noise-added) with equivalent SNR.The fact that CR cannot predict as well as OLS (noiseless) is not surprising because, unlike OLS (noiseless), the CR b vectors are derived from noise-added data.The OLS (noiseless) result is considered the fundamental limit, given the SNR of the ¯eld spectra and the accuracy of reference concentration measurements.When the only error source is noise in the ¯eld spectra, the ultimate prediction error for OLS is inversely proportional to both the SNR and an overlap factor which quanti¯es spectral overlap among constituents.When other sources of error such as concentration error exist, the total prediction error can be estimated by the error terms in quadrature. 25,26To closely simulate the volunteer data in the numerical studies, 6 mg/dL was used as the reference concentration error, and apparently it dominates the OLS (noiseless) prediction error, 6.7 mg/dL.However, such a limit can be further improved with more accurately measured reference concentrations.
In the numerical studies, the size of the calibration sample set (25) was chosen to match the average number of samples obtained from each individual in the volunteer study.This makes the sample-to-PLSfactor ratio $ 3, smaller than the standard recommended value between 5 and 10. 27 To con¯rm that 9 is an appropriate number of factors for PLS, a numerical study with larger calibration sample size (50) was conducted.The SEP's for PLS, CR with glucose as the constraint, and CR with b OLS as the constraint, are 8.9, 7.9 and 7.8 mg/dL, respectively.Since more samples in the calibration set e®ectively increases the SNR for the calibration and reduces the degree of ill-posedness, the larger sample set gives lower prediction errors for all methods.However, the results suggest that CR with proper constraints can achieve the same quality calibration with fewer calibration samples as compared to PLS.
Finally, it appears that CR results in a much greater improvement over PLS when applied to the numerical studies than when applied to the volunteer study.One possible explanation is that the spectral constraints used for the volunteer study are obtained from clear samples, e.g., glucose in water or other individual pure constituents, whereas, transcutaneous Raman spectra are distorted by tissue turbidity.Therefore, the constraints employed in the volunteer study may not be the best choices.Incorporating tissue turbidity into constraints is currently under investigation in our laboratory.Another possible explanation is that the intense slowly varying background, not present in the spectral constraints, becomes a nonlinear confounding factor in the volunteer study.We are presently seeking a better way to model this background.

Conclusion
CR enjoys all advantages of implicit calibration while e®ectively avoids potential pitfalls such as spurious correlations among constituent concentrations.As has been demonstrated, there is °exibility in the choice of constraints, a convenient one being the spectrum of the analyte of interest itself.This °exibility is crucial because it is di±cult, if not impossible, to quantify high-¯delity spectra and the relative \importance" of individual constituents in living tissues.We have shown that CR signi¯cantly outperforms PLS when analyzing Raman spectra from human subjects, as well as simulated data with or without spurious correlations among constituents.In the volunteer study, CR reduced the SEV by $ 20% on average compared to PLS.In the numerical study, the error reduction is $ 73:5% without spurious correlations, and $ 61% in the presence of spurious correlations.
Further, there is no particular condition that must be satis¯ed in order for CR to function (apart from system linearity), which suggests that it can serve as a general method suitable for any scenario in which implicit calibration is performed.Its application to transcutaneous measurement of glucose concentration using near infrared absorption re°ectance/transmission spectroscopy is presently under investigation in our laboratory.

Fig. 3 .
Fig. 3. Comparison of the glucose spectrum (dashed line) with the smoothed PLS b vector (dotted line) and CR b vector (solid line) for data from one volunteer.

Table 2 .
Comparison of di®erent calibration methods with simulated data.Column 1 indicates the constraint used, columns 2 and 3 list the SEV and the SEP values obtained, with corresponding error reduction compared to PLS in parentheses, using the OLS (noise-added) and SEP (7.2) as the best achievable values.OLS (noiseless) indicates OLS in the absence of spectral noise; OLS (noise-added) denotes OLS with the same amount of spectral noise as for the PLS/CR data.Both OLS results include the same amount of concentration error as for PLS/CR data.

Table 3 .
Comparison of di®erent calibration methods using simulated data with substantial spurious correlations.See text for details.Column designations are the same as for Table2.