Monte Carlo lookup table-based inverse model for extracting optical properties from tissue-simulating phantoms using diffuse reflectance spectroscopy

Abstract. We present a Monte Carlo lookup table (MCLUT)-based inverse model for extracting optical properties from tissue-simulating phantoms. This model is valid for close source-detector separation and highly absorbing tissues. The MCLUT is based entirely on Monte Carlo simulation, which was implemented using a graphics processing unit. We used tissue-simulating phantoms to determine the accuracy of the MCLUT inverse model. Our results show strong agreement between extracted and expected optical properties, with errors rate of 1.74% for extracted reduced scattering values, 0.74% for extracted absorption values, and 2.42% for extracted hemoglobin concentration values.


Introduction
][3][4] Typically, DRS uses a fiber to deliver light to tissue.The delivered light is both scattered and absorbed by the tissue and is then recollected by another fiber a short distance from the source fiber.The collected light, or diffuse reflectance, contains quantitative information about tissue structure and composition.While the instrumentation for a DRS system is very simple, the accurate extraction of optical properties from the collected signal is a significant challenge.An accurate tissue model is needed to relate the collected signal to tissue optical properties.One method for analyzing diffuse reflectance spectra relies on the solution to the diffuse approximation of the radiative transport equation.However, the diffusion approximation is not valid for short source-detector separations or in highly absorbing tissues. 5ecause most cancers originate in the epithelial layer, it is necessary to use probes with a short source-detector distance in order to sample photons that travel through the epithelial layer. 6In addition, angiogenesis, an indicator of early cancer, can cause a significant increase in absorption due to blood.Unfortunately, the diffuse approximation is not valid in these regimes.Recently, reflectance lookup tables have been used to analyze diffuse reflectance spectra. 7,8These lookup tables are created in two different ways: experimental measurements of phantoms with known optical properties 7 or Monte Carlo simulations. 9Creating a lookup table (LUT) with experimental measurements has the advantage of incorporating unknown system responses into the LUT.However, the creation of an experimental LUT is time-consuming, and the accuracy is dependent on the skill and experience of the investigator.The Monte Carlo method is especially useful for creating an LUT, because it provides the ability to model complex probe geometries and tissue structures.However, intensive computation is required to achieve results with desirable variance, which can make it extremely time-consuming to populate an LUT containing thousands of values.
Much prior work has been undertaken to improve the speed and efficiency of the Monte Carlo method for modeling light transport in turbid media.These methods can be separated into three groups: 1. methods that use the information from a small set of Monte Carlo simulations and scale the results to calculate a wide range of optical properties 10 2. methods that use geometry splitting techniques to increase the fraction of useful photons 11 3. methods that parallelize the Monte Carlo simulations. 12e first set of methods has the advantage of not requiring a large number of simulations to create an LUT.However, errors arise because the photon trajectory information necessary to perform the scaling operation can be recorded only in several depth intervals with finite widths. 13Geometry splitting techniques decrease the number of photons needed for a simulation to converge, but their implementation is very difficult for complex probe geometries.The third set of methods uses classical Monte Carlo simulations and therefore does not make any sacrifices in accuracy, flexibility, or implementation difficulty.Speedup is achieved by simulating multiple photons simultaneously on different processors.Because each photon is independent of every other photon, this problem is considered "embarrassingly parallel."Alerstam et al. have shown that general-purpose graphic processing units (GPGPUs) can increase the speed of Monte Carlo simulations of photon transport by three orders of magnitude on a relatively inexpensive GPU when compared to the sequential implementation. 12We present, for the first time, a Monte Carlo LUT-based model where all values in the LUT were created by independent Monte Carlo simulations by using a parallel implementation on a GPGPU.

Creation of the Lookup Table
A two-dimensional Monte Carlo code written in ANSI C implemented on an NVIDIA GTX 560 Ti GPU with NVIDIA's Compute Unified Device Architecture (CUDA) was used to simulate photon reflectance in a single-layer tissue model on 386 parallel threads. 12The multiply with carry random number generator was used.The refractive index above the tissue was set to 1.452 to match the refractive index of the fiber, and the refractive index of the tissue was set to 1.33 to match the refractive index of the phantoms.To test the effect of errors in the refractive index, we created LUTs with different refractive indices and repeated the extraction of optical properties with the different LUTs.We found a 5% error rate in the refractive index corresponds to error increases of 1.3% for extracted μ 0 s values and 0.8% for extracted μ a values.To prevent photons from exiting the tissue volume, the radius and width of the tissue volume were set to 3 cm.A total of 1 × 10 6 photons were launched to obtain the impulse response.To ensure stochastic noise would be sufficiently low in LUT locations with high albedo, 100 separate MC simulations were performed with the optical properties in the LUT that would give the lowest value of reflectance.We found that using 1 × 10 6 photons reduced the standard deviation of the 100 different reflectance values to less than 0.5% of the mean.The diffuse reflectance for our specific probe geometry was then calculated by convolving the impulse response with the beam profile. 14Our probe was modeled using a Gaussian shaped beam profile of collimated light with a diameter of 200 μm, a detector diameter of 200 μm, and a sourcedetector separation of 250 μm.The diffuse reflectance values for all physiological realistic combinations of scattering and absorption were calculated using the GPGPU Monte Carlo implementation with the tissue and probe geometry described above.Twenty increments (N) were used for both scattering (0 to 50 cm −1 ) and absorption (1 to 50 cm −1 ), meaning a total of 400 separate Monte Carlo simulations were needed to create the LUT.The Henyey-Greenstein phase function was used for sampling scattering angles.The scattering anisotropy (g) was set to 0.85 for all simulations.For the range of g values present in human tissue (g > 0.8), it has been shown 10 that the diffuse reflectance will be the same for any values of μ s and g that generate the same μ 0 s .The resulting MCLUT is shown in Fig. 1.It took 2 min to run the 400 separate Monte Carlo simulations.

Forward and Inverse Models
For the forward model used to generate diffuse reflectance spectra, the reduced scattering coefficient was constrained to the form where λ 0 ¼ 630 nm.The absorption coefficient was calculated using where ε i ðλÞ is the wavelength-dependent extinction coefficient of a chromophore, C i is the concentration of that chromophore, and N is the number of chromophores.Depending on the type of tissue sampled and the wavelength range of interest, any number of chromophores can be used to calculate μ a ðλÞ.Once μ 0 s ðλÞ and μ a ðλÞ are calculated, the MCLUT can be used to generate a modeled reflectance spectrum.Cubic splines were used to interpolate between values in the LUT. Figure 2(a) shows the forward model of diffuse reflectance.
Fig. 1 The resulting lookup table [Rðμ 0 s ; μ a Þ] created using 400 separate Monte Carlo simulations.Each Monte Carlo simulation was used to calculate a reflectance value for a given scattering coefficient and absorption coefficient combination.shows the inverse model used to fit our diffuse reflectance spectra.First, we made an initial guess for the optical properties, and then the forward model was used to generate a spectrum.Next, the sum of squares error between the predicted reflectance and the measured reflectance was calculated using where δ is the sum of squares error, K is the number of wavelength points, R m is the measured spectrum, R e is the modeled spectrum, and λ i is the wavelength.The parameters are then iteratively updated until the sum of squares error is minimized.An interior-point nonlinear optimization routine provided in the MATLAB optimization toolbox (Mathworks, Nattick, Massachusetts) was used as the optimization algorithm.The average fit time was 1.3 s, with the number of iterations limited to 5000 and the termination tolerance on the error function set to 1 × 10 −3 .Because the modeled spectra are in absolute units (photons counted) and the measured spectra are in units relative to a baseline calibration measurement, it is necessary to perform a calibration so that the modeled and measured spectra can be compared.Additionally, the calibration corrects for wavelengthdependent responses in the experiment that are not accounted for in our forward model.We performed the calibration by taking the ratio of a modeled spectrum and a measured spectrum with the same optical properties.Then, to make the measured spectra equivalent to the modeled spectra, all measured spectra were multiplied by this ratio.To ensure that the choice of optical properties used in the calibration step did not bias the results, the mean of the ratios for all spectra used in the validation set was used for calibration.The mean calibration ratio is shown in Fig. 3(a).

Validation and Results
To test the performance of our MCLUT-based inverse model, we created 21 tissue phantoms with hemoglobin (Hb) (Sigma-Aldrich) as the absorber and polystyrene beads (diameter ¼ 1 μm) as the scatterer.Hb concentration [(Hb)] ranged from 0 to 3 mg∕ml, and the reduced scattering coefficient [μ 0 s ðλ 0 ¼ 630 nmÞ] ranged from 6.4 to 27.5 cm −1 .We used Mie theory to calculate the μ 0 s of the tissue phantoms.For our inverse model, we assumed the absorption in the visible range was due to oxy-hemoglobin.We measured the optical density for the HbO 2 solution using a spectrophotometer and calculated the absorption spectrum using Beer's law.Because the addition of HbO 2 dilutes the solution, a small change in μ 0 s was accounted for when calculating the known values for μ 0 s .The DRS system consisted of a xenon flash lamp (Model: E6611, Hamamatsu) as the light source, a spectrograph (Model SP2150i, Princeton Intruments) and camera (Cool-SNAP, Photometrics) as the spectrometer; and a fiber optic probe with the geometry described above (FiberTechOptica, Ontario, Canada).The diffuse reflectance spectrum and its associated fit can be seen in Fig. 3(b) and shows that the inverse model can accurately fit the experimental data.Figure 4(a) and 4(b) shows scatter plots of the extracted versus expected for μ 0 s ðλ 0 Þ and [Hb], respectively.The solid line in each graph is the line of perfect agreement.The results indicate there is excellent agreement between the extracted and expected optical properties.The MCLUT inverse model estimated the optical properties over a wide range with average root-meansquare percent errors of 1.74% for μ 0 s , 0.74% for μ a , and 2.42% for [Hb].We compared the performance of our MCLUT-based model to an experimental LUT-based model.The MCLUT model was able to estimate μ 0 s and μ a with decreases in percent error of 3.16% and 10.86%, respectively, when compared to the experimental LUT model. 7ur inverse model was then tested on previously collected data from a clinical feasibility study 15 to illustrate the application of the model for noninvasive detection of skin cancer.Figure 4(c) shows representative spectra from two groups: clinically normal and basal cell carcinoma (BCC).The plot shows good agreement between the MCLUT fit and the measured in vivo spectra.For this analysis, the absorption coefficients calculated in the forward model were determined by melanin, deoxyhemoglobin, and oxy-hemoglobin.

Conclusions
Although other inverse models have recently been developed for extracting optical properties from diffuse reflectance spectra, this work represents the first Monte Carlo lookup table-based inverse model where the LUT was generated entirely by Monte Carlo simulation.Previously, the amount of time required to generate a LUT entirely by Monte Carlo simulation made this technique infeasible.However, recent advances in GPGPU computing have allowed parallel Monte Carlo implementations capable of running three orders of magnitude faster than traditional, serial implementations of Monte Carlo simulations.By creating an LUT entirely by Monte Carlo simulation, our method is not subject to the errors that arise from using either the diffusion approximation or the Monte Carlo scaling method.When compared to an experimental LUT, our method was more accurate, but, more importantly, it has the advantages of being easier to implement and repeatable.This model can also be adapted to more complex probe and tissue geometries.

Fig. 2
Fig.2Flowcharts for the forward (a) and inverse (b) models of diffuse reflectance used to create the modeled spectrum and to fit the MCLUT model to the reflectance data.

Figure 2 (
Figure 2(b) shows the inverse model used to fit our diffuse reflectance spectra.First, we made an initial guess for the optical properties, and then the forward model was used to generate a spectrum.Next, the sum of squares error between the predicted reflectance and the measured reflectance was calculated using

Fig. 3 (
Fig.3(a) A diffuse reflectance spectrum (dashed) [μ 0 s ðλ 0 Þ ¼ 25.4 cm −1 and ½Hb ¼ 1.5 mg∕ml] and the associated MCLUT-fit (solid).(b) Mean of all ratios of measured and modeled diffuse reflectance spectra with the same optical properties.This ratio was used for calibration.