|
|
1.IntroductionOptical properties of biological tissues, such as absorption and scattering, can be utilized to study the morphological and biochemical status of tissues.1, 2 Among many noninvasive optical techniques, diffuse reflectance spectroscopy (DRS) has been widely applied to study the optical properties of various in vivo or ex vivo samples, for its information richness and relatively low system cost.3, 4 DRS is a model based technique whose efficiency in recovering the optical properties of samples mainly depends on the photon transport model employed. Most DRS systems employ the standard diffusion equation (SDE), which is derived from the radiative transport equation, as the forward model for its simplicity, accuracy, and computational efficiency.5 To ensure that the photon propagation can be properly described by the SDE, the ratio of μ ′ s to μ a of samples has to be much greater than unity and the photon travel length has to be beyond a certain distance [generally 5 to 10l t, where l t = 1/(μ a+ μ ′ s)].6 Therefore, the DRS techniques that work with the SDE are usually applied to study highly scattering biological deep tissues.1, 7 However, with proper modification to the SDE or the measurement geometry, the DRS techniques can work with the SDE to investigate the properties of superficial tissues. For instance, Kim proposed an semi-empirical model derived from the SDE to work with their DRS system to determine the optical properties of samples at depths in the range from 0.5 to 1.5 mm.8 In addition, we developed a DRS probe design, called the diffusing probe, which enables the SDE to work at source-detector separations shorter than 5 l t and it was used to determine the optical properties of skin in vivo.9 We have shown that the average interrogation depths of the diffusing probe are in range from 100 to 800 μm depending on the sample's optical properties.10 The diffusing probe has to work with a two-layer diffusion model that is slower than a typical semi-infinite diffusion model by three orders of magnitude.11 To enhance the model computational efficiency, we proposed a new measurement geometry, which is called the diffusing probe with planar source (DPPS) geometry as shown in Fig. 1. We have demonstrated that the computational efficiency of the diffusion model for the DPPS geometry is 1159-fold faster than the diffusion model for the diffusing probe geometry.11 The efficiency of the DPPS diffusion model is comparable to a standard diffusion model in semi-infinite geometry. In this study, we carried out a phantom study and demonstrates that the error of the DPPS method in recovering optical properties was about 6%, which is similar to the diffusing probe method.11 Since the DPPS method proposed in our previous study employs a diffusion model that assumes a semi-infinite sample geometry, it is not suitable for precisely quantifying the properties of the superficial volume of some biological samples, such as the 50- to 100-μm thick epidermis of skin. To investigate such superficial tissue volumes using the DPPS method, employing a multilayer diffusion model in the DPPS geometry is necessary. To this end, we employ a multilayer diffusion model in the DPPS geometry adopted from the classical planar source illumination (PSI) diffusion model. The PSI geometry and diffusion model were proposed by Pham to interrogate the properties of layered superficial samples.12 In this work, we use MC simulations as the reference to quantify the accuracy of the multilayer diffusion models in the DPPS and PSI geometries. We found that the PSI diffusion model does not always agree well with the Monte Carlo simulation results especially when the sample has high absorption. Based on our MC simulation data, we further explain the possible causes of inaccuracy of the diffusion model in the PSI geometry. The inverse problem, i.e., the capability of the DPPS and PSI probing geometries in recovering the optical properties of samples, will also be studied. In the simulations, the geometry and the optical properties of samples are designed to mimic in vivo human skin composed of epidermis and dermis. The recovery results for various upper layer thickness and optical properties settings will be demonstrated. Finally, the applicability of utilizing the DPPS probe to monitor the invasion depth of cutaneous melanoma will be discussed. The performance comparison between the DPPS probe and the classical PSI probe will be made in this regard. In addition, the spatial resolution of the DPPS probe is estimated using MC simulations. 2.Models and Methods2.1.Multilayer Diffusion EquationUtilizing the DPPS geometry to determine the properties of layered samples requires a proper multilayer diffusion model. Many different methods for solving the diffusion equation in layered media have been proposed.12, 13 Pham developed a diffusion model for the geometry very similar to the DPPS geometry.12 In Pham's model, a planar source irradiates onto a multilayered sample. The source and the sample have infinite lateral dimensions and the diffuse reflectance is collected at the top surface of the sample. In this study, we adopt the multi-layer diffusion model proposed by Pham and make two modifications to fit the DPPS geometry: 1. The first layer is set to be a man-made layer with high scattering and low absorption properties. 2. The detector is placed at the center of the interface between the first layer and the second layer (sample layer). Pham have provided a detailed derivation for the solution of the two-layer diffusion equation. Similarly, the solutions to the diffusion equation for media of three or more layers can be obtained. We briefly describe the derivation of the frequency domain solution of the three-layer diffusion model in the following. Note that the first layer is a man-made high scattering layer, while the second and third layers are the sample layers. Starting from the time-domain diffusion equation Eq. 1[TeX:] \documentclass[12pt]{minimal}\begin{document}\begin{eqnarray} &&\hskip-5pt\left[ {\nabla ^2 - 3\mu _{ai} \mu _{ti} - \frac{{3\mu _{ti} }}{c}\left({1 + \frac{{\mu _{ai} }}{{\mu _{ti} }}} \right)\frac{\partial }{{\partial t}} - \frac{3}{{c^2 }}\frac{{\partial ^2 }}{{\partial t^2 }}} \right]\phi _i\nonumber\\ &&\hskip-5pt\quad = \left({ - 3\mu _{ti} - \frac{3}{c}\frac{\partial }{{\partial t}}} \right)\mu '_{si} q_i, \end{eqnarray}\end{document}Eq. 2a[TeX:] \documentclass[12pt]{minimal}\begin{document}\begin{equation} q_1 (z) = P_0 \exp \left[ { - \left({\mu _{t1} + i\frac{\omega }{c}} \right)z} \right],\quad 0 < z \le \ell _1, \end{equation}\end{document}Eq. 2b[TeX:] \documentclass[12pt]{minimal}\begin{document}\begin{eqnarray} q_2 (z) &=& P_0 \exp \left[ { - \left({\mu _{t1} + i\frac{\omega }{c}} \right)\ell _1 } \right]\nonumber\\ &&\times\,\exp \left[ { - \left({\mu _{t2} + i\frac{\omega }{c}} \right)(z - \ell _1)} \right],\enspace \ell _1 < z < \ell _1 + \ell _2,\nonumber\\ \end{eqnarray}\end{document}Eq. 2c[TeX:] \documentclass[12pt]{minimal}\begin{document}\begin{eqnarray} q_3 (z) &=& P_0 \exp \left[ { - \left({\mu _{t1} + i\frac{\omega }{c}} \right)\ell _1 } \right]\exp \left[ { - \left({\mu _{t2} + i\frac{\omega }{c}} \right)\ell _2 } \right]\nonumber\\ &&\times\,\exp \left[ { - \left({\mu _{t3} + i\frac{\omega }{c}} \right)(z - \ell _1 - \ell _2)} \right],\quad \ell _1 + \ell _2 < z,\nonumber\\ \end{eqnarray}\end{document}3.Eq. 3a[TeX:] \documentclass[12pt]{minimal}\begin{document}\begin{eqnarray} \phi _1 (z) &=& \psi _1 \exp \left[ { - \left({\mu _{t1} + i\frac{\omega }{c}} \right)z} \right] + A_1 \exp \left[ { - \frac{z}{{\delta _1 }}} \right]\nonumber\\ && +\, A_2 \exp \left[ {\frac{z}{{\delta _1 }}} \right],\quad 0 < z \le \ell _1, \end{eqnarray}\end{document}Eq. 3b[TeX:] \documentclass[12pt]{minimal}\begin{document}\begin{eqnarray} \phi _2 (z) &=& \psi _2 \exp \left[ { - \left({\mu _{t2} + i\frac{\omega }{c}} \right)(z - \ell _1)} \right] + A_3 \exp \left[ { - \frac{z}{{\delta _2 }}} \right]\nonumber\\ && +\, A_4 \exp \left[ {\frac{z}{{\delta _2 }}} \right],\quad \ell _1 < z < \ell _1 + \ell _2, \end{eqnarray}\end{document}Eq. 3c[TeX:] \documentclass[12pt]{minimal}\begin{document}\begin{eqnarray} \phi _3 (z) &=& \psi _3 \exp \left[ { - \left({\mu _{t3} + i\frac{\omega }{c}} \right)(z - \ell _1 - \ell _2)} \right]\nonumber\\ && +\, A_5 \exp \left[ { - \frac{z}{{\delta _3 }}} \right],\quad \ell _1 + \ell _2 < z, \end{eqnarray}\end{document}Eq. 4[TeX:] \documentclass[12pt]{minimal}\begin{document}\begin{equation} \delta _i = \left[ {3\mu _{ai} \mu _{ti} - \frac{{3\omega ^2 }}{{c^2 }} + i\left({1 + \frac{{\mu _{ai} }}{{\mu _{ti} }}} \right) \frac{{3\mu _{ti} \omega }}{c}} \right]^{ - 1/2}, \end{equation}\end{document}Eq. 5a[TeX:] \documentclass[12pt]{minimal}\begin{document}\begin{equation} \psi _1 = \frac{{3\delta _1 ^2 \mu _{s1} \mu _{t1} \big(1 + i\frac{\omega }{{c\mu _{t1} }}\big)}}{{1 - \delta _1 ^2 \big(\mu _{t1} + i\frac{\omega }{c}\big)^2 }}, \end{equation}\end{document}Eq. 5b[TeX:] \documentclass[12pt]{minimal}\begin{document}\begin{equation} \psi _2 = \frac{{3\delta _2 ^2 \mu _{s2} \mu _{t2} \big(1 + i\frac{\omega }{{c\mu _{t2} }}\big)}}{{1 - \delta _2 ^2 \big(\mu _{t2} + i\frac{\omega }{c}\big)^2 }}\exp \left[ { - \left({\mu _{t1} + i\frac{\omega }{c}} \right)\ell _1 } \right], \end{equation}\end{document}Eq. 5c[TeX:] \documentclass[12pt]{minimal}\begin{document}\begin{eqnarray} \psi _3 &=& \frac{{3\delta _3 ^2 \mu _{s3} \mu _{t3} \big(1 + i\frac{\omega }{{c\mu _{t3} }}\big)}}{{1 - \delta _3 ^2 \big(\mu _{t3} + i\frac{\omega }{c}\big)^2 }}\nonumber\\ &&\times\,\exp \left[ { - \left({\mu _{t1} + i\frac{\omega }{c}} \right)\ell _1 - \left({\mu _{t2} + i\frac{\omega }{c}} \right)\ell _2 } \right]. \end{eqnarray}\end{document}6.Eq. 6a[TeX:] \documentclass[12pt]{minimal}\begin{document}\begin{equation} j_1 (0) = \frac{{1 - R_{\it eff} }}{{2(1 + R_{\it eff})}}\phi _1 (0), \end{equation}\vspace*{-6pt}\end{document}Eq. 6b[TeX:] \documentclass[12pt]{minimal}\begin{document}\begin{equation} \frac{{\phi _1 (\ell _1)}}{{\phi _2 (\ell _1)}} = \frac{{n_1 ^2 }}{{n_2 ^2 }}, \end{equation}\vspace*{-6pt}\end{document}Eq. 6c[TeX:] \documentclass[12pt]{minimal}\begin{document}\begin{equation} \frac{{\phi _2 (\ell _1 + \ell _2)}}{{\phi _3 (\ell _1 + \ell _2)}} = \frac{{n_2 ^2 }}{{n_3 ^2 }}, \end{equation}\vspace*{-10pt}\end{document}2.2.Monte Carlo SimulationThe MC model is employed here as a benchmark method to evaluate the performance of the multilayer diffusion models, either in the DPPS or PSI geometries. Our MC computer code was adopted from the one developed by Wang 14 In the MC model, the source and the sample are defined in the cylindrical coordinate. The detector is located at the origin of the system and its diameter and refractive index are set to 100 μm and 1.4, respectively. To reduce the simulation time, an acceptance angle of 180° is used for the detector. The radius of the sample is 106 cm to simulate infinite lateral sample dimension while the circular source is uniformly distributed within the radius r s. Photon packets propagating in the media according to the Henyey-Greenstein phase function with the anisotropy factor g = 0.8. Each simulation continuously ran until a total number of 10,000 photon packets arrived at the detector. In most cases in this study, this program termination condition required at least 108 photon packets being launched into the system. For the DPPS geometry, a high scattering layer of n = 1.35, μ a = 10−6/mm, μ ′ s = 20/mm, and thickness = 0.5 mm was used in the simulations.11 Refractive indices of the sample layers were kept at 1.4. In all MC simulations, the frequency domain reflectance was generated in the range from 0 to 1500 MHz. As the MC simulated frequency domain reflectance was employed in the inverse fitting, random noise of 0.3% in amplitude demodulation, and 0.08° in phase delay were added to the data to simulate measurement noise. 2.3.InstrumentationA frequency domain diffuse reflectance measurement system is employed to understand the performance of the two measurement geometries discussed in this study. Similar systems have been used to quantify the optical properties of various samples.15, 16 Our system is composed of a 658-nm laser diode, a laser diode mount (TCLDM9, THORLABS, New Jersey), a laser diode controller (LDC-3908, ILX Lightwave, Montana), a network analyzer (N5230C, Agilent Technologies, California), and an avalanche photodiode (C5658, Hamamatsu, New Jersey). We fabricated two optical probes that had PSI and DPPS geometries, respectively. A circular planar light of 10 mm radius was created as the light source for both geometries. A 1-mm thick Spectralon slab (LabSphere, New Hampshire) (μ a = 10−6/mm, μ ′ s = 50/mm @658 nm) was used as the high scattering layer in the DPPS geometry. The optical fibers used to fabricate the optical probes had a core diameter of 1 mm and numerical aperture of 0.48. The laser diode was modulated from 50 to 450 MHz and the average output power measured at the distal end of the optical probe was 7 mW. The optical probes were placed on the surface of a solid phantom (μ a = 0.081/mm, μ ′ s = 2.532/mm @ 658 nm), which was made by mixing silicone, Titanium dioxide (scatterer), and India ink (absorber). The radius and thickness of the phantom were 5 and 8 cm, respectively. In the measurement, the diffuse reflectance was detected by the avalanche photodiode, and the network analyzer received the signal from the avalanche photodiode and determined the phase delay and the amplitude demodulation versus the modulation frequency. The data were fit to an appropriate photon diffusion model using the MATLAB least square fitting method “lsqcurvefit” (MathWorks, Massachusetts) to calculate absorption and reduced scattering coefficients of the phantom. 3.Results and Discussion3.1.Analysis of the Diffusion Model in the Planar Source Illumination GeometryTo understand the performance of this diffusion model, we first study a simple case where the sample is semi-infinite and the planar light source is directly irradiated on the sample. In such a case, the optical properties of the second and third layers in three-layer diffusion model derived in Sec. 2 are set the same and the diffuse reflectance is calculated at the sample's surface. The frequency domain reflectance determined using the MC and diffusion models in the PSI geometry are shown in Fig. 2 as symbols and lines, respectively. To observe the effect of the sample optical properties on the reflectance, two sets of sample optical properties were employed to simulate light (μ a = 0.03/mm, μ ′ s = 3.0/mm) and dark (μ a = 0.3/mm, μ ′ s = 3.0/mm,) complexion human skin.4, 17 In Fig. 2, the maximum percent deviations of amplitude demodulation and phase delay between the diffusion and MC models are 0.003%, −7.45% and 0.03%, −16.73% for light skin and dark skin, respectively. It can be seen that the discrepancy between the diffusion and MC models increases as the sample absorption coefficient increases. This phenomenon is expected since the diffusion approximation becomes more tenuous when the μ ′ s to μ a ratio gets smaller. Note that the phase delay of the diffusion model is always smaller than that of the MC model. This can be understood by observing the time-resolved reflectance simulation results demonstrated by various independent groups.18, 19 The time resolved diffusion model is not able to accurately predict the photon behavior at a very short time frame. In general, the standard diffusion model predicts more short time-of-flight of photons than the Monte Carlo model does.18, 19 Thus, the average travel time of photons estimated by the diffusion model is shorter than that estimated by the MC model. The shorter travel time in the time domain translates to a smaller phase delay in the frequency domain. Therefore, the phase delay of the frequency domain reflectance determined using the diffusion model is smaller than that determined using the MC model in the PSI geometry. Fig. 2(a) Amplitude demodulation and (b) phase delay versus source modulation. The absorption and reduced scattering coefficients of the semi-infinite light and dark skin are μ a = 0.03/mm, μ ′ s = 3.0/mm and μ a = 0.3/mm, μ ′ s = 3.0/mm, respectively. Lines and symbols represent reflectance determined using the diffusion model and the MC model, respectively.  To study the inverse problem for the sample optical properties recovery, we used the diffusion model as a forward model to fit the MC model generated reflectance. The MATLAB least square curve fitting routine lsqcurvefit was employed to perform the data fitting. The recovery results and percent errors are listed in Table 1. In summary, the recovery errors are within 26.5% for light skin and increase to within 84.8% for dark skin. The increase in the recovery error as the sample's absorption increases is not surprising since we have seen this trend in the raw data illustrated in Fig. 2. The results shown here suggest that the PSI measurement geometry with the diffusion model cannot reliably determine the optical properties of samples. To further investigate the reasons that induce the results reported above, we employed MC simulations to look into the photon behavior in the PSI geometry. Fig. 8(a) Recovered upper layer thickness versus the tumor radius used in the Monte Carlo simulations. Squares represent the tumor radii used in the MC simulations. Circles and triangles represent the upper layer thicknesses recovered in the DPPS and PSI geometries, respectively. (b) Monte Carlo simulated amplitude demodulation (squares) and phase delay (triangles) versus lateral offset between a tumor with 0.5 mm radius and the detector.  Table 1Optical properties of semi-infinite samples recovered with the PSI and the DPPS geometries.
We recorded the incident positions of all detected photon packets in the MC simulation and calculated the detected photon packet number versus the incident radial position. Figure 3 shows the normalized detected photon packets number versus the incident radial position for a light skin sample. In the MC simulation, a total amount of ten thousand photon packets was captured at the detector. The detected photon packet number originating from a certain radial position was normalized with respect to the total amount of captured photon packets. Since the planar source in the PSI diffusion model is infinitely wide, the source radius r s in the MC model has to be large enough to satisfy the infinity assumption. From our MC simulation results, the source radius r s has to be at least 10 mm to produce a stable detected reflectance. The difference between the reflectance determined at r s = 10 and 20 mm is within 0.3%. Fig. 3Normalized detected photon number versus incident radial position in a semi-infinite medium with μ a = 0.03/mm and μ ′ s = 3.0/mm. The cumulative distribution function of the normalized detected photon number is also depicted (dotted line).  In addition, Fig. 3 illustrates the cumulative distribution function (CDF) for the normalized detected photon packet number. From the CDF, it can be observed that about 50% of detected photon packets are launched from the radial positions within 1.2 mm (3.64 l t) for a light skin sample (as shown in Fig. 3) and 0.9 mm (2.97 l t) for a dark skin sample (data not shown). These results indicated that the PSI geometry is analogous to a classical point source, point detector DRS measurement scheme at a very short source-detector separation. Furthermore, our MC code is capable of calculating the cumulative zigzag path length and the average interrogation depth of photon packets in the PSI measurement geometry.10 The average cumulative travel length and the average interrogation depth were 3.87 mm (11.73 l t) and 424 μm for light skin, and 1.37 mm (4.52 l t) and 204 μm dark skin samples, respectively. We can observe that although the PSI geometry employ an infinitely wide planar light source, the photons detected in this geometry mostly originate from the region near the detector and their travel lengths are short. Since the standard diffusion equation generally does not work well under these conditions, the performance of the PSI method is poor as demonstrated in Table 1. 3.2.Analysis of the Diffusion Model in the Diffusing Probe with Planar Source GeometryThe accuracy of the diffusion model is greatly improved when it is used in the DPPS geometry. The frequency domain diffuse reflectance in the DPPS geometry determined using the diffusion and MC models are illustrated in Fig. 4. In both models, the first layer is a high scattering layer and the second layer is a semi-infinite sample layer that has optical properties of light or dark skin as used in Sec. 3.1. The reflectance is collected at the center of the boundary between the first and second layers. As shown in Fig. 4, the frequency domain reflectance in the DPPS geometry determined using the diffusion and MC models have excellent agreement. The differences between the reflectance generated using the two models, either for light or dark skin, are within 0.58 and 2.04% for the amplitude demodulation and the phase delay, respectively. Note that the refractive index of the detector might affect the measurement as indicated by Bargo 20 As we modified the refractive index of the detector 1.4 to 1.458 in the MC simulations, the determined reflectance did not show any differences. This infers that the presence of the detector fiber in a real measurement system has minimal influence on the reflectance and, therefore, boundary conditions (6a)–(6e) are reasonable. Using the diffusion model to fit the MC generated data, the recovered optical properties of the samples are listed in Table 1. The DPPS geometry can accurately recover the sample optical properties with percent errors smaller than those obtained in the PSI geometry generally by an order of magnitude. Fig. 4(a) Amplitude demodulation and (b) phase delay versus source modulation frequency in the DPPS geometry with the first layer optical properties of μ a = 10−6/mm, μ ′ s = 20/mm. The absorption and reduced scattering coefficients of the semi-infinite light and dark skin samples are μ a = 0.03/mm, μ ′ s = 3.0/mm and μ a = 0.3/mm, μ ′ s = 3.0/mm, respectively. Lines and symbols represent reflectance determined using the diffusion model and the MC model, respectively.  It is worth noting that the reflectance for light skin sample shown in Fig. 4 is almost identical to the reflectance illustrated in our previous work [Fig. 2a in Ref. 11]. The light source employed in our previous work was a delta function located in the sample, while an exponentially decaying light source from the sample surface is assumed in this work. Almost identical results demonstrated in this work and our previous work using the same sample optical parameters infers that the light source setting would not introduce significant differences in the determined reflectance as the diffusion models are used in the DPPS geometry. The excellent agreement between the DPPS diffusion and MC models comes from the great number of scattering events introduced by the highly scattering slab. The long travel length of photons in the high scattering layer facilitates larger light fluence ϕ to flux j ratio, and, thus, diffusion approximation is better satisfied. This phenomenon has been studied in detail in our previous work.10 The normalized detected photon number and CDF versus the incident radial position is illustrated in Fig. 5. In Fig. 5, the radii for the CDF to reach 50 and 90% are 2.7 and 8.8 mm, respectively. For dark skin, the radii for the CDF to reach 50 and 90% are 1.3 and 3.4 mm, respectively (data not shown). The average cumulative travel lengths and the average interrogation depths were 9.69 mm (29.36 l t) and 460 μm for light skin, and 1.52 mm (5.04 l t) and 162 μm for dark skin, respectively. In addition, the average cumulative travel lengths of photon packets in the high scattering layer were 4.04 mm (80.8 l t) and 3.11 mm (62.2 l t) for light and dark skin samples, respectively. Fig. 5Normalized detected photon number and CDF versus incident radial position in a semi-infinite medium with μ a = 0.03/mm and μ ′ s = 3.0/mm.  It can be seen that the photon travel lengths in the sample are longer in the DPPS geometry than in the PSI geometry. The longer sample travel length in the DPPS geometry results from the fact that photons in the DPPS geometry generally have larger incident radii than in the PSI geometry, as shown in Figs. 3, 5. In addition to the long travel length of photons provided by the high scattering layer in the DPPS geometry, the longer sample travel length of photons in the DPPS geometry further ensures the validation of diffusion approximation. Although the average sample travel length is longer in the DPPS geometry than in the PSI geometry, the average interrogation depths of the two geometries are comparable. Our MC simulations indicated that most photons have oblique sample incident angle in the DPPS geometry. The interrogation depths of oblique incident photons are usually shallower than the interrogation depths of the normal incident photons, such as the photons in the PSI geometry, at a same incident radius. The interplay between the facts that the photons have larger incident radii and oblique incident angles in the DPPS geometry than in the PSI geometry results in similar interrogation depths in both geometries. The difference between the average cumulative path lengths for photon packets traveling in light and dark skin is 24.32 l t for the DPPS and 7.21 l t for the PSI geometries. Consequently, the difference between the amplitude demodulation or the phase delay of light and dark skin is larger in the DPPS geometry than in the PSI geometry, as can be seen by comparing Figs. 2, 4. This observation infers that the sensitivity of the DPPS method in the variation of sample optical properties is not impaired by the presence of the high scattering layer; rather, it is more sensitive than the PSI method. Comparing Figs. 3, 5, the slope of the CDF is smaller in the DPPS geometry than in the PSI geometry. The slope of the CDF affects the actual size of the probe. To rigorously investigate the effect of the finite probe radius on the detected reflectance and determine the minimal probe size of the DPPS geometry, we carried out two additional MC simulations in which the probe radii were 1 and 0.5 cm, respectively. The results are displayed in Fig. 6. The reflectance data for the light skin sample illustrated in Fig. 4 are also included in Fig. 6 for comparison (r s = 2 cm). It can be seen that as the radius of the probe decreases, the deviations between the results determined from the MC and diffusion models become more substantial. Quantitatively, the maximum deviations of the amplitude demodulation and phase delay generated using the MC and diffusion models for r s = 1 and 0.5 cm are 0.7%, 1.2% and 6.4%, 6.6%, respectively. Using the DPPS diffusion model to fit the MC generated data, the recovered sample optical properties are μ a = 0.03/mm, μ ′ s = 3.261/mm, for r s = 1 cm, and μ a = 0.043/mm, μ ′ s = 1.803/mm, for r s = 0.5 cm. The percent errors of the revered absorption and scattering coefficients are 3.1%, 8.7%, for r s = 1 cm, and 43.3%, 39.9%, for r s = 0.5 cm, respectively. At r s = 0.5 or 1 cm, the recovery results for dark skin are more accurate than those for light skin due to smaller photon incident radii for dark skin as discussed earlier. Fig. 6(a) Amplitude demodulation and (b) phase delay versus source modulation frequency in the DPPS geometry with the first layer optical properties of μ a = 10−6/mm, μ ′ s = 20/mm and the second layer optical properties of μ a = 0.03/mm, μ ′ s = 3/mm. Solid lines and symbols represent the reflectance determined using a two-layer diffusion model and the MC model, respectively. Squares, circles, and triangles represent MC results for the probe radius of 2, 1, and 0.5 cm, respectively.  The simulation results shown here suggest that when measuring light skin samples, the probe radius of 1 cm is sufficient to meet the infinitely wide source assumption in the diffusion model within 1.2% error. Since the minimal DPPS probe radius decreases as the absorption coefficient of the sample increases, we estimate that the minimal probe radius for investigating all skin phototypes is about 1 cm. Furthermore, our simulation results reveal that the validity of the DPPS diffusion model holds for a wide range of optical properties. The DPPS diffusion model generally deviates from the MC model (with a probe radius of 1 cm) less than 3% for samples that possess biological tissue relevant optical properties. For example, we used μ a = 0.04/mm, μ ′ s = 0.36/mm to simulate cervix tissue,21 the deviations between the reflectance generated using the DPPS diffusion and MC models with a probe radius of 1 cm were 0.29 and 1.98% for amplitude demodulation and phase delay, respectively. Fitting the MC data with the DPPS diffusion model, the recovered optical properties had errors within 4.3%. We carried out a preliminary phantom measurement to study the performance of the DPPS geometry. The DPPS probe was gently placed on the surface of the phantom. The frequency domain diffuse reflectance in terms of amplitude demodulation and phase delay are demonstrated in Fig. 7. In Fig. 7, the filled squares and the corresponding lines represent the measurement data and the fitting results, respectively, for the DPPS geometry. The recovered optical properties of the phantom were μ a = 0.076/mm and μ ′ s = 2.364/mm, representing −6.2% and −6.6% deviations from the reference optical properties of the phantom. Moreover, the phantom measurement results for the PSI geometry are also illustrated in Fig. 7 as open squares. The recovery results for the PSI geometry were μ a = 0.040/mm and μ ′ s = 1.281/mm, representing −50.6% and −49.4% deviations from the reference values. The percent errors of the recovery results of the DPPS and PSI geometries are comparable to the simulations results listed in Table 1. 3.3.Extraction of the Properties of Multilayer SamplesIn Secs. 3.1, 3.2, we have demonstrated that the DPPS method is more accurate than the PSI method in determining the optical properties of semi-infinite samples. In this section, the ability of the DPPS and PSI methods in recovering the properties of multilayer samples will be studied and compared. MC simulations were employed to generate the reference reflectance, which was then fit by the DPPS and PSI diffusion models to recover the sample parameters. In the MC simulations, the samples were designed to simulate human skin. The epidermis and dermis were approximated as a plane-parallel slab of variable thickness and a semi-infinite medium, respectively. Thus, there are five variables in the diffusion model to characterize the sample: absorption and scattering coefficients of the upper and lower layers, and the thickness of the upper layer of the sample. Pham indicated that the frequency domain diffuse reflectance could be used to robustly derive two out of five variables of such two-layer samples.12 They also pointed out that fitting more than two variables often resulted in convergence failures and the outcome was very sensitive to the initial guesses. We have reached the same conclusion with our simulation results, either for the PSI or the DPPS geometries. Therefore, in this section, we will only show the results of employing the diffusion models to fit the MC data to recover the following three sets of variables: 1. absorption and scattering coefficients of the upper sample layer, 2. absorption and scattering coefficients of the lower sample layer, and 3. absorption coefficients of the upper and lower sample layers. In recovering the values of the two variables of each set, other three variables are assumed to be known a priori. The optical properties of the human epidermis and dermis in the 600 to 1000 nm wavelength range are primarily induced by melanin, collagen, blood, and water, and have been reported by various groups.22, 23, 24 We referred to the optical properties of epidermis and dermis reported by other groups and chose μ a = 0.5/mm, μ ′ s = 3.0/mm and μ a = 0.03/mm, μ ′ s = 3.0/mm to simulate the absorption and the scattering coefficients of epidermis and dermis of light skin at 650 nm, respectively. The thickness of the epidermis was 50 or 100 μm in the simulations. The results of the optical property recovery are listed in Tables 2, 3. Table 2Optical properties of two-layer samples recovered with the PSI and the DPPS geometries.
Table 3Absorption coefficients of two-layer samples recovered with the PSI and the DPPS geometries.
It can be seen in Table 2 that the epidermal optical properties recovered using the PSI method have errors that are much larger than those recovered with the DPPS method for either epidermal thickness. Note that the epidermal optical property recovery errors of both methods reduce as the epidermis thickness increases from 50 to 100 μm. As the epidermis thickness increases, the detected photon packets travel longer in the epidermis and the variation of the epidermal properties affect more on the reflectance. Therefore, the recovered epidermal optical properties are more faithful as the thickness of the epidermis increases. Moreover, as shown in Table 2, the DPPS method is more accurate than the PSI method not only for epidermal interrogation but also for the determination of the dermal optical properties. However, from the third and sixth columns of Table 2, it is observed that, for the DPPS geometry, the recovery errors increase as the target of optical property recovery changed from epidermis to dermis. Conversely, for the PSI geometry, the recovery errors decrease as the recovery target changed from epidermis to dermis. It can be inferred that the DPPS method is very sensitive to the superficial volumes and suitable for epidermis interrogation applications. Based on the simulation results, we have found that, in order to accurately recover the epidermal optical properties in vivo using the DPPS method, the epidermal thickness and the dermal optical properties have to be determined in advance. Practically, the thickness of the epidermis of a certain region can be noninvasively determined using alternative techniques such as ultrasonography or optical coherence tomography (OCT). In addition, the dermal optical properties have been summarized by several groups25, 26 and the values generally do not significantly vary between normal subjects. Thus, by keeping the dermal optical properties and the epidermis thickness as known parameters, the DPPS method could be used to evaluate and monitor the epidermal optical properties in real time. Furthermore, van Gemert indicated that the scattering coefficients of the epidermis and the dermis are very similar in the visible-NIR range.26 By taking the thickness of the epidermis (50 or 100 μm) and the scattering coefficients of the epidermis and the dermis (μ ′ s = 3.0/mm) as known parameters, we recovered the absorption coefficients of the lower and upper layers of the samples and the results are listed in Table 3. It is worth noting that the recovery errors of the epidermal absorption coefficient are in the range from 9.1 to 17.4% for the DPPS method, while those for the PSI method are in the range from 77.5 to 99.9%. The dermis absorption coefficient recovery errors for both methods are very similar. As the scattering coefficients of the epidermis and dermis are known, the DPPS method greatly excels the PSI method in recovering the epidermis absorption coefficient. The pigment in the epidermis and the vascular network in the upper dermis affect the complexion of skin. Simultaneously, monitoring UV-induced pigmentation and vascular changes remains a challenging task. Our results suggest that the DPPS method could be used to study the variation of the epidermal pigmentation and the dermal hemodynamics of in vivo human skin introduced by UV exposure. 3.4.Estimation of the Cutaneous Tumor Invasion DepthDetermining and monitoring the depth of the malignant cells invasion is crucial for skin melanoma prognosis. The depth that the tumor cells have invaded has been an effective prognostic factor for skin melanoma diagnoses. As the tumor depth (Breslow's depth) is less than 2 mm, the 5-year survival rate is approximately 80 to 96%.27 Thus, being able to detect the tumor growth as its size is smaller than 2 mm is critical for skin tumor prognosis. In this section, we employed a series of MC simulations to investigate the feasibility of the DPPS and PSI methods for cutaneous melanoma prognostic applications. In the MC simulations, we placed a spherical tumor of radius in the range from 0 to 1.5 mm and μ a = 0.5/mm, μ ′ s = 3/mm in a semi-infinite sample of light skin optical properties (μ a = 0.03/mm, μ ′ s = 3/mm) to simulate skin melanoma. The tumor surface was flush with the detector and the center of the tumor was vertically aligned with the detector. The MC simulations were carried out in the DPPS and PSI geometries with a 1-cm probe radius. The MC simulated data were then fit to the DPPS or PSI diffusion models to calculate the thickness of the homogeneous, infinitely wide upper sample layer. In the diffusion models, the optical properties of the upper and lower layers of the sample were kept constant at μ a1 = 0.5, μ ′ s1 = 3/mm and μ a2 = 0.03/mm, μ ′ s2 = 3/mm, respectively. The determined thickness of the upper layer of the sample is illustrated in Fig. 8a. It can be seen that the trend of the upper layer thickness recovered using DPPS method follows well with the true tumor radius variation. Due the limitation of the current forward model, we cannot determine the exact size of the tumor since the tumor is treated as a slab in our forward model. However, our simulation results, shown in Fig. 8, suggest that the DPPS method is effective in monitoring the growth of the spherical tumor of radius as small as 0.25 mm, which is critical for skin tumor prognosis. In contrast, the PSI geometry is not sensitive to the presence of tumor of radius smaller than 0.25 mm. Moreover, with tumor radius smaller than 0.75 mm, the upper layer thickness recovered in the PSI geometry does not share the same trend as the variation of the true tumor radius. Please note that Fig. 8a shows that sometimes the PSI method determines the “thickness of an infinitely wide slab-shape tumor” that is closer to the “radius of tumor” than the DPPS method does. This does not mean that the PSI method is more accurate than the DPPS method, since the thickness of a slab cannot be directly translated to the radius of a sphere. It is important to observe in Fig. 8a that the trend of the recovery results of the DPPS method is similar to the trend of tumor radius variation (or tumor growth rate), while the trends of the results of the PSI and DPPS methods are distinct. This observation suggests that the DPPS method could be useful for early tumor development monitoring. The setting of the upper layer optical properties in the DPPS diffusion model is flexible. We have tested many upper layer setting combinations in the diffusion model that have μ a higher than 0.03/mm and arbitrary μ ′ s, and have found similar trends as demonstrated in Fig. 8a.11 Although the current forward model does not allow us to determine the exact diameter of a tumor, our preliminary simulation results suggest that the DPPS probe could be an effective and efficient tool to monitor the vertical progression of melanotic lesions in vivo and could facilitate the prognosis of skin melanoma. Moreover, the spatial resolution of the DPPS probe was characterized by another set of MC simulations. In this set of MC simulations, the radius of the tumor (μ a = 0.5/mm, μ ′ s = 3/mm) was fixed at 0.5 mm (volume = 0.52 mm3) and the tumor surface was flush with the semi-infinite sample's surface (μ a = 0.03/mm, μ ′ s = 3/mm). The probe radius was 10 mm and the offset between the centers of the tumor and the detector was adjusted from 0 to 10 mm. The amplitude demodulation and phase delay of the frequency domain reflectance at 1500 MHz versus the tumor lateral offset are illustrated in Fig. 8b. It can be seen that the DPPS probe is generally not sensitive to the inhomogeneity located more than 4 mm apart from the detector. The tumor offset that introduces half maximum of the phase delay curve in Fig. 8b is approximately 1.5 mm. 4.ConclusionIn our previous study, we have shown that the DPPS geometry with its accompanying diffusion model could be used to determine the optical properties of superficial tissues in real time. In this study, we employed the MC simulation results as the reference to characterize the performance of the DPPS and the PSI methods in detail. First, we used semi-infinite sample geometry to investigate the accuracy of the DPPS and PSI diffusion models. We found that both the PSI and DPPS geometries had shallow interrogation depths; however, while the DPPS diffusion model could determine the phase delay of the frequency domain reflectance within 2.04% error for the sample with dark skin optical properties, the PSI diffusion model deviated from the MC model by 16.73% under the same conditions. Therefore, the recovery results of the two measurement geometries showed notable differences: while the DPPS geometry demonstrated optical property recovery errors within 6.43% for the dark skin sample, the optical property recovery errors for the PSI geometry were within 84.77% for the same setting. We have carried out preliminary phantom measurements to support our simulation results. Second, we employed a two-layer sample geometry to mimic the epidermis and dermis. We found that using the DPPS geometry in conjunction with the diffusion model could precisely determine the optical properties of 50-μm thick epidermis with maximal error of 14.7% as the thickness of epidermis and the optical properties of the dermis were known. The recovery error was reduced to 10.8% as the epidermis thickness is increased to 100 μm. The PSI method, on the contrary, showed very large recovery errors and was, therefore, not suitable for epidermis interrogation. Our simulation results also indicated that the DPPS method is able to simultaneously determine the epidermal and,dermal absorption coefficients. The DPPS probe could be a very useful tool to study the epidermal (ex: skin immediate pigment darkening) and the dermal (ex: erythema) reaction to UV radiation. Finally, to study the applicability of the DPPS method for skin melanoma prognosis and monitoring, we employed MC simulations and found that the DPPS method could discern the presence of a cutaneous tumor as the radius of the tumor is as small as 0.25 mm and has a lateral resolution of 1.5 mm. AcknowledgmentsSheng-Hao Tseng would like to acknowledge the support provided by the National Science Council of Taiwan under Grant No. NSC-99-2628-E-006-012. ReferencesA. Cerussi, N. Shah, D. Hsiang, A. Durkin, J. Butler, and
B. J. Tromberg,
“In vivo absorption, scattering, and physiologic properties of 58 malignant breast tumors determined by broadband diffuse optical spectroscopy,”
J. Biomed. Opt., 11 044005
(2006). https://doi.org/10.1117/1.2337546 Google Scholar
H. Ou-Yang, G. Stamatas, and
N. Kollias,
“Spectral responses of melanin to ultraviolet A irradiation,”
J. Invest. Dermatol., 122 492
–496
(2004). https://doi.org/10.1046/j.0022-202X.2004.22247.x Google Scholar
R. Zhang, W. Verkruysse, B. Choi, J. A. Viator, R. Jung, L. O. Svaasand, G. Aguilar, and
J. S. Nelson,
“Determination of human skin optical properties from spectrophotometric measurements based on optimization by genetic algorithms,”
J. Biomed. Opt., 10 024030
(2005). https://doi.org/10.1117/1.1891147 Google Scholar
S. H. Tseng, P. Bargo, A. Durkin, and
N. Kollias,
“Chromophore concentrations, absorption and scattering properties of human skin in-vivo,”
Opt. Express, 17 14599
–14617
(2009). https://doi.org/10.1364/OE.17.014599 Google Scholar
R. C. Haskell, L. O. Svaasand, T. T. Tsay, T. C. Feng, and
M. S. Mcadams,
“Boundary-conditions for the diffusion equation in radiative-transfer,”
J. Opt. Soc. Am. A, 11 2727
–2741
(1994). https://doi.org/10.1364/JOSAA.11.002727 Google Scholar
K. M. Yoo, F. Liu, and
R. R. Alfano,
“When does the diffusion-approximation fail to describe photon transport in random-media,”
Phys. Rev. Lett., 64 2647
–2650
(1990). https://doi.org/10.1103/PhysRevLett.64.2647 Google Scholar
E. Gratton, V. Toronov, U. Wolf, M. Wolf, and
A. Webb,
“Measurement of brain activity by near-infrared light,”
J. Biomed. Opt., 10 011008
(2005). https://doi.org/10.1117/1.1854673 Google Scholar
A. Kim, M. Roy, F. Dadani, and
B. C. Wilson,
“A fiber optic reflectance probe with multiple source-collector separations to increase the dynamic range of derived tissue optical absorption and scattering coefficients,”
Opt. Express, 18 5580
–5594
(2010). https://doi.org/10.1364/OE.18.005580 Google Scholar
S. H. Tseng, C. Hayakawa, B. J. Tromberg, J. Spanier, and
A. J. Durkin,
“Quantitative spectroscopy of superficial turbid media,”
Opt. Lett., 30 3165
–3167
(2005). https://doi.org/10.1364/OL.30.003165 Google Scholar
S. H. Tseng, C. Hayakawa, J. Spanier, and
A. J. Durkin,
“Investigation of a probe design for facilitating the uses of the standard photon diffusion equation at short source-detector separations: Monte Carlo simulations,”
J. Biomed. Opt., 14 054043
(2009). https://doi.org/10.1117/1.3253386 Google Scholar
S. H. Tseng and
M. F. Hou,
“Analysis of a diffusion-model-based approach for efficient quantification of superficial tissue properties,”
Opt. Lett., 35 3739
–3741
(2010). https://doi.org/10.1364/OL.35.003739 Google Scholar
T. H. Pham, T. Spott, L. O. Svaasand, and
B. J. Tromberg,
“Quantifying the properties of two-layer turbid media with frequency-domain diffuse reflectance,”
Appl. Opt., 39 4733
–4745
(2000). https://doi.org/10.1364/AO.39.004733 Google Scholar
A. Liemert and
A. Kienle,
“Light diffusion in N-layered turbid media: frequency and time domains,”
J. Biomed. Opt., 15 025002
(2010). https://doi.org/10.1117/1.3368682 Google Scholar
L. H. Wang, S. L. Jacques, and
L. Q. Zheng,
“Mcml – Monte-Carlo modeling of light transport in multilayered tissues,”
Comput. Methods Programs Biomed., 47 131
–146
(1995). https://doi.org/10.1016/0169-2607(95)01640-F Google Scholar
J. B. Fishkin, S. Fantini, M. J. VandeVen, and
E. Gratton,
“Gigahertz photon density waves in a turbid medium: theory and experiments,”
Phys. Rev. E, 53 2307
–2319
(1996). https://doi.org/10.1103/PhysRevE.53.2307 Google Scholar
T. H. Pham, O. Coquoz, J. B. Fishkin, E. Anderson, and
B. J. Tromberg,
“Broad bandwidth frequency domain instrument for quantitative tissue optical spectroscopy,”
Rev. Sci. Instrum., 71 2500
–2513
(2000). https://doi.org/10.1063/1.1150665 Google Scholar
S. H. Tseng, A. Grant, and
A. J. Durkin,
“In vivo determination of skin near-infrared optical properties using diffuse optical spectroscopy,”
J. Biomed. Opt., 13 014016
(2008). https://doi.org/10.1117/1.2829772 Google Scholar
E. Alerstam, S. Andersson-Engels, and
T. Svensson,
“White Monte Carlo for time-resolved photon migration,”
J. Biomed. Opt., 13 041304
(2008). https://doi.org/10.1117/1.2950319 Google Scholar
A. Kienle and
M. S. Patterson,
“Improved solutions of the steady-state and the time-resolved diffusion equations for reflectance from a semi-infinite turbid medium,”
J. Opt. Soc. Am. A, 14 246
–254
(1997). https://doi.org/10.1364/JOSAA.14.000246 Google Scholar
P. R. Bargo, S. A. Prahl, and
S. L. Jacques,
“Collection efficiency of a single optical fiber in turbid media,”
Appl. Opt., 42 3187
–3197
(2003). https://doi.org/10.1364/AO.42.003187 Google Scholar
R. Hornung, T. H. Pham, K. A. Keefe, M. W. Berns, Y. Tadir, and
B. J. Tromberg,
“Quantitative near-infrared spectroscopy of cervical dysplasia in vivo,”
Hum. Reprod., 14 2908
–2916
(1999). https://doi.org/10.1093/humrep/14.11.2908 Google Scholar
S. L. Jacques and
D. J. Mcauliffe,
“The melanosome – threshold temperature for explosive vaporization and internal absorption-coefficient during pulsed laser irradiation,”
Photochem. Photobiol., 53 769
–775
(1991). Google Scholar
C. R. Simpson, M. Kohl, M. Essenpreis, and
M. Cope,
“Near-infrared optical properties of ex vivo human skin and subcutaneous tissues measured using the Monte Carlo inversion technique,”
Phys. Med. Biol., 43 2465
–2478
(1998). https://doi.org/10.1088/0031-9155/43/9/003 Google Scholar
D. Yudovsky and
L. Pilon,
“Rapid and accurate estimation of blood saturation, melanin content, and epidermis thickness from spectral diffuse reflectance,”
Appl. Opt., 49 1707
–1719
(2010). https://doi.org/10.1364/AO.49.001707 Google Scholar
R. R. Anderson and
J. A. Parrish,
“The optics of human skin,”
J. Invest. Dermatol., 77 13
–19
(1981). https://doi.org/10.1111/1523-1747.ep12479191 Google Scholar
M. J. van Gemert, S. L. Jacques, H. J. Sterenborg, and
W. M. Star,
“Skin optics,”
IEEE Trans. Biomed. Eng., 36 1146
–1154
(1989). https://doi.org/10.1109/10.42108 Google Scholar
C. M. Balch, S. J. Soong, J. E. Gershenwald, J. F. Thompson, D. S. Reintgen, N. Cascinelli, M. Urist, K. M. McMasters, M. I. Ross, J. M. Kirkwood, M. B. Atkins, J. A. Thompson, D. G. Coit, D. Byrd, R. Desmond, Y. Zhang, P. Y. Liu, G. H. Lyman, and
A. Morabito,
“Prognostic factors analysis of 17,600 melanoma patients: validation of the American Joint Committee on Cancer melanoma staging system,”
J. Clin. Oncol., 19 3622
–3634
(2001). Google Scholar
|
||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||
