Perturbation and differential Monte Carlo methods for measurement of optical properties in a layered epithelial tissue model

. The use of perturbation and differential Monte Carlo (cid:1) pMC/ dMC (cid:2) methods in conjunction with nonlinear optimization algorithms were proposed recently as a means to solve inverse photon migration problems in regionwise heterogeneous turbid media. We demonstrate the application of pMC/dMC methods for the recovery of optical properties in a two-layer extended epithelial tissue model from experimental measurements of spatially resolved diffuse reﬂectance. The results demonstrate that pMC/dMC methods provide a rapid and accurate approach to solve two-region inverse photon migration problems in the transport regime, that is, on spatial scales smaller than a transport mean free path and in media where optical scattering need not domi-nate absorption. The pMC/dMC approach is found to be effective over a broad range of absorption (cid:1) 50 to 400% (cid:2) and scattering (cid:1) 70 to 130% (cid:2) perturbations. The recovery of optical properties from spatially resolved diffuse reﬂectance measurements is examined for different sets of source-detector separation. These results provide some guidance for the design of compact ﬁber-based probes to determine and isolate optical properties from both epithelial and stromal layers of superﬁcial tissues. ©


Introduction
There has been significant interest in the noninvasive determination of the morphology and biochemical composition of biological tissues through the measurement of optical properties. Epithelial tissues are a major target because they represent the site from which more than 85% of all cancers originate. 1 Epithelial tissues line organ surfaces with a single or stratified layer of closely packed cells and are supported by underlying stromal tissue composed primarily of collagen and other extracellular matrix proteins. It is well known that the vast majority of epithelial cancers are preceded by distinct biochemical and morphological changes in both the epithelia and stroma. 2 The thickness of the cellular epithelial layer is typically Շ500 m and lies within the transport regime with respect to optical scattering. 3,4 On one hand, the thickness of the cellular epithelium is of the same order or larger than the single scattering mean free path l s ͑=1 / s ͒ of photons in the visible and near-IR, typically l s ϳ 5 to 100 m. Thus, when using light to interrogate both the epithelium and the underlying stroma, one will invariably detect photons that have multiply scattered within the tissue; often interacting with both these layers. On the other hand, the epithelial layer thickness is typically smaller than the transport mean free path l * ͓=1 / ͑ a + s Ј͔͒ and precludes the use of models based on the diffusion approximation to radiative transport to analyze the optical signals and extract optical properties. As a result, the use of optical approaches to detect and separate the optical properties of these two tissue layers continues to present a major challenge in biomedical optics.
The accurate characterization of both epithelial and stromal tissues is critically important because as epithelial cancer progresses, the optical properties of the epithelial and stromal layers can be altered in markedly different ways. For example, the development of dysplasia and epithelial neoplasm often results in an increase in epithelial scattering due to the enlargement of nuclear size, increased nuclear/cytoplasmic ratio and hyperchromasia. Simultaneously, stromal scattering can decrease due to degradation of collagen fibers resulting from the release of matrix metalloproteases ͑MMPs͒ and stromal absorption can increase due to the release of angiogenic factors such as vascular endothelial growth factor 1,2 ͑VEGF͒.
To address these issues, many investigators have examined the sensitivity of reflectance measurements made at specific source-detector separations to optimize their sensitivity to changes in optical properties in homogeneous media. [5][6][7] This knowledge has spurred the design of optical probes to isolate or differentially sample the superficial epithelial tissue versus the supporting stromal tissue layer or vice versa. Several studies have examined the impact of many fiber probe design characteristics including fiber diameter, numerical aperture, and source-detector ͑s-d͒ separation on the sensitivity of fluorescence and reflectance measurements in these two layers. [8][9][10][11] Recently, many have proposed novel fiber probe configurations that potentially isolate contributions of the epithelium to the detected signal 12,13 or provide depth-selective reflectance spectra to separate epithelial scattering from stromal scattering. 14,15 Others, including the groups of Backman and Sokolov, have used polarization gating techniques to detect photons that travel "short path lengths" within the tissue from those that travel "long path lengths" to provide differential sensitivity to the superficial epithelial layer versus the deeper stromal layer. [15][16][17][18] While such approaches no doubt enhance the contribution of either the epithelial or stromal layer to the detected optical signal, the fact that the epithelial tissue thickness is typically of the same order or larger than the single scattering mean free path l s prevents these probe designs from fully separating and isolating the contributions of both the epithelial and stromal tissue layers to the detected signal.
These innovative developments in probe design to better interrogate epithelial and stromal tissues would be complemented greatly by accurate modeling of radiative transport in heterogeneous tissues on submillimeter length scales. If available, such models could assess the epithelial and stromal contributions to spatially resolved reflectance measurements and perhaps recover the distinct optical properties of the two layers. Such a capability would improve greatly the potency of diffuse reflectance measurements to detect and characterize the physiological changes occurring in both epithelial and stromal layers during dysplastic and neoplastic transformation. As noted earlier, the epithelial thickness is typically intermediate between l s and l * . Moreover, the epithelium and sometimes the underlying stroma may not be strongly scattering. Both these features demand the use of models that can provide both forward and inverse predictions of radiative transfer in the transport regime, i.e., outside the range of the diffusion approximation. 4 The development of such radiative transport models is difficult. Optical diagnostic modalities based on multiple light scattering to recover optical properties of spatially heterogeneous tissues have relied on methods that either invoke the diffusion approximation to radiative transport or utilize hybrid diffusion/Monte Carlo approaches. Unfortunately, the validity of these approaches are limited to highly scattering tissues ͓i.e., ͑ s Ј/ a ͒ տ 20͔ and/or spatial scales significantly larger than a transport mean free path l * ͓ϵ1/͑ a + s Ј͔͒. Specifically, the numerous methods that have been developed to determine optical properties of layered media from spatially resolved, time domain and/or frequency-domain diffuse reflectance measurements cannot be applied to epithelial tissues due to the submillimeter thicknesses of the tissue layers involved. [19][20][21][22][23][24][25][26][27][28][29][30] More recently, several investigators have used Monte Carlo methods to examine the impact of epithelial and/or stromal optical property changes as well as changes in epithelial thickness on the resulting photon migration signal. [31][32][33] However, these studies have not addressed the inverse problem of optical property determination from the photon migration measurements. Thus, the availability of methods to enable the recovery of optical properties in low/ moderately scattering tissues on spatial scales comparable to l * from the signals measured by the new generation optical probes aimed to examine epithelial tissues would be quite valuable for the accurate determination of layered optical properties.
These facts motivated us to develop a computational platform to address both forward and inverse problems in a framework faithful to the radiative transport equation. This platform combines perturbation and differential Monte Carlo ͑pMC/dMC͒ methods to provide transport rigorous estimates of optical signals as well as extraction of optical properties. These features make it applicable to systems with arbitrary optical properties possessing heterogeneities of characteristic dimension approaching l s . Our group proposed and computationally validated this approach in an earlier study. 3 This study demonstrated the use of pMC/dMC methods in conjunction with a gradient-based optimization algorithm to process MCsimulated spatially resolved frequency-domain reflectance measurements to provide for the efficient and rapid recovery of optical properties in a layered epithelial tissue model. In this technique, the optimization determined the optical properties of the heterogenous tissue regions that result in a radiative transport prediction that best matched the simulated measurements. The radiative transport predictions required for the optimization procedure are provided using a pMC algorithm. The pMC algorithm provides these estimates by rescaling the weights of photon biographies stored from the MC simulation of a related "background" medium in a way consistent with the assumed optical properties of the heterogeneous tissue. In addition a dMC algorithm provides a sensitivity analysis of the measured optical signals relative to the optical properties of each heterogeneous region to determine how the optical properties in these regions should be changed for each iteration of the gradient-based optimization. The advantage of using this combined pMC/dMC approach was demonstrated by rapidly determining optical properties in a two-layer turbid medium under conditions in which diffusion-approximationbased approaches are not tenable.
Since our initial introduction of the pMC/dMC technique, other groups have adopted its use for the examination of other problems in biomedical optics. For example, the studies of Kumar and Vasu 34 and Yalavarthy et al. 35 applied this technique to both simulated and measured data sets of timeresolved diffuse reflectance/transmittance for optical property image reconstruction of low-scattering turbid media. These studies demonstrated the superior performance of the pMC/ dMC technique relative to diffusion-approximation-based reconstruction methods.
Given these recent advances in the application of the pMC/ dMC technique to optical diagnostics, the objectives of this paper are twofold. First, we wish to demonstrate the ability of the pMC/dMC approach to recover optical properties from experimental measurements of steady-state ͑CW͒ spatially resolved diffuse reflectance ͑SRDR͒ of two-layer turbid media. The successful utilization of pMC/dMC techniques to recover optical properties from CW reflectance data is not evident because this data does not contain information on the dynam-ics of optical transport in the medium. The ability to recover optical properties using CW data is of practical significance due to the reduced expense and ease of use of CW measurement systems relative to time-domain or frequency-domain systems. Second, we wish to explore the impact of different sets of s-d separation on the ability to recover optical properties of the superficial epithelial layer versus the deeper stromal layer. The experimental system utilizes tissue phantoms whose optical absorption and scattering properties are similar to those of epithelial tissues with layer thicknesses smaller than a transport mean free path. In what follows, we provide details of the theory that forms the basis of pMc/dMc methods, describe the methodology for determination of optical properties and the experimental system in which the measurements were performed. We then provide the results of our optical property recovery and discuss the general applicability of this approach.

pMC/dMC Methods
The pMC algorithm implemented in this paper employs a conventional MC simulation for a fixed s-d geometry with homogeneous background optical absorption and scattering coefficients a and s . We launch 5 ϫ 10 6 photons at the origin of semi-infinite media and save the biographies ͑path length and scattering locations͒ of the detected photons. For a narrow beam normally incident onto a tissue, intercollision distances s are chosen from the probability density function t exp ͑− t S͒ as the photons propagate within the tissue, where t = a + s . At each collision point, a new scattering angle is sampled from the Henyey-Greenstein phase function.
The weight of each photon is reduced by a factor of ͑ s / t ͒ at each collision. Thus, a photon that has undergone k collisions within a homogeneous tissue prior to detection carries a weight W = A͑ s / s ͒ k , where A is a factor that accommodates the numerical aperture of the detector and the refractive index mismatch between the detector and the tissue.
A valuable feature of pMC is that it provides for the rapid determination of detected photon weights for a system that can be considered perturbed relative to a background system. Specifically, let us consider a heterogeneous tissue consisting of N perturbed regions each with optical properties a,i and s,i . If the optical properties of each of these regions can be expressed in terms of a perturbation to the background optical properties, i.e., a,i = a + ⑀ a,i and s,i = s + ⑀ s,i , ͑see Fig. 1͒ the detected photon weight in the second system can be estimated by where j i is the number of collisions the photon experiences in the ith perturbed region, and S i is the path length of photon travel within the ith perturbed region. 3 The pMC method enables the simulation of several different, but related transport problems at a fraction of the computational time and cost required for separate MC simulations. 36,37 In our application, Eq. ͑1͒ is used to obtain estimates for the change in the detected photon signal in a perturbed system using the saved photon biographies from a background simulation. Note that the background system is not required to be spatially homogeneous but can instead be composed of several structures containing different optical properties. The expected value of Eq. ͑1͒ summed over all photons is the signal produced by the perturbed system. Thus, the implementation of pMC requires only a single MC simulation for the background system that can serve as a basis to determine the optical signal of any number of related systems so long as these systems can be related in terms of regions whose optical properties that can be expressed in terms of the background system.
To find the optical properties of the perturbed regions of interest, derivative information must be supplied to the gradient-based optimization. This information is obtained by taking the derivative of the Eq. ͑1͒ with respect to either the absorption perturbation ⑀ a,i or scattering perturbation ⑀ s,i of the ith heterogeneity. These sensitivities can be derived for all N heterogeneous regions in the perturbed geometry. The closed form expressions for these derivatives follow and constitute dMC.

Epithelial Tissue Model
The SRDR measurements are performed using turbid liquid phantoms with 10% Intralipid ͑B. Brown, Irvine, California͒ solution providing the scattering and nigrosin ͑Sigma-Aldrich͒ as the optical absorber. The optical phantoms were designed to reproduce optical properties of normal cervical tissue measured by Hornung et al. 38 We consider the epithelial tissue to be composed of a thin ͑500 m͒ cellular layer supported by a thicker stromal layer. The mean optical properties found by Hornung et al. at = 849 nm: a = 0.034 mm −1 and s Ј=0.611 mm −1 were used. This set of parameters corresponds to l * = 1.55 mm, ͑ s Ј/ a ͒ =18, and an epithelial layer The consistent preparation of a two-layer liquid phantom with a 500-m top-layer thickness proved to be difficult. As a result, we constructed an epithelial tissue phantom with spatial dimensions that are roughly double in size. To preserve optical similarity in the phantom system, the value of ͑ s Ј/ a ͒ was kept constant in this "extended" tissue model. Thus, in the extended epithelial tissue model, the thickness of the epithelial layer is 1 mm and the optical properties of the background were chosen to be a Ј=0.016 mm −1 and s Ј = 0.284 mm −1 . These optical properties ensure that the extended tissue model has the same ratio of reduced scattering to absorption ͓͑ s Ј/ a ͒ =18͔, while doubling the transport mean free path to l * = 3.33 mm.
Based on these background optical properties, we prepared phantoms with perturbed optical absorption or scattering properties in either the top or bottom layer. Note that although perturbations in both a and s can be examined simultaneously using the proposed pMC/dMC method, we chose to examine perturbations in only one property in either of the layers to test the effectiveness of the two parameter inversion algorithm and its susceptibility to "crosstalk" between absorption and scattering properties. This led to the examination of four cases: perturbation in the ͑1͒ top layer absorption, ͑2͒ bottom layer absorption, ͑3͒ top layer scattering, and ͑4͒ bottom layer scattering. We examined perturbations to the background optical properties in the range of 50 to 400% for a and 50 to 150% for s . Depiction of these four cases is presented in Fig. 2 along with the range of optical properties examined. The detailed optical properties for each absorption and scattering perturbation case considered are listed in Tables 1 and 2. Note that collectively the range of optical properties tested span a broad range of ͑ s Ј/ a ͒Ϸ4 to 40. Figure 3 provides a schematic of the experimental setup for the SRDR measurements. A cylindrical container of 80-mm diameter and 100-mm height was filled with an optically turbid medium composed of 10% Intralipid, nigrosin, and deionized water. The container was composed of two compartments, with top and bottom layers separated by a thin ͑Ͻ100 m͒ sheet of plastic film. The bottom compartment was filled with the phantom solution with optical properties corresponding to the stromal layer while a premeasured volume of the phantom solution for the epithelial layer was poured into the top compartment to provide a thickness of 1 mm. A He-Ne laser ͑JDS Uniphase, Mountain View, Cali-fornia͒ emitting at = 632.8 nm was coupled to a 200m-diam multimode optical fiber using a collimating lens. A neutral density filter was placed between the laser and fiber to attenuate the intensity of the light source delivered to the sample. To ensure that specular reflection was not captured, a fiber optic illuminator ͑Fiberguide, Stirling, New Jersey͒ was designed and fabricated to deliver light to the medium and block specular reflection from reaching the camera. The illuminator was composed of a right-angle microprism 0.5 mm in size coupled to a fiber with numerical aperture of 0.12. The prism was used for deflecting light from the fiber for perpendicular illumination of the medium. The fiber and prism were placed in the center of a black, rectangularly shaped housing ͑30 mm L ϫ 4 mm Wϫ 5 mm H͒. In making the measurement, the microprism was placed flush with the tissue phantom surfaces, as shown in Fig. 3. The use of this illuminator enabled the acquisition of quality images at locations proxi- mal to laser source without saturation and enabled the utilization of the full dynamic range of the camera.

Measurement of SRDR
The SRDR was measured by placing the sample surface in the focal plane of a 35-mm photographic lens ͑Nikon, Japan͒ coupled to a thermoelectrically cooled 16-bit CCD camera ͑Photometrics, Tucson, Arizona͒. A 3.7ϫ 3.7-cm region on the sample surface was captured by the camera and the resulting image was binned in groups of 2 ϫ 2 pixels to improve the SNR. Control of the CCD readout, image storage, and data processing were provided by a personal computer ͑Xi computer, San Clemente, California͒. The spatial resolution was determined to be Ϸ70 m using a U.S. Air Force resolution standard ͑Edmund Scientific͒. Pixel intensity as a function of radial position away from the source was acquired by the camera 15 times to provide an average and standard deviation for the SRDR measurement. Figure 4 provides a contour plot from a captured image of a liquid optical phantom and also provides the raw reflectance data obtained from a line scan of the image.
To convert the measured intensity to an absolute SRDR measurement, we must account for factors such as the optical point spread function of the system, sensitivity of the CCD, spatial variations in the CCD response, and the instrument response of the system. These factors were accounted for using a calibration procedure described previously by Pham et al. 39 This procedure consists of two steps. First, images of the diffuse reflectance were acquired from a reference phantom with known optical properties. These calibration images were taken using the same settings as those used for measurements of the layered tissue phantoms. The instrument response of the system was then calculated using the data obtained from the reference phantom and a prediction of the diffuse reflectance generated by a MC simulation using the same optical properties of the reference phantom. The raw image measured by the CCD ͑I R ͒ is a convolution of the true image ͑I T ͒ and the instrument response ͑I IR ͒, which can be described as When Eq. ͑4͒ is Fourier transformed into the spatial frequency domain, the convolution relation is converted into a simple multiplication, so that the instrument response can be determined by the division of the Fourier transform of the raw image by the Fourier transform of the "true" image, as provided by the MC simulation. Using this instrument response, any measured intensity data can be converted into absolute reflectance ͑R abs ͒ using the following relationship: where F and F −1 represent the Fourier and inverse Fourier transforms, respectively. This reflectance data is then fed into

Source-Detector Configurations
To study the quality of the optical property recovery with the range and spacing of detector locations used in the SRDR measurements, we analyzed the measured data using specific sets of s-d separations. The four configurations chosen are illustrated in Fig. 5 and listed in Table 3. All sets consisted of 10 detectors each, DS1, DS2, and DS3 have detectors spaced at uniform intervals with the minimum s-d separation at r min = 0.5 mm and the maximum s-d separation set at r max = 2.75, 5.00, and 9.50 mm, respectively, This range in the maximum s-d separation corresponds to a range of approximately l to 3l * of the background medium. The fourth configuration DS4 spans the same spatial range as the DS3, except that the spacing between the detectors is not constant but scaled exponentially. We examined this case to provide greater sampling in the region closest to the source where the intensity of the diffuse reflectance falls sharply and sparser sampling at more distal locations where the intensity is falling less rapidly.

Determination of Optical Properties
Our approach to the determination of optical properties is depicted schematically in Fig. 6. The capital letters shown in the square brackets below refer to the individual steps shown within this schematic. We use a conventional MC simulation to generate photon biographies using the s-d geometry used in the measurement. In this study, the background optical properties used were those found by Hornung et al., as mentioned in Sec. 3.1. However, in general, the background optical properties could either be properties of the superficial epithelial layer determined through a technique such as differentialpath-length spectroscopy 12,13 or that of the supporting stromal tissue determined using SRDR measurements at s-d separations larger than a transport mean free path. Moreover, in this paper we assume that the thickness of the epithelial layer is known either from anatomical information, another imaging modality such as high-resolution ultrasound, or using optical measurement using the method proposed recently by Gandjbakhche's group. 31 For post-processing of the photon biographies, we save only information needed to compute the perturbation in the detected photon weight and the sensitivities as specified by Eqs. ͑2͒ and ͑3͒. This requires only the storage of both the number of collisions and total photon path length within each tissue region as well as the angle of photon exit from the tissue upon detection. With this photon biography information generated and stored ͓A͔, pMC is used to calculate a predicted diffuse reflectance signal I p at each detector location r i using an initial set of optical properties in the perturbed region that are identical to the background ͑i.e., ⑀ a = ⑀ s =0͒ ͓B͔. Note that in this case, N =1 as we are perturbing only the properties of the stromal layer or the epithelial layer. The algorithm then compares the actual experimental measurement I m ͑r i ͒ to the radiative transport prediction provided by pMC using a least-squares metric 2 ͓C͔; where M is the number of measurement locations, and i is the standard deviation of the diffuse reflectance measurement taken at location r = r i . The 2 metric is tested for its convergence to the minimum within a Levenberg-Marquardt algorithm as described by Press et al. 40 ͓D͔. The convergence criteria of the two-parameter fit requires that the vector distance between two successive iterates is less than a predefined value= 0.001. If the convergence criteria is not met, a new set optical properties for the perturbed region ͑ a , s ͒ is determined using the derivative information provided by dMC using Eqs. ͑2͒ and ͑3͒ ͓E͔. We then used pMC again to generate a new estimate for the predicted reflectance signal I p ͓B͔. For recovery of the top layer optical properties, we supplied the measured SRDR data at the s-d separations specified by the detector sets described above ͑see Fig. 5 and Table 3͒ to the pMC/dMC inversion algorithm to determine the toplayer optical properties. The inversion results for each of the detector sets DS1 to DS4 are shown in Figs. 8͑a͒-8͑d͒ and provide the recovered top layer absorption ‫͒ؠ͑‬ and scattering ͑•͒ coefficients as a function of the absorption perturbation in the top layer. For each detector set, we provide the results of optical property recovery from SRDR measurements made in 10 different phantoms where s is held fixed and the top layer a is varied between 50 and 400% of the background. The diagonal line indicates the true values of a , and the horizontal line indicates the true value of s . All the detector sets successfully recover the true absorption and scattering properties in the perturbed top layer and decouple the constant Table 3 Listing of the absolute detector locations ͑r͒ and their relative location based on transport mean free path ͑l * ͒ of the four detector configurations used for the recovery of optical properties of two-layer extended epithelial tissue model.    value of s , Detector sets DS2, DS3, and DS4 all provide data to the pMC/dMC algorithm of sufficient quality to recover the true absorption and scattering properties of the bottom layer over the full range of absorption perturbation. Specifically, for DS2, DS3, and DS4 detector sets, the bottom layer a is predicted with a relative error of less than 20% throughout the full range of the 50 to 400% a perturbation. The unperturbed s is of the bottom layer is recovered within 15% relative error. The results obtained using data provided by the DS2 detector set are particularly notable in that the accuracy in optical property recovery is similar to the DS3 and DS4 detector sets even though the span of the DS2 detector set is Ͻ5 mm, about half the span used in DS3 and DS4. However, the reduced detector span of DS2 results in slightly higher uncertainties in the recovered a values. Figure 11 provides the SRDR measurements for a series of phantoms in which the optical scattering in the top layer is varied in the range 50 to 150% of the background value.

Top Layer Scattering Perturbation
These results show that increases in the top layer scattering result in increases in the SRDR, as indicated by the direction of the arrow in the figure. Moreover, increases in the SRDR are most noticeable at detector locations proximal to the source ͑r Ͻ 6 mm or Շ2l * ͒ and becomes insignificant at more distal locations. This occurs because detector locations proximal to the source receive photons that have interacted much more with the top layer and thus experience a greater probability of being backscattered. However, as the s-d separation increases, the photons that are detected have spent an increasing amount of time in the lower layer and thus the probability of being detected is not affected greatly by changes in top layer scattering. This lack of sensitivity of the SRDR at larger s-d separations to top layer scattering changes indicates that SRDR measurements at these locations possess little informa- Fig. 9 Experimentally measured SRDR data from extended epithelial tissue model with an absorption perturbation in the bottom layer. Data points shown correspond to detector set 4 ͑DS4͒. The relative error in the measured reflectance is smaller than the symbol size and less than 5% in all cases. The direction of the arrow indicates an increase of the absorption perturbation. These effects are clearly evident in the results for optical property recovery by the pMC/dMC algorithm, as shown in Figs. 12͑a͒-12͑d͒, which provide the recovered top layer absorption ‫͒ؠ͑‬ and scattering ͑•͒ as a function of the scattering perturbation in the top layer. For each detector set, we provide the results of optical property recovery from SRDR measurements made in 10 different phantoms where only the top layer s is varied between 50 and 150% of the background. The diagonal line indicates the true values of s , and the horizontal line indicates the true value of a . While the top layer s ͑•͒ is recovered accurately by all detector sets for scattering perturbations in the range 70 to 130%, the recovery of a ‫͒ؠ͑‬ is not decoupled from the changes in s . For scattering perturbations of 50 to 70%, a is successfully decoupled from s and predicted with a relative error of Ͻ20%, but the top layer s is overestimated with a relative error as large as 40%.
The best predictions are provided by detector set DS1, which has all 10 detectors placed within 1l * of the source where the SRDR is most sensitive to the top layer scattering perturbation ͑see Fig. 11͒. However, the overall success in the recovery of top layer optical properties when subject to changes in scattering using this pMC/dMC is poor. Figure 13 provides the SRDR measurements for a series of phantoms in which optical scattering of the bottom layer is varied in the range 50 to 150% of the background value. For this range of scattering perturbation, an increase in bottom layer scattering results in a broad increase in the SRDR. An interesting feature is that the relative difference between the SRDR measurements is greatest at intermediate s-d separations of 1 to 3 mm. This indicates that the highest sensitivity to the scattering changes in the bottom layer of superficial Fig. 11 Experimentally measured SRDR data from extended epithelial tissue model with a scattering perturbation in the top layer. Data points shown correspond to detector set 4 ͑DS4͒. The relative error in the measured reflectance is smaller than the symbol size and less than 5% in all cases. The direction of the arrow indicates an increase of the scattering perturbation.  in recovering the bottom layer optical properties when using DS2, DS3, and DS4 detector sets for scattering perturbations in the range of 60 to 150%. The DS1 detector set does poorly and does not provide data sufficient to decouple the changes in the bottom layer s from a . Again we find that the DS3 and DS4 detector sets result in the best recovery of optical properties and the smallest confidence intervals for a . Overall the DS3 detector set provides accurate recovery of both a and s with relative error of less than 20 and 5%, respectively, for scattering perturbations in the range of 60 to 150%.

Sensitivity Characteristics of SRDR Measurements
To understand the performance characteristics of the pMC/ dMC algorithm, we conducted independent two-region MC simulations with optical properties given by Tables 1 and 2. Seven million photons were simulated for each perturbation case and the resulting SRDR is added with 2% noise to consider the potential impact of experimental error. Figures 15͑a͒-15͑d͒ present these MC simulated SRDR measurements for absorption ͓Figs. 15͑a͒ and 15͑b͔͒ and scattering ͓Figs. 15͑c͒ and 15͑d͔͒ perturbations in the top ͓Figs. 15͑a͒ and 15͑c͔͒ and bottom ͓Figs. 15͑b͒ and 15͑d͔͒ layers, respectively.  Below each SRDR plot, we also present the percent change in the SRDR for each of the perturbation cases relative to the background. These latter figures clearly show the percentage difference of the SRDR measurement relative to the background for each of the perturbation cases considered.  unperturbed case show that increases in either top or bottom layer absorption results in broad reductions in the SRDR over all s-d separations, the amount of reduction generally increasing with s-d separation. This is consistent with the inversion results shown in Figs. 8 and 10 where good inversion results were obtained for all detector sets except DS1 which has a maximum span of less than 3 mm ͑Ͻ1l * ͒. This indicates that the pMC/dMC algorithm benefits greatly from the information content of the SRDR measurement at larger s-d separations. Figures 15͑c͒ and 15͑d͒ provide the MC simulation results for scattering perturbations in the top and bottom layers, respectively. Again, the MC results for the SRDR compare favorably to the experimental measurements shown in Figs. 11 and 13, respectively. As compared to the absorption perturbation cases considered above, scattering perturbations of the top or bottom layer provide relative changes in the SRDR with a very different character. Increases in top layer scattering result in a significant increase in the SRDR only at s-d separations Շ5 mm. At larger s-d separations, the scattering perturbation does not provide much of a change in the SRDR. For s-d separations տ8 mm, increases in top layer scattering result in a slight decrease of the SRDR. By contrast, increases in bottom layer scattering result in a broad increase in the SRDR over the range of s-d separations considered.
These sensitivity features are consistent with the quality of our inversion results shown in Figs. 12 and 14. While our attempted recovery of optical properties for top layer s perturbations was not wholly successful, the best results were obtained for the DS1 detector set which has all detectors within 1l * of the source. As shown in Fig. 11, this is in fact the region where the SRDR is most sensitive to changes in top layer scattering. By contrast, our success in the recovery of bottom layer scattering perturbations is consistent with Fig.  13, which shows broad sensitivity to such perturbations over the entire range of s-d separation examined. As a result, detector sets DS3 and DS4, which span the entire range of s-d separations considered, provide the best-quality recovery of optical properties for bottom layer scattering perturbations.

Optical Property Recovery
Overall, optical property recovery using the pMC/dMC approach is best achieved using data provided by the DS3 detector configuration. This set consisted of 10 equispaced detectors located in the range of 0.5 to 9.5 mm from the source. This is equivalent to having 10 equispaced detectors spanning a distance of Ϸ3l * from the source. This detector set provided excellent optical property recovery for absorption perturbations in both the epithelial and stromal layers as well as for scattering perturbations in the stromal layer. However, optical property recovery proved difficult for scattering perturbations in the top layer. In this case, we believe that the difficulty arises not due to a limitation of the pMC/dMC technique but due to the s-d configurations employed in the experimental study. As Fig. 15͑c͒ demonstrates, sensitivity to changes in top layer scattering are concentrated at s-d separations of Շ2 mm from the source. It is quite possible that processing of SRDR data from exclusively this region would result in better optical property estimates. Moreover, as suggested by recent studies, 9-11 using a s-d detector geometry wherein the source is oriented not perpendicular to the tissue surface but angled at a 45-deg angle to the surface normal, may result in better interrogation of the epithelial layer and greater sensitivity of the SRDR measurement to scattering changes in the epithelial layer. We are currently investigating both these approaches.

Practical Considerations and Future Directions
As presented, the use of the pMC/dMC algorithm requires knowledge not only of the epithelial layer thickness but also of the optical properties of one of the two layers. Practically, the optical properties of the stromal layer can be found by processing a conventional SRDR measurement using only detector locations larger than approximately 2l * . Data from these locations are not sensitive to the optical properties of the superficial tissue layer. 19,22,23 The pMC/dMC method can then be used to determine the properties of the epithelial layer. Alternatively, one can use approaches such as elastic scattering or differential path length methods 12,13,16,17 to determine the optical properties of the epithelial tissue and use the pMC/ dMC approach to determine the optical properties of the stromal layer. From an instrumentation perspective, this latter approach is particularly attractive as it is very easy to incorporate elastic scattering and/or differential path length methodology with an SRDR measurement capability within a single fiber-based probe.
Note also that all of the optical property recovery results shown here were accomplished with only one background MC simulation. This one simulation provided the information necessary for the pMC/dMC algorithm to estimate the optical properties from 160 different SRDR measurement sets shown. Optical property recovery in each case required less than 10 iterations of the optimization algorithm and required less than 3 min on 2.1-GHz Pentium IV processor. Ultimately, we would like to move away from requiring a priori information for the epithelial layer thickness as well as the singlescattering asymmetry coefficient g, both of which we held fixed in this study. We are currently working on more generalized formulations of the pMC/dMC algorithm to make possible the recovery of these parameters, in addition to the optical properties, from a SRDR measurement.

Conclusions
We developed and validated a novel and efficient method for optical property determination of layered superficial tissues. The novelty of this method lies in the combined use of pMC and dMC methods. We used pMC to provide rapid estimates of the SRDR produced by changes in the optical properties of a heterogeneous tissue system. We used dMC to provide a sensitivity analysis of the SRDR measurement with respect to changes of the optical properties in each of the tissue regions. These two methods are used within an iterative nonlinear optimization algorithm to determine optical properties of heterogeneous structures.
We demonstrated the accuracy of this method by determining optical properties of two-layered extended epithelial tissue systems using SRDR measurements acquired using a CCD camera platform. Specifically, our results demonstrate that this combined pMC/dMC approach can successfully recover optical properties for both absorption and scattering perturba-tions in a stromal tissue layer and absorption perturbations in an epithelial tissue layer. Since MC methods are used, the inversion technique can be applied in the transport regime, that is, on spatial scales comparable to the transport mean free path and in media where optical scattering is not required to be dominant over absorption. Our results also show that the location of detectors in the SRDR measurement has a significant impact on the accuracy of the results. Thus, an opportunity exists to optimize probe designs and spatial sampling procedures of diffuse reflectance measurements for increased sensitivity to specific physiological or structural features of superficial tissues.