Ultraviolet-visible diffuse reflectance spectroscopy has been explored for the early detection of epithelial cancers for decades.¹ In this technique, an accurate model of light transport is essential to quantitatively extract optical properties from measured diffuse reflectance spectra. Epithelial tissues with tumors are frequently considered as a multilayered tissue model with tumor-like heterogeneities.² In this situation, the Monte Carlo (MC) method provides a flexible tool to model light transport. Since the MC method can solve a radiative transport equation with any accuracy³ for a complex tissue model and probe geometry, it is considered the gold standard for modeling light transport in turbid media. However, the main drawback of the MC method is the requirement of intensive computation to achieve results with desirable accuracy, which makes it extremely time consuming. Several methods have been proposed to speed up the MC method for modeling light transport in complex tissue models. Liu et al.⁴ presented a scaling method for fast MC simulation of diffuse reflectance spectra from multilayered turbid media. Hayakawa et al.⁵ proposed a perturbation Monte Carlo (pMC) method to solve inverse photon migration problems in a two-layered tissue model based on spatially resolved diffuse reflectance and validated this method experimentally.⁶ To our best knowledge, there has been no effort in the literature to speed up the MC method in a multilayered tissue model with finite-size tumor-like heterogeneities. Theoretically, the pMC method may be used in this case, but the applicable range of optical properties in the tissue model and the heterogeneity will be limited. Moreover, the probe configuration is always fixed in the previous pMC methods, which would cause significant inconvenience for a fiber-optics probe involving multiple configurations such as multiple source-detector separations.

In this letter, we present a hybrid method that combines the multilayered scaling method⁴ and the pMC method⁵ mentioned above for fast MC simulation of diffuse reflectance from a multilayered tissue model with tumor-like heterogeneities. Our method consists of two steps as shown in Fig. 1. The first step applies the multilayered scaling method on a set of photon trajectory information including the exit weight and the $x$ and $y$ offsets in each random walk step of all survival photons escaping from the top surface of the tissue model, generated from a single baseline simulation to scale the exit weight and exit distance of photons for the multilayered tissue model without heterogeneities. In the second step, a convolution scheme is used first to determine the probability of a survival photon collected by the fiber-optic probe geometry of interest.⁷ Then the second set of photon trajectory information including the locations of all collision events for each collected photon, which is generated from the same baseline MC simulation, is processed to determine the path length and the number of collisions of photon spent in the tumor. Finally, the scaling result, that is, the exit weight of collected photons, as well as the path length and the number of collisions spent in the tumor, is utilized by the pMC method to compute the diffuse reflectance for the given probe configuration.

A previous MC code⁴ was modified to create a photon trajectory database for scaling and perturbation. A single simulation was run for a homogeneous baseline tissue model, in which ${\mu }_a=0{\mathrm{cm}}^{-1}$, ${\mu }_s=100{\mathrm{cm}}^{-1}$, and the anisotropy factor $g=0.8$. The refractive indices of the medium above the tissue model, the tissue model, and the medium below the tissue model were set at 1.47, 1.4, and 1.4, respectively. The thickness of the tissue model was set at 4 cm to mimic a semi-infinite medium. A total of ${10}^7$ photons were launched at the origin of a Cartesian coordinate system to obtain the impulse response of the tissue model in diffuse reflectance. When a photon exits from the top surface of the tissue model, its exit angle relative to the $z$-axis is calculated. If the exit angle is smaller than the cut-off angle defined by an numerical aperture (NA) of 0.22, the relevant trajectory information of this photon is stored in a numerical array. Approximately $2.4\times {10}^5$ photons were detected in this manner, and a total memory of 10 gigabytes (GB) was needed for the storage of the trajectory data. Then based on the stored trajectory data, the multilayered scaling method⁶ and the pMC method⁷ are sequentially carried out as described previously to estimate diffuse reflectance from the multilayered tissue model with finite-size heterogeneities. Both the scaling and pMC methods were coded and run in Matlab 10 (Mathworks, Natick, Massachusetts, US).

A basal cell carcinoma (BCC) skin tissue model was used to evaluate the effectiveness of the hybrid method. The BCC usually originates from the basal layer of the epidermis and frequently grows downward deeply into the dermis,⁸ thus it is induced into the dermis in our theoretical tissue model as shown by the cross-sectional view in Fig. 2. The thickness of the epidermis was set at 80 µm and the thickness of the dermis was set at 4 cm to mimic a thick skin tissue. The epidermal thickness is representative of that on the neck and back.⁹ The length, width, and thickness of the BCC tumor were all set at 400 µm. The optical properties of the epidermis and dermis were selected from the literature¹⁰ and listed in Table 1. A refractive index of 1.4 and an anisotropy factor of 0.8 were used in the entire tissue model including the BCC tumor. The absorption and scattering coefficients of the tumor were varied sequentially to investigate the valid range of the hybrid method. The fiber configuration in Fig. 2 shows the source and detector fibers placed side by side, with both perpendicular to the tissue surface. The bisecting line between the two fibers overlaps with the middle line of the tumor in the width dimension, that is, the $x$ dimension in Fig. 2. The two fibers both had a core diameter of 200 µm and NA value of 0.22. The refractive indices of the fibers were set at 1.47. In total, two sets of tests were performed. In the first set, the absorption coefficient, ${\mu }_a$, of the tumor was varied from 1% to 400% of that of the dermis, while the scattering coefficient, ${\mu }_s$, of the tumor was kept identical to that of the dermis. In the second set, the scattering coefficient, ${\mu }_s$, of the tumor was varied from 25% to 190% of that of the dermis, while the absorption coefficient, ${\mu }_a$, of the tumor was kept identical to that of the dermis. Diffuse reflectance values, which refer to the ratio between detected and incident powers in this paper, calculated by the hybrid method were compared to those from independent MC simulations (no scaling or perturbation methods were used), which were run by using a previously validated MC code,¹¹ to evaluate the effectiveness of the hybrid method.

Table 1

Optical properties of BCC tissue model at 500 nm.10

Each independent simulation was run five times and 10 million photons were used. The percent deviation in diffuse reflectance between results calculated by the hybrid method and those simulated independently was calculated by Eq. (1) to quantify the accuracy of the calculated results.

where “Hybrid” refers to the diffuse reflectance value calculated by the hybrid method and “Simulated” refers to the mean of simulated diffuse reflectance values from five runs of the independent simulation on the same tissue model. The percent deviations of five individually simulated diffuse reflectance values relative to their mean were also calculated in the same manner. The 95% confidence interval (CI) of the percent deviation of simulated diffuse reflectance values relative to their mean was then estimated by Eq. (2). where $m$ is the number of runs ($m=5$), and “mean” and “std” refer to the mean and standard deviation of the percent deviations for simulated diffuse reflectance values, respectively. It should be noted that the mean of the percent deviations calculated in this manner is always zero as indicated by the circles in Figs. 3(b) and 4(b).

Figure 3(a) shows the comparison in diffuse reflectance values between the hybrid method and independent MC simulations for ${\mu }_a$ varying from 1% to 400% of the dermal value. The two sets of symbols completely overlap at nearly every point, which indicates the high accuracy of the hybrid method. Figure 3(b) shows the percent deviations of diffuse reflectance calculated by the hybrid method according to Eq. (1) and the 95% CIs of the percent deviations of simulated diffuse reflectance calculated by Eq. (2) as indicated by the error bars. The percent deviations for the hybrid method are always smaller than 5% when ${\mu }_a$ is varied from 1% to 390% of the dermal value. Moreover, they are all close to or within the 95% CIs of the percent deviations of simulated diffuse reflectance values. Figure 4(a) compares diffuse reflectance values between the hybrid method and independent MC simulations for ${\mu }_s$ varying from 25% to 190% of the dermal value. Figure 4(b) compares the percent deviations of diffuse reflectance values obtained by the hybrid method with the 95% CIs of the percent deviations for independently simulated diffuse reflectance values as indicated by the error bars. The percent deviations for the hybrid method are all less than 5% when ${\mu }_s$ is varied from 30% to 180% of the dermal value. Moreover, they all fall within or close to the 95% CIs of the percent deviations of simulated reflectance values. The following interesting trends have been observed in Figs. 3 and 4. The hybrid method overestimates the reflectance value when the absorption coefficient of the tumor is larger than that of the dermis and underestimates the reflectance value when the absorption coefficient of the tumor is smaller than that of the dermis. The trend for the scattering coefficient is opposite. Similar trends have been observed for the pMC method.^5,6 Thus the trends are most likely caused by the perturbation part of the hybrid method.

All tests were performed on a laptop computer with an Intel Core i5 CPU and 4 GB memory. Ten hours were needed to run the baseline simulation and generate the photon trajectory information. It should be noted that the major portion of time in the baseline simulation was spent in saving data. It took about 1.5 min for the scaling step, 15 s for the convolution scheme to calculate the probability of each survival photon being collected, and about 6 min to determine the path length and the number of collisions spent in the tumor. The pMC step took about 90 ms for each set of optical properties to yield the final diffuse reflectance. In contrast, it took about 30 min to run one independent MC simulation to get the same value.

In summary, we have developed a hybrid method for fast MC simulation of diffuse reflectance from a multilayered tissue model with tumor-like heterogeneities. In the hybrid method, one needs to perform the multilayered scaling first to find the exit weight and exit distance for every survival photon in the multilayered tissue model without considering the tumor-like heterogeneity, then use the photon trajectory information recorded for the entire baseline medium to determine the portions of the optical path and number of collisions in the tumor region and scale them to perform perturbation for the multilayered tissue model with an arbitrary tumor-like heterogeneity. The latter step is different from all the perturbation methods that have been reported in the literature to our best knowledge, which is the key step that distinguishes the hybrid method from the simple addition of the two existing methods. The hybrid method takes advantage of the proven high accuracy of the multilayered scaling method in the simulation of a multilayered tissue model so that the perturbation method is applied to only the heterogeneities instead of the entire tissue model. This feature significantly expands the range of applicable tissue models compared to pMC methods alone. Moreover, the hybrid method only requires a single baseline simulation to generate the photon trajectory information and applies it to a multilayered tissue model embedded with tumor-like heterogeneities, in which all tissue layers and tumors could have arbitrary absorption and scattering coefficients and dimensions. This advantage will speed up the computation by several orders of magnitude. Therefore the method is suitable for simulating diffuse reflectance spectra or creating a MC database to extract optical properties of a multilayered tissue model with tumor-like heterogeneities from a diffuse reflectance measurement. Such a method can be useful in the early diagnosis of epithelial cancer using diffuse reflectance spectroscopy.