2.1.

Combined Monte Carlo/Finite-Difference Time-Domain Algorithm

In MC modeling, each tissue layer is described by a thickness $D$ and several optical properties, including the refractive index $n$ , absorption coefficient ${\mu }_a$ , scattering coefficient ${\mu }_s$ , and scattering phase function $p(\widehat{\mathbf{s}},{\widehat{\mathbf{s}}}^{\prime })$ that characterizes probability of scattering from unit direction ${\widehat{\mathbf{s}}}^{\prime }$ to unit direction $\widehat{\mathbf{s}}$ . If the tissue layer is isotropic, scattering depends only on the deflection angle $\theta$ between unit vector directions $\widehat{\mathbf{s}}$ and ${\widehat{\mathbf{s}}}^{\prime }$ .^{15, 16} The Henyey-Greenstein (HG) function has evolved as a simple and convenient expression that is frequently used to approximate the scattering patterns of tissue scatterers. This function describes the probability density for $\cos \theta$ and is specified fully by only a single model parameter as¹⁵

FDTD modeling can be used to compute scattering patterns of tissue constituents. These scattering patterns describe the intensity of scattered light at discrete angles and cannot be described analytically. Hence, a discrete random-variate generation algorithm is needed to allow sampling of the scattering direction from scattering patterns obtained using the FDTD method. Assume that the scattering pattern is denoted by $I\left(\theta \right)$ , where $\theta ∊\{0,1,\dots ,180\}$ is the scattering angle in degrees. The scattering phase function is then given by the probability mass function $p_{\mathrm{FDTD}}\left(\theta \right)$ that describes the probability that the scattering angle $\Theta =\theta$ :

Note that in defining the probability mass function in Eq. 5, the scattering pattern $I\left(\theta \right)$ has been multiplied by $\sin \theta$ .²⁸ The angular interval $\Delta \theta$ equals 1 and cancels out in Eq. 5, but it has been included for completeness. The inverse-transform technique can be used to sample from the FDTD phase function $p_{\mathrm{FDTD}}\left(\theta \right)$ .^{17, 19, 20, 38, 39, 40} First, the cumulative distribution function $F\left(\theta \right)$ is computed using

Note that Eq. 7 represents an interpolation scheme that allows generation of continuous scattering angles from an FDTD phase function. In general, the inverse-transform technique may not be highly efficient, since it requires a search on $F\left(\theta \right)$ . However, tissue scatterers are mostly forward scattering and scattering angles are usually restricted to small angles. The search on $F\left(\theta \right)$ is most of the time expected to terminate without having to scan a large number of entries. Therefore, the presented sampling method does not suffer from significant computational burden.

2.2.

Model Input

2.2.2.

Combined Monte Carlo/finite-difference time-domain simulation parameters

Epithelial tissue was modeled as consisting of two layers, the top thin layer being the epithelium with a thickness of $300\mu \mathrm{m}$ ,⁴⁷ and the second layer being the semi-infinite stroma. Table 2 lists the absorption and scattering coefficients of both the epithelium and the stroma for normal and precancerous epithelial tissue. The values listed for normal tissue were based on measurements from cervical tissue or other epithelial tissues having a similar architecture, and were reported in a recent paper by Arifler ⁵ Optical parameters used in modeling precancer were obtained by increasing the epithelial scattering coefficient by a factor of 3, increasing the stromal absorption coefficient by a factor of 2, and decreasing the stromal scattering coefficient by 25%. These changes were adopted to represent increased epithelial scattering due to structural and morphological alterations in cell nuclei, angiogenic activity, and hence increased hemoglobin content in the stroma, and degradation of collagen fibers in the stroma. The coefficients listed in Table 2 for precancerous tissue have been observed to generate reflectance spectra consistent with measurements from highly dysplastic cervical tissue.⁵ The three wavelengths chosen for simulations sample the visible part of the electromagnetic spectrum, where $\lambda =420\mathrm{nm}$ corresponds to the hemoglobin absorption peak. The refractive index of both tissue layers was set to $n=1.4$ . Note that even though the refractive index of the outside medium was assumed to be 1.35 for FDTD computations, rounding this number to 1.4 in MC modeling does not cause a significant deviation in simulation results as previously demonstrated.⁵ HG phase function with $g=0.88$ was used to describe angular scattering probability distribution in the stromal layer. For the epithelial layer, phase functions computed using the FDTD method were employed, while HG phase functions with identical $g$ values as the FDTD functions were also tested to present a comparative analysis.

Table 2

Absorption and scattering coefficients for normal and precancerous epithelial tissue.5

Layer	Normal		Precancer
	μa (cm−1)	μs (cm−1)	μa (cm−1)	μs (cm−1)
$\lambda =420\mathrm{nm}$
Epithelium	3.0	42.4	3.0	127.2
Stroma	9.09	266.3	18.18	199.7
$\lambda =500\mathrm{nm}$
Epithelium	2.0	35.6	2.0	106.8
Stroma	2.41	223.7	4.82	167.8
$\lambda =650\mathrm{nm}$
Epithelium	1.2	27.4	1.2	82.2
Stroma	0.97	172.0	1.94	129.0

Optical probing depth, or the extent of the tissue volume to which the optical probe is most sensitive, depends on the source-detector geometry. We model three fiber optic probe geometries to determine the extent of influence of using FDTD phase functions for different designs. The first design represents the simplest arrangement, where the source and detector fibers are oriented side by side and perpendicular to the surface, and is widely used in research and clinical settings.³³ Extensive efforts have recently been made to design probe geometries with increased selectivity for optical signals from the top epithelial layer, since early identification of changes in epithelial cells is expected to be of significant importance in cancer diagnostics. The second and third designs are based on previous studies that suggest using tilted source and detector fibers to limit the penetration depth of photons, ^{9, 10, 11, 12, 14, 34, 35, 36} or coupling fibers to a half-ball lens such that the flat surface of the lens faces the tissue.¹⁰ The three different probe geometries modeled are illustrated in Fig. 1 . In all cases, both the source fiber and the detector fiber were assigned a $100\mu \mathrm{m}$ diameter and a numerical aperture of 0.11 (in air), and their refractive indices were set to 1.5. The material between the fiber tips was assumed to be highly absorptive, suppressing any internal reflection of incoming photons. For the first two designs [Figs. 1a and 1b], center-to-center source-detector separation $s$ was $150\mu \mathrm{m}$ , whereas for the third design [Fig. 1c], $s$ was set to $900\mu \mathrm{m}$ . The tilt angle $\beta$ for the second design was chosen to be $30\mathrm{deg}$ . The third design involved a sapphire half-ball lens of radius $R=500\mu \mathrm{m}$ and a flat sapphire window of thickness $w=300\mu \mathrm{m}$ , while the thickness of the air gap between the fiber tips and the top of the lens was $d=200\mu \mathrm{m}$ . The refractive index of sapphire is known to be wavelength dependent and was set to 1.78 at $\lambda =420\mathrm{nm}$ , 1.77 at $\lambda =500\mathrm{nm}$ , and 1.76 at $\lambda =650\mathrm{nm}$ .

Fig. 1

Fiber optic probe geometries with (a) perpendicularly oriented, (b) tilted, and (c) half-ball lens coupled source and detector fibers. In (a) and (b), fiber tips are polished parallel to the tissue surface, and in (c), the source and detector fibers are symmetrically oriented about the central axis of the lens.

The fixed-weight Monte Carlo code used for the simulations presented in this work was implemented in C/C++ and has been previously described.¹⁰ The output of each simulation includes the total number of reflected photons detected, maximum penetration depths of detected photons, and scattering angles that detected photons undergo in each tissue layer. Three simulations, each with ${10}^8$ input photons, were carried out for each set of optical properties and phase function, and the results presented represent averages over these three simulations.

3.1.

Finite-Difference Time-Domain Simulation Results

Figure 2 shows the phase functions obtained using the FDTD method at $\lambda =420$ , 500, and $650\mathrm{nm}$ . These functions, denoted by $p\left(\theta \right)$ , describe the probability of scattering at an angle $\theta ∊\{0,1,\dots ,180\}$ and have been calculated using Eq. 5, where the descriptive subsript FDTD in $p_{\mathrm{FDTD}}\left(\theta \right)$ has been dropped for notational simplicity. The plots on the left side of the figure [Figs. 2a, 2b, 2c] present the phase functions for normal epithelial cells, and the plots on the right side [Figs. 2d, 2e, 2f] present the phase functions for precancerous epithelial cells. Note again that the results shown for each wavelength and for each diagnostic category correspond to the averages over three simulations with different cell parameters. The angular resolution for the scattering angle is $1\mathrm{deg}$ for each plot. The probability of scattering at 0 and $180\mathrm{deg}$ is zero due to the inclusion of the $\sin \theta$ factor, and these data points have been excluded from the semilog plots.

Fig. 2

Comparison of FDTD and HG phase functions for normal epithelial cells at (a) $\lambda =420\mathrm{nm}$ , $g=0.936$ ; (b) $\lambda =500\mathrm{nm}$ , $g=0.941$ ; and (c) $\lambda =650\mathrm{nm}$ , $g=0.945$ ; and for precancerous epithelial cells at (d) $\lambda =420\mathrm{nm}$ , $g=0.917$ ; (e) $\lambda =500\mathrm{nm}$ , $g=0.934$ ; and (f) $\lambda =650\mathrm{nm}$ , $g=0.934$ .

To facilitate comparison of the FDTD phase functions to analytical HG phase functions with identical anisotropy factors, Fig. 2 also shows the HG functions computed using a discretized version of Eq. 3 and plotted over the FDTD results. The subscript HG in $p_{\mathrm{HG}}\left(\theta \right)$ has again been dropped for notational simplicity. It is observed that HG phase functions are mainly good approximations of the FDTD results that show quite frequent variations over the entire angular range $\left(0\text{to}180\mathrm{deg}\right)$ . However, there exist subtle but distinct differences between FDTD and HG functions. FDTD functions are characterized by relatively higher forward scattering probability and higher high-angle $(>\sim 160\mathrm{deg})$ scattering probability for all six cases considered. It is also apparent that FDTD functions tend to suggest decreased probability of scattering for the angular range $\sim 2\text{to}20\mathrm{deg}$ , increased probability of scattering for $\sim 20\text{to}60\mathrm{deg}$ , and decreased probability of scattering for $\sim 60\text{to}110\mathrm{deg}$ compared to HG functions. For normal epithelial cells, the behavior of phase functions over the angular range $\sim 110\text{to}160\mathrm{deg}$ is wavelength dependent: for $\lambda =420\mathrm{nm}$ , scattering probability for this range is higher for the FDTD phase function compared to the HG phase function [Fig. 2a]; for $\lambda =500\mathrm{nm}$ , FDTD and HG phase functions overlap [Fig. 2b]; and for $\lambda =650\mathrm{nm}$ , scattering probability is higher for the HG phase function compared to the FDTD phase function [Fig. 2c]. Phase functions for precancerous epithelial cells demonstrate a more consistent trend over this same angular range: for all the wavelengths considered, scattering probability is higher for HG phase functions, and the largest differences between FDTD and HG functions arise for $\lambda =650\mathrm{nm}$ .

Figure 3 combines the FDTD or HG results for normal and precancerous cells on the same plot for each of the wavelengths simulated. As in the case of Fig. 2, the data points for 0 and $180\mathrm{deg}$ have been excluded from the semilog plots. FDTD results on the left side of the figure [Figs. 3a, 3b, 3c] highlight the changes in phase functions of epithelial cells due to precancer. At $\lambda =420\mathrm{nm}$ , the FDTD phase function for precancerous cells is well above the phase function for normal cells over the angular range $\sim 5\text{to}100\mathrm{deg}$ , indicating increased probability of scattering for this range. The trend reverses for $\sim 100\text{to}160\mathrm{deg}$ , where the probability of scattering for precancerous cells appears to be lower compared to the normal case. For angles $>\sim 160\mathrm{deg}$ , the scattering probability for precancerous cells exceeds the scattering probability for normal cells. At $\lambda =500$ and $650\mathrm{nm}$ , the phase functions for precancerous cells tend to be slightly higher compared to the phase functions for normal cells up to $100\mathrm{deg}$ . For the angular range $100\text{to}160\mathrm{deg}$ , the phase functions for normal and precancerous cells overlap, and for $>\sim 160\mathrm{deg}$ , the scattering probability for precancerous cells again exceeds the scattering probability for normal cells at both wavelengths. The plots on the right side of Fig. 3 [Figs. 3d, 3e, 3f] show the corresponding HG phase functions with identical anisotropy factors as the FDTD phase functions. It is apparent that differences between HG phase functions for normal and precancerous cases are very similar for all wavelengths, and scattering probability is almost always higher for precancerous cells.

Fig. 3

FDTD phase functions for normal and precancerous epithelial cells at (a) $\lambda =420\mathrm{nm}$ , (b) $\lambda =500\mathrm{nm}$ , and (c) $\lambda =650\mathrm{nm}$ . The corresponding HG phase functions with identical anisotropy factors as the FDTD phase functions are shown on the right for (d) $\lambda =420\mathrm{nm}$ , (e) $\lambda =500\mathrm{nm}$ , and (f) $\lambda =650\mathrm{nm}$ .

3.2.

Verification of the Combined Monte Carlo/Finite-Difference Time-Domain Algorithm

It is possible to keep track of the scattering angles sampled from a given probability distribution $p\left(\theta \right)$ during a Monte Carlo run. A normalized histogram generated from the sampled values simply describes the relative frequency of scattering angles and can be directly compared to the input function $p\left(\theta \right)$ . Figure 4 shows the results for a standard HG phase function with $g=0.95$ [Fig. 4a] and for an FDTD phase function with $g=0.945$ [Fig. 4b]. In both cases, the original HG or FDTD phase function used for model input is represented by a solid line. Histograms of sampled scattering angles have been generated with a bin width of $1\mathrm{deg}$ and are plotted using a circular symbol at each data point. The values for 0 and $180\mathrm{deg}$ have again been excluded from the plots.

Fig. 4

Verification of the sampling algorithm used to generate scattering angles from a given probability distribution $p\left(\theta \right)$ . (a) Sampling from a standard HG phase function with $g=0.95$ , and (b) sampling from an FDTD phase function with $g=0.945$ . In both cases, the histogram of sampled scattering angles generated during a Monte Carlo run (●) is compared to the original HG or FDTD function used as model input (—).

The results in Fig. 4 verify that the sampling algorithm can regenerate $p\left(\theta \right)$ used as input for Monte Carlo simulations. There is extremely good agreement between the original HG or FDTD phase functions and the corresponding histograms generated from the sampled scattering angles. Some deviations occur for large scattering angles close to $180\mathrm{deg}$ , but considering that the probability functions cover a dynamic range of about five orders of magnitude, it is not surprising to see slight discrepancies for large angles where the probability of scattering is usually the lowest.

Note that when scattering angles are obtained using the HG phase function in a standard MC implementation, the expression given in Eq. 1 can easily be inverted to obtain an analytical function to be used in conjunction with a random number generator.¹⁵ Here, a discretized version of the function has been created to be used as model input only for verification purposes. An angular distribution identical to that shown in Fig. 4a is obtained when the conventional inverse-transform technique¹⁵ is applied to an HG phase function with $g=0.95$ .

3.3.

Combined Monte Carlo/Finite-Difference Time-Domain Simulation Results

Monte Carlo simulations were performed to assess the influence of using FDTD phase functions on modeled tissue reflectance. To this extent, simulations with FDTD phase functions were repeated using standard HG phase functions with identical anisotropy factors and keeping all the other optical properties constant. Parameters analyzed include detected reflectance intensity that corresponds to the number of photons remitted and collected out of ${10}^8$ input photons, average maximum penetration depth of detected photons, which provides information about the optical probing depth of a given source-detector fiber geometry, and percentage of photons collected from epithelium that reveals the layer selectivity associated with the source-detector geometry in question. As indicated before, three simulations were carried out for each case considered, and the results presented correspond to the averages over three different simulations. The standard errors computed were negligibly small compared to the averages $(<1\%)$ , indicating convergence of Monte Carlo results.

3.3.1.

Detected reflectance intensity

Figure 5 shows the detected reflectance intensity at $\lambda =420$ , 500, and $650\mathrm{nm}$ for three different fiber optic probe designs illustrated in Fig. 1. Simulation results for HG and FDTD phase functions are plotted separately for both normal and precancerous epithelial tissue. Error bars corresponding to the standard error values over three simulations are the same size as or smaller than the symbols shown.

Fig. 5

Number of detected photons for normal and precancerous epithelial tissue at $\lambda =420$ , 500, and $650\mathrm{nm}$ for three different probe designs, where the source and detector fibers are (a) oriented perpendicular to the tissue surface, (b) tilted at $30\mathrm{deg}$ with respect to the tissue surface, and (c) half-ball lens coupled.

The results in Fig. 5 demonstrate that for the probe design with perpendicularly oriented source and detector fibers, the detected reflectance intensity does not depend on the phase function used, since the total number of photons detected is the same for both HG and FDTD phase functions [Fig. 5a]. For probe designs with tilted and lens coupled fibers, however, the exact form of the phase function used strongly affects the detected reflectance intensity, especially for precancerous tissue [Figs. 5b and 5c]. In almost all cases, the number of photons detected decreases when FDTD phase functions are used. The only exception to this trend occurs for normal tissue at $\lambda =420\mathrm{nm}$ when tilted fibers are used, and the number of photons detected in this case increases when the FDTD phase function is used. The probe design where the source and detector fibers are tilted at $30\mathrm{deg}$ with respect to the tissue surface appears to be the most sensitive to the form of the phase function; HG and FDTD simulation results for precancerous tissue differ by $\sim 20\text{to}30\%$ , where the high end of this percentage range corresponds to the values at $\lambda =650\mathrm{nm}$ . Note that precancer leads to a decrease in the number of detected photons for perpendicularly oriented fibers, whereas it causes a significant increase in reflectance intensity for the other two probe designs.

3.3.2.

Optical probing depth profiles

Maximum penetration depths of all the photons collected have been recorded during the MC simulations. Note that maximum penetration depth for a given photon is defined to be the maximum depth in tissue at which the photon undergoes a scattering event. Average maximum penetration depth corresponds to the maximum penetration depth averaged over all detected photons, and this parameter gives valuable information regarding the optical probing depth of a given probe design. Figure 6 shows the average maximum penetration depths computed using HG and FDTD phase functions at $\lambda =420$ , 500, and $650\mathrm{nm}$ . Standard error bars for three simulations are again the same size as or smaller than the symbols shown. Horizontal lines in Fig. 6 represent the boundary between the epithelium and the stroma.

Fig. 6

Average maximum penetration depths for normal and precancerous epithelial tissue at $\lambda =420$ , 500, and $650\mathrm{nm}$ for three different probe designs, where the source and detector fibers are (a) oriented perpendicular to the tissue surface, (b) tilted at $30\mathrm{deg}$ with respect to the tissue surface, and (c) half-ball lens coupled. Horizontal lines represent the boundary between the epithelium and the stroma.

The results in Fig. 6 show that for all three probe designs simulated, the average maximum penetration depths are not strongly sensitive to the exact form of the phase function used. This is in contrast to the sensitivity observed in Fig. 5 for detected reflectance intensity. In most cases, the use of FDTD phase functions leads to a slight increase in average maximum penetration depth, and the most significant difference $(\sim 15\%)$ occurs at $\lambda =650\mathrm{nm}$ for normal tissue when tilted fibers are used. It is important to note that for the probe design with perpendicularly oriented fibers, the average maximum penetration depth is always greater than the thickness of the epithelial layer [Fig. 6a], suggesting that photons penetrate deep into the stroma. Probe geometries with tilted and lens coupled source and detector fibers, on the other hand, are most of the time preferentially sensitive to the epithelium, as evidenced by average maximum penetration depths that are smaller than the thickness of the epithelial layer [Figs. 6b and 6c]. Also note that the penetration depths calculated for precancerous tissue are smaller compared to the depths calculated for normal tissue, and this suggests that the same geometry probes more superficial tissue depths when precancer is present.

3.3.3.

Layer selectivity

Figure 7 shows the percentage of detected photons that only scatter in the epithelium without penetrating into the stroma. This parameter is expected to provide information about the layer selectivity of a given probe design. Simulation results for HG and FDTD phase functions are plotted separately at $\lambda =420$ , 500, and $650\mathrm{nm}$ . Error bars corresponding to the standard error values over three simulations are negligibly small and have been omitted from the plots.

Fig. 7

Percentage of photons collected from epithelium for normal and precancerous epithelial tissue at $\lambda =420$ , 500, and $650\mathrm{nm}$ for three different probe designs, where the source and detector fibers are (a) oriented perpendicular to the tissue surface, (b) tilted at $30\mathrm{deg}$ with respect to the tissue surface, and (c) half-ball lens coupled.

As in the case of average maximum penetration depths, the results for percentage of photons detected from epithelium are not very sensitive to the exact form of the phase function used. The most significant differences are observed for tilted fibers [Fig. 7b], where FDTD results are characterized by up to a 5% decrease in fraction of photons that only scatter in the epithelium and do not penetrate into the stroma before being remitted and collected. Figure 7 also shows that sensitivity to the epithelial layer increases with precancer, and layer selectivity is best achieved with half-ball lens coupled fibers, where the fraction of photons collected from epithelium is always $>\sim 90\%$ .

4.1.

Finite-Difference Time-Domain Phase Functions

Since the HG phase function can be described by a simple analytical expression, it is commonly used to describe angular scattering properties of tissues when developing photon propagation models or when implementing inverse algorithms to extract tissue optical properties from measurements. FDTD phase functions and the corresponding HG functions presented in Fig. 2 demonstrate that HG approximation is indeed representative of the overall characteristics of more realistic FDTD functions that incorporate complex interior morphology and dielectric structure of epithelial cells. There are, however, subtle differences between HG and FDTD phase functions, and the main goal of this study was to determine whether these subtle differences lead to significant changes in spectroscopic reflectance signals. In all cases, forward scattering is higher for FDTD functions, and inclusion of nucleus in FDTD modeling is likely to account for this behavior.²⁴ High-angle scattering is sensitive to smaller inclusions such as organelles or nuclear refractive index heterogeneities in the form of chromatin clumps.^{24, 26, 27} Therefore, increased high-angle scattering probability observed for FDTD functions is again a result of using more realistic cell geometries in FDTD computations. FDTD phase functions are also characterized by a sharp increase in scattering probability at angles very close to $180\mathrm{deg}$ . This behavior is consistent with earlier measurements by Flock, Wilson, and Patterson,⁴⁸ and Mourant, ⁴⁹ and also with recent findings by Mishchenko ⁵⁰ about scattering by random particulate media. It is apparent that the differences between HG and FDTD phase functions at intermediate angles depend on wavelength as well as on diagnostic category.

Phase functions computed using the FDTD method and presented in Figs. 2 and 3 reflect morphological and structural changes associated with precancer. Increased probability of scattering up to $\sim 100\mathrm{deg}$ for precancerous cells is most likely due to increased nuclear size and nuclear refractive index.^{24, 26, 27} The most significant increase tends to occur at $420\mathrm{nm}$ , the smallest wavelength considered. This is consistent with the fact that for small wavelengths, changes in size and refractive index of the scattering structure are expected to have a more dramatic effect on its small-angle scattering properties. Since the phase functions presented are by definition normalized, a large increase in small-angle scattering probability at $420\mathrm{nm}$ for precancerous cells has resulted in a relative decrease of scattering probability over $\sim 100\text{to}160\mathrm{deg}$ . Increased scattering probability of precancerous cells at all wavelengths for $>\sim 160\mathrm{deg}$ is most likely due to aggregation of chromatin clumps and hence increased variation in nuclear refractive index profile.^{24, 26, 27} Increase in high-angle scattering from tumorigenic cells has been previously reported by Ramachandran, ⁵¹ who carried out angle-dependent light-scattering measurements on tumorigenic and nontumorigenic cells.

4.2.

Influence of Using Finite-Difference Time-Domain Phase Functions on Modeled Reflectance

The results presented indicate that when the optical illumination and detection system is sensitive to the epithelial layer, detected reflectance intensity is highly dependent on the exact form of the phase function used to describe scattering from epithelial cells. MC simulations carried out using HG phase functions with identical anisotropy factors as the FDTD phase functions, and keeping all other optical properties constant, illustrate that there can be extensive differences in the number of photons collected from tissue, and the extent of these differences can vary depending on the probing depth of the probe geometry, as well as on changes in optical properties due to precancer.

When perpendicularly oriented source and detector fibers are used to illuminate tissue and to collect remitted photons, the form of the phase function employed in simulations does not have any influence on the number of detected photons, as evidenced in Fig. 5. The probe geometry with perpendicularly oriented source and detector fibers is largely sensitive to deep portions of the tissue; most of the photons collected enter the stroma and scatter multiple times before being remitted at the tissue surface. This is verified later in Figs. 6 and 7 that provide information about the probing depth of the probe geometry and the percentage of photons that only scatter in epithelium and do not enter the stroma before being detected. The results shown indicate that this geometry is characterized by a probing depth that is much greater than the epithelial thickness and cannot selectively target the epithelial layer. In this case, it is not surprising at all that the form of the epithelial phase function does not have any influence on detected reflectance. Note that increased epithelial scattering, increased stromal absorption, and decreased stromal scattering that accompany precancer progression lead to an overall decrease in reflectance intensity, again consistent with the fact that perpendicularly oriented source and detector fibers are sensitive to stromal optical properties. The most significant difference in probing depth and percentage of photons collected from epithelium for normal and precancerous tissue occurs at $420\mathrm{nm}$ , where stromal absorption is the greatest.

The probe geometry with tilted source and detector fibers shows increased selectivity for epithelium and hence increased sensitivity to the exact form of the phase function used. Moreover, when the phase functions presented in Fig. 2 are examined in detail, it is possible to find a one-to-one correspondence between differences in HG and FDTD phase functions and the resulting differences in the number of detected photons shown in Fig. 5. For example, it is obvious that tilted source and detector fibers are very sensitive to the phase function over the angular range $\sim 100\text{to}150\mathrm{deg}$ . At $420\mathrm{nm}$ , the FDTD phase function exceeds the HG phase function for normal epithelial cells, resulting in a greater number of detected photons when the FDTD phase function is used in Monte Carlo simulations. At $500\mathrm{nm}$ , HG and FDTD phase functions for normal epithelial cells overlap, and no significant phase function dependency can be observed in detected reflectance intensity. Finally, at $650\mathrm{nm}$ , the FDTD phase function for normal epithelial cells is below the corresponding HG phase function, and this directly leads to a decrease in detected reflectance intensity when the FDTD function is used in MC simulations. The observation that tilted fibers are sensitive to the phase function over $\sim 100\text{to}150\mathrm{deg}$ holds true for precancerous tissue as well. The most significant difference in detected reflectance intensity occurs at $650\mathrm{nm}$ , where the extent of differences between HG and FDTD phase functions over this angular range is the greatest. In contrast to the case with perpendicularly oriented fibers, precancer leads to an increase in detected reflectance intensity, since increased epithelial scattering is the decisive factor and precancerous changes in stromal optical properties do not contribute extensively to alterations in reflectance profiles. It is important to note that the form of the phase function used has a greater influence on detected reflectance when precancerous tissue is modeled. This is in line with increased epithelial scattering associated with precancer, and hence a greater sensitivity to any changes in angular scattering properties of epithelial cells.

The probe design with lens coupled fibers shows the greatest selectivity for the epithelial layer, as evidenced in Figs. 6 and 7. The use of FDTD phase functions leads to a consistent decrease in detected reflectance intensity. This behavior suggests that lens coupled fibers are likely to be most sensitive to the phase function over $\sim 60\text{to}110\mathrm{deg}$ . Precancer leads to an increase in the number of detected photons and an increase in sensitivity to the phase function, as was the case with tilted fibers.

In addition to basing the relationship between the phase functions and the corresponding Monte Carlo modeling results on visual inspection of these functions, we can carry out a more quantitative analysis to identify specific angular ranges to which each probe design is most sensitive. This can be achieved computationally by keeping track of all the angles each detected photon scatters through before getting detected. Figure 8 shows typical distributions of angles through which detected photons scatter in the epithelium for probe designs where the source and detector fibers are tilted and lens coupled. These results have been obtained using HG or FDTD phase function and optical properties for precancerous epithelial tissue at $500\mathrm{nm}$ . The normalized histograms generated with a bin width of $1\mathrm{deg}$ validate our observations that tilted fibers are sensitive to the phase function over $\sim 100\text{to}150\mathrm{deg}$ , whereas lens coupled fibers are sensitive to the phase function over $\sim 60\text{to}110\mathrm{deg}$ . Similar distributions (not shown) are obtained for other combinations of optical properties and HG or FDTD phase functions. It is important to note that when each fiber is tilted at $30\mathrm{deg}$ , most photons have to undergo scattering through $120\mathrm{deg}$ to reach the detector fiber. Therefore, the peak in angular distribution at around $120\mathrm{deg}$ is not surprising for the case of tilted fibers [Figs. 8a and 8c]. The distribution of scattering angles peaks at around $90\mathrm{deg}$ for lens coupled fibers [Figs. 8b and 8d], indicating that this geometry can be considered to be roughly equivalent to a design where the source and detector fibers are tilted at $45\mathrm{deg}$ . Also note that the distributions shown in Fig. 8 represent averages over three simulations with ${10}^7$ input photons. Keeping track of all scattering angles during MC simulations is memory intensive, and a smaller number of input photons were used to generate the histograms presented. Even though the distributions are characterized by sudden jumps, they provide a useful insight into the basic trends observed. In our simulations, we also recorded the number of scattering events each detected photon undergoes in each tissue layer. The results for both tilted and lens coupled fiber geometries suggest that the majority of detected photons either undergo a single intermediate-angle scattering event in the epithelium, or experience one to three near-forward or small-angle scattering events, followed by an intermediate-angle scattering event that reverses their path toward the detector before they even reach the stroma. This is consistent with the distributions in Fig. 8, where small-angle scattering probability is observed to be highly pronounced, along with the apparent enhancement of scattering probability over an intermediate angular interval identified earlier as the sensitivity range for each probe design.

Fig. 8

Distribution of angles through which detected photons scatter in the epithelium for two probe designs, where the source and detector fibers are (a) and (c) tilted at $30\mathrm{deg}$ with respect to the tissue surface, and (b) and (d) half-ball lens coupled. The normalized histograms shown correspond to simulation results obtained using HG or FDTD phase function and optical properties for precancerous epithelial tissue at $\lambda =500\mathrm{nm}$ .

It should be emphasized that probe designs used in this study were characterized by small-sized fibers with a low numerical aperture. This is especially important for tilted and lens coupled fiber geometries that are optimized to demonstrate significant sensitivity to the top epithelial layer and hence to the epithelial phase function. If we increase the size or the numerical aperture of the source and detector fibers, or if we change the source-detector separation, photons will no longer be effectively localized to the epithelium and the contribution of stromal scattering will increase, diminishing the influence of epithelial phase function on detected reflectance. Using small-sized fibers with a low numerical aperture surely results in collection of fewer photons, but a decrease in optical signal strength down to an acceptable level is a tradeoff to ensure selective probing of the epithelial layer, which is of extreme importance in cancer diagnostics.¹⁰

4.3.

Implications

The work described here combines MC modeling with the powerful FDTD technique that can be used to generate realistic tissue phase functions. The results presented provide strong evidence that the form of the phase function can play an important role in determining the reflectance profile of tissues. MC simulations indicate that when the probe geometry is sensitive to the epithelium, subtle differences in the shape of epithelial phase functions translate into significant changes in the number of photons detected; FDTD and HG phase functions with the same anisotropy factor but different shapes can lead to significantly different results for detected reflectance intensity. It is interesting to note, however, that these differences do not lead to major changes in the probing depth of the geometry considered. This suggests that even though the form of the phase function does not affect the depth range to which the probe geometry is sensitive, it may lead to the collection of more or fewer photons from the same tissue volume, resulting in changes in the intensity of the optical signals acquired. The form of the phase function may also alter the wavelength dependence of the measured light intensity.

The dependence of tissue reflectance on the exact form of the phase function needs to be taken into account when developing photon propagation models to computationally predict the optical response of tissues. It is also necessary to consider the influence of the phase function when implementing inverse algorithms to analyze optical measurements and to estimate tissue optical parameters from these measurements. Most studies tend to focus on extracting scattering and absorption coefficients from optical signals by assuming a simple analytical form for the phase function. This may lead to possible errors in reconstructed optical parameters and hence to misinterpretation of the measurements, especially when small source-detector separations are used.

It should be noted that this study focused on analyzing the effect of angular scattering properties of epithelial cells. The FDTD method has recently been applied to study micro-optical properties of collagen fibers that constitute the stromal layer.⁵² The work presented here can be extended to use combined MC/FDTD modeling to assess the influence of using realistic phase functions for the stroma as well. In that case, probe designs that target the stromal layer are expected to demonstrate significant sensitivity to the exact form of the phase function used.