2.1.

Monte Carlo Probability of Visitation and Detection $P_{\mathrm{V}\cap \mathrm{D}}$

We wish to characterize the spatial distribution of only those photons that are launched by a specified source and subsequently captured by a detector of interest. We accomplish this by computing a “contributon” response function that was first developed in the nuclear engineering community for the solution of deep-penetration nuclear transmission problems using MC methods.^28–34 In the transmission problems, a surface between the source and detector is first specified. A forward simulation for photon transport from the source is then matched with an adjoint simulation of photon transport from the detector over the midway surface using the contributon response function. In a reflectance geometry, we use this idea to determine the probability that detected photon packets have visited a depth $d$ within the tissue by defining the surface at the selected depth as our midway surface. We determine the probability that photon packets from the source “visit” $d$, $P(\mathrm{V})$, and then determine the probability that they will subsequently be detected, $P(\mathrm{D}|\mathrm{V})$. Bayes theorem³⁵ is used to determine the probability of visitation and detection, $P_{\mathrm{V}\cap \mathrm{D}}$

One method to determine only those photon packets that originate from the source, travel to the midway surface, and subsequently arrive at the detector, would be to match the radiance determined by a forward simulation from the source with an adjoint simulation from the detector over the midway surface. This approach works well when the source and detector are “small” relative to the midway surface.^14,36 In the biomedical optics community, the use of such coupled forward-adjoint simulations has been used to address fluorescence excitation and detection,³⁷ “photon hitting density” maps,⁶ and tomographic sensitivity analysis for spatially resolved reflectance.^14,36,38

In this work, we produce $P_{\mathrm{V}\cap \mathrm{D}}$ distributions for SFD measurements using a single conventional MC simulation. We take this approach because MC simulations of SFD methods utilize a “small” source and “large” detector, i.e., light is injected into the medium at a single point while detection occurs at all locations on the tissue surface. In such a scenario, a conventional MC approach provides better computational efficiency relative to coupled forward-adjoint methods.³⁸ Note that the SFD situation is unique relative to other diffuse optical methods where both “small” sources and “small” detectors are typically used. Specifically, we create depth-dependent SFD $P_{\mathrm{V}\cap \mathrm{D}}(z)$ distributions by defining parallel $x-y$ planes, placed at regular intervals of 0.01 mm below the tissue surface, as our midway surfaces. We choose this interval size to provide submillimeter resolution for our probability distributions. We then perform an MC simulation that tracks those photon packets that both visit these depths and are subsequently detected. Figure 1 shows a schematic of the tissue segmented into layers using planes located at depths $z_0=0,z_1,z_2,\cdots$ separated by height $\mathrm{\Delta }z$. These planes define midway surfaces at depths $z_i$ at which $P_{\mathrm{V}\cap \mathrm{D}}(z_i)$ is determined. We calculate $P_{\mathrm{V}\cap \mathrm{D}}(z)$ by tracking each photon packet from the source to these various midway surfaces and from those surfaces to the detector. In doing so, these calculations inherently combine the probability of a photon packet visiting the midway surface after being launched by the source with the probability of detection after visiting the midway surface.

Fig. 1

Schematic of an SFD $P_{\mathrm{V}\cap \mathrm{D}}$ MC simulation within a tissue subdivided into layered surfaces.

The $P_{\mathrm{V}\cap \mathrm{D}}(z_i)$ tallies are computed using an MC simulation of radiative transport by determining the midway surfaces visited by each detected photon packet. As mentioned earlier, we utilize an MC simulation using the MSM²⁷ that tallies the photon packet weights directly in the SFD. A narrow, collimated beam normally incident on the tissue surface is used as the source within the simulation and the photon packets are transported using conventional MC propagation.^39,40 The complex weights specified by the MSM for the specified spatial frequency provide both the amplitude change and phase shifts that result in photon packet propagation from source to detector. When utilizing illumination modulated along the $x$-axis at spatial frequency $f_x$, the SFD complex tally is²⁷

For each detected photon packet, the SFD complex random variable [Eq. (2)] is tallied for each midway surface that the photon packet crossed and notated as $\xi (z_i)$. For example, if the photon packet reached a maximum depth $z_{\max }$ prior to detection, the photon packet crossed all midway surfaces residing at depths less than $z_{\max }$. For the photon packet trajectory shown in Fig. 1, the detected photon packet weight is tallied to the midway surface $z_3$ as well as all midway surfaces residing at shallower depths because the photon packet also crossed those surfaces. Once $N$ photon packet trajectories are simulated, the expected value of $\xi$ for depth $z_i$, $E[\xi (z_i)]$ produces the probability that the trajectory of a detected photon packet crossed depth $z_i$, $P_{\mathrm{V}\cap \mathrm{D}}(z=z_i)$

2.2.

Monte Carlo Maximum Depth of Penetration $P_{z_{\text{max}}}$

The $P_{\mathrm{V}\cap \mathrm{D}}(z)$ distribution over all depths $z$ does not result in a directly measurable quantity because photon packets that contribute to the $P_{\mathrm{V}\cap \mathrm{D}}$ tally at a surface $z_i$ also contribute to $P_{\mathrm{V}\cap \mathrm{D}}$ tallies at all locations shallower than $z_i$. However, we can use $P_{\mathrm{V}\cap \mathrm{D}}(z)$ to derive a distribution that isolates the contribution of each detected photon packet to a single bin corresponding to the maximum depth visited by that photon packet trajectory. In doing so, we ensure that each detected photon packet is tallied to only a single depth bin within the tissue. This distribution is formed by taking differences of the $P_{\mathrm{V}\cap \mathrm{D}}(z)$ tally in successive bins

We call $P_{z_{\max }}(z)$ the “$z_{\max }$” distribution which isolates those detected photon packets that crossed into depth $z_i$ but did not cross into $z_{i+1}$. Because each photon packet is tallied only once at its maximum depth of propagation “$z_{\max }$,” $P_{z_{\max }}(z)$ has the property that its integral over all depths $z$ recovers all the weight of all the detected photon packets, which is equivalent to the total diffuse reflectance, $R_{\mathrm{d}}$

If the upper limit of integration on the right-hand side of the above equation is taken instead to some finite depth $d$, the result would tally only those photon packets that contribute to the reflectance and whose trajectories were restricted to depths less than or equal to $d$. We notate this as $P_{z_{\max }}(z\le d)$ or

Division of $P_{z_{\max }}(z\le d)$ by $R_{\mathrm{d}}$ produces a probability distribution function that describes the fraction $(X)$ of the detected light that visited tissue depths $d$ or less

2.3.

Experimental Validation

We used experimental SFD measurements to validate our computational results. These measurements were taken in a two-layer phantom designed to determine the optical sampling depth as a function of spatial frequency. The two-layer phantom was housed in a container with ($L\times W\times H$) dimensions of $7.2\mathrm{cm}\times 10.8\mathrm{cm}\times 6.1\mathrm{cm}$. The top layer of the phantom was a liquid composed of water, nigrosin, and ${\mathrm{TiO}}_2$ particles with optical properties ${\mu }_{\mathrm{s}}^{\prime }/{\mu }_{\mathrm{a}}=100$ and $l^{*}=1/({\mu }_{\mathrm{a}}+{\mu }_{\mathrm{s}}^{\prime })=2\mathrm{mm}$ at $\lambda =731\mathrm{nm}$. The bottom layer was a highly absorbing solid phantom composed of agar, water, and nigrosin, and occupied a total volume of 350 mL in the container. The top layer thickness $d$ was varied from $[0-7.5]l^{*}$. Figure 2 provides a schematic of the setup.

Fig. 2

Experimental setup and schematic of two-layer phantom with top layer thickness $d$ varying from $[0-7.5]l^{*}$ where $l^{*}=2\mathrm{mm}$. Theta ($\theta$) is 15 deg. This experimental setup was used to acquire reflectance as a function of spatial frequency $(f_x)$ and layer thickness $d$. All measurements were performed at $\lambda =731\mathrm{nm}$.

We used the OxImager RS SFDI system (Modulated Imaging Inc., Irvine, California) to measure the two-layer phantom. The SFDI system utilizes crossed linear polarizers in front of the projection and detection lenses to minimize the effect of specular reflection and select for diffuse reflection. The projection field of view (FOV) was $20\mathrm{cm}\times 15\mathrm{cm}$ and directed to the tissue surface at an angle of 15 deg relative to the surface normal. Detection was performed perpendicular to the surface of the phantom with an FOV of $8\mathrm{cm}\times 6\mathrm{cm}$ (effective NA = 0.253). A vertical translation stage was used to support and adjust the height of the container with the two-layer phantom. After taking the first measurement of the highly absorbing solid phantom (i.e., $d=0\mathrm{mm}$), incremental amounts of the liquid phantom were added successively with a predetermined volume such that the top liquid layer thickness $d$ above the solid phantom increased by 0.5 mm between each measurement following the $d=0\mathrm{mm}$ measurement. The micrometer on the translation stage was used to lower the two-layer phantom system by 0.5 mm following each measurement such that the top surface of the two-layer phantom system remained at a constant image plane for all top layer thicknesses measured. This was done to avoid the need for height correction during data processing. All measurements were performed at a wavelength $\lambda =731\mathrm{nm}$ with spatial frequencies $f_x=0$, 0.0125, 0.025, 0.0375, 0.05, 0.0625, 0.075, 0.0875, 0.1, 0.125, 0.15, 0.175, 0.2, 0.25, and $0.3{\mathrm{mm}}^{-1}$. At each spatial frequency (in this case, along one spatial dimension $x$), raw reflectance images at three different phases (0, $2\pi /3$, and $4\pi /3$ radians) were sequentially projected onto the phantom using a digital micromirror device. The resulting reflected light was imaged with a camera. The images were then demodulated to extract the amplitude envelope for each spatial frequency measurement using an established amplitude demodulation algorithm.^42,43 A separate reference measurement at the same spatial frequencies was made on a calibration phantom with known optical properties for calibration of the SFDI source intensity and instrument response. This calibration enables the measured reflectance to be converted to absolute reflectance. This is achieved by comparing the measured reflectance from the calibration phantom with its predicted diffuse reflectance from an MC-based forward model using the phantom’s known optical properties. The region of interest (ROI) chosen for data analysis was centered in the detection FOV to avoid edge effects and measured $7.5\mathrm{cm}\times 2\mathrm{cm}$. The reported experimental data are average values taken over the ROI.

3.1.

Probability of Visitation and Detection $P_{\mathrm{V}\cap \mathrm{D}}$

We first consider a highly scattering tissue system with refractive index $n=1.4$ with optical properties providing $({\mu }_{\mathrm{s}}^{\prime }/{\mu }_{\mathrm{a}})=100$ and $l^{*}=1\mathrm{mm}$ and single-scattering anisotropy $g=0.8$. This corresponds to optical absorption and scattering coefficients ${\mu }_{\mathrm{a}}=0.00990099{\mathrm{mm}}^{-1}$ and ${\mu }_{\mathrm{s}}=4.95049505{\mathrm{mm}}^{-1}$, respectively. In the MC simulation, we utilize a narrow collimated beam normally incident into the tissue. Using the method described in Sec. 2.1, we launched $N={10}^8$ photon packets to obtain $P_{\mathrm{V}\cap \mathrm{D}}(z)$ with $z$ bins incremented at 0.01 mm and a set of spatial frequencies spanning 0 and $0.5{\mathrm{mm}}^{-1}$.

Figure 3(a) shows the $P_{\mathrm{V}\cap \mathrm{D}}(z)$ distributions obtained for spatial frequencies $f_x=0$, 0.025, 0.05, 0.075, 0.1, 0.125, 0.15, 0.175, 0.2, 0.250, 0.3, and $0.5\hspace{0.17em}\mathrm{m}{\mathrm{m}}^{-1}$. The value of the $P_{\mathrm{V}\cap \mathrm{D}}$ at the first surface $z=0\mathrm{mm}$, $P_{\mathrm{V}\cap \mathrm{D}}(z=0)$, is equivalent to the total diffuse reflectance for each spatial frequency, $R_{\mathrm{d}}(f_x)$. This is because every diffusely reflected photon packet will traverse this surface. The plots are shown with $1-\sigma$ error bars which indicate that 68% of independent simulations will produce results, which lie within this interval. The $f_x=0{\mathrm{mm}}^{-1}$ plot has a value of 0.62 at $z=0\mathrm{mm}$ representing total diffuse reflectance and decays to $4.7\times {10}^{-4}$ by depth $z=20\mathrm{mm}$, representing approximately a 3 order of magnitude reduction. The $f_x=0.5{\mathrm{mm}}^{-1}$ plot has a value of 0.03 at $z=0$ and decays to $3.6\times {10}^{-7}$ at $z=20\mathrm{mm}$, nearly a 6 order of magnitude reduction. These data illustrate the low-pass optical transport characteristics of scattering tissues.

Fig. 3

(a) $P_{\mathrm{V}\cap \mathrm{D}}(z)$ generated for media with optical properties ${\mu }_{\mathrm{s}}^{\prime }/{\mu }_{\mathrm{a}}=100$, $l^{*}=1\mathrm{mm}$ using $f_x=0,0.025,0.05,0.075,0.1,0.125,0.15,0.175,0.2,0.250,0.3,0.5{\mathrm{mm}}^{-1}$ with $1-\sigma$ error bars. Data were obtained at a $\mathrm{\Delta }z=0.01\mathrm{mm}$. Data symbols are shown on this plot at $z$-intervals of 1 mm. (b) $P_{z_{\max }}(z)$ derived from $P_{\mathrm{V}\cap \mathrm{D}}(z)$ distributions shown in (a).

Figure 3(b) shows the $P_{z_{\max }}(z)$ derived from $P_{\mathrm{V}\cap \mathrm{D}}(z)$ using Eq. (4). As described in Sec. 2.2, the benefit of transforming the $P_{\mathrm{V}\cap \mathrm{D}}(z)$ results to the max $z_i$ formulation $P_{z_{\max }}(z=z_i)$ is that (a) $P_{z_{\max }}(z=z_i)$, when integrated over all $z_i$, produces total diffuse reflectance $R_{\mathrm{d}}$ and (b) it effectively tracks the maximum depth sampled by each photon packet and provides group statistics on sampling depth for each spatial frequency.

3.2.

Experimental Validation

The experimental setup described in Sec. 2.3 consists of measurements taken from a two-layer phantom with a highly absorbing bottom layer placed at various depths $d$ that extinguishes any photons that propagate to that depth. The resulting measured reflectance is composed of only photons that never reach depths $z>d$, i.e., the photons detected possess trajectories with a maximum $z\le d$. The analogous computational result is given by Eq. (6).

Figure 4(a) shows the calibrated experimental measurements of diffuse reflectance versus top layer thickness $d$ and plots of $P_{z_{\max }}(z\le d)$ and (b) their difference. The plot shows a subset of the measured $f_x$ values for clarity. The depths $d$ and the spatial frequencies have been normalized to $l^{*}=1\mathrm{mm}$. For normalized spatial frequency $f_x{l}^{*}=0$, the plot rises from 0 and reaches a value of 0.595 for top layer thickness of $7.5\hspace{0.17em}d/l^{*}$. The plot rises monotonically as the top layer thickness increases because increasing numbers of photons fail to be extinguished by the bottom layer and are able to return to the surface to contribute to reflectance. As the spatial frequency increases, the measured diffuse reflectance flattens at a certain depth indicating that spatially modulated light for this frequency does not interrogate the tissue below that depth. For example, for $f_x{l}^{*}=0.3$, the spatially modulated reflectance rises from 0 and rises to 0.05 for top layer thickness of $1.0\hspace{0.17em}d/l^{*}$ without further increases for larger top layer thicknesses. The absolute difference between the experimental measurements and the computational predictions ranges between [$-0.012,0.025$].

Fig. 4

(a) Experimentally measured and calibrated $R_d$ (solid lines) and $P_{z_{\max }}(z\le d/l^{*})$ from simulation (dashed lines) and (b) their difference.

3.3.

Metrics for Optical Penetration Depth

We can determine the depth-dependent variation of the detected photon packets by determining the tissue depth that contains the complete photon packet trajectories corresponding to a certain fraction of the total measured diffuse reflectance at a given spatial frequency. For example, the spatial frequency dependence of sampling depth that encloses all the photon packet trajectories corresponding to only 50% of the detected reflectance can be determined by finding the value $d_{50}$ that results in a value of $X=0.5$ using Eq. (7). Similarly, by substituting alternate values for $X=0.1$, 0.25, 0.75, and 0.9 into Eq. (7), we computed depths $d_{10}$, $d_{25}$, $d_{75}$, and $d_{90}$ that correspond to the tissue depths that enclose the photon packet trajectories responsible for 10%, 25%, 75%, and 90% of the detected reflectance, respectively. These are shown in Fig. 5 with specific numerical values provided in Table 3. The span of the gray rectangles [25 to 75]% and vertical-capped lines [10 to 90]% provides range of depths sampled by these portions of the detected reflectance for these optical properties. In Appendix B, we provide similar plots and tables for a several ${\mu }_{\mathrm{s}}^{\prime }/{\mu }_{\mathrm{a}}$ values to show how these spans vary with optical properties.

Fig. 5

Median sampling depth ($d_{50}$) with [25 to 75]% (gray rectangle) and [10 to 90]% (vertical-capped line) intervals versus $f_x$ for media with optical properties ${\mu }_{\mathrm{s}}^{\prime }/{\mu }_{\mathrm{a}}=100$ and $l^{*}=1\mathrm{mm}$. Table of plot values are provided in Table 3.

These depth metrics are determined from the $z_{\max }$ distribution and provide a direct correspondence between the maximum tissue depths sampled by portions of the detected reflectance. Unlike the photon hitting density,⁶ these depths do not tally the pathlengths that the detected photon packet traverses within the tissue volumes considered and do not map directly to absorption sensitivity. Rather they provide tissue depths that correspond to portions of the detected reflectance. These allow the SFD user to determine what fraction of the measured reflectance is available to interact with the tissue beyond a certain depth. For example, if a given tissue has a $d_{90}=5\mathrm{mm}$, the user is certain that 90% of the measured reflectance is restricted to tissue depth $\le 5\mathrm{mm}$ whereas only 10% of the measured reflectance has the opportunity to sample tissue depth $>5\mathrm{mm}$.

3.4.

Sampling Depth Estimates for Real Tissue Types

The results provided in Fig. 5 can be generated for any set of optical properties. We apply our methodology to various tissue types to determine the depths that correspond to detection of 50%, 75%, and 90% of the total reflectance at a given spatial frequency. Table 1 shows the tissue types and wavelengths considered and the corresponding optical properties.^44–47 We performed independent MC simulations for each tissue type and wavelength pair to determine depth estimates for each. Figure 6 shows the predicted sampling depths versus spatial frequency for the various pairs of tissue type and wavelength considered.

Table 1

Tissue optical properties44–47 used in our sampling depth analysis.

Fig. 6

(a) $d_{50}$, (b) $d_{75}$, (c) $d_{90}$ sampling depths for human brain, mouse skin, human skin and human breast tissues at $\lambda =851\mathrm{nm}$ (solid lines) and 731 nm (dashed lines).

The tissue type, in which spatially modulated illumination can probe most deeply at the smaller spatial frequencies, is human breast at $\lambda =851$ and 731 nm. This is due to the higher ${\mu }_{\mathrm{s}}^{\prime }/{\mu }_{\mathrm{a}}$ properties for this tissue, 256 and 145, respectively, along with moderate transport mean-free path values of $l^{*}=0.89\mathrm{mm}$ and 1.04 mm, respectively. The median sampling depth of the human brain is 20% smaller than for human breast. While human brain tissue at $\lambda =731\mathrm{nm}$ has the equivalent ${\mu }_{\mathrm{s}}^{\prime }/{\mu }_{\mathrm{a}}$ value as human breast at $\lambda =851\mathrm{nm}$, the transport mean-free path in human brain is only $l^{*}=0.76\mathrm{mm}$ as compared to $l^{*}=1.04\mathrm{mm}$ in human breast. While mouse skin has the lowest ${\mu }_{\mathrm{s}}^{\prime }/{\mu }_{\mathrm{a}}$ values of the four tissue types considered, which would suggest more superficial optical sampling depths, the transport mean-free paths are the largest of the tissues considered resulting in optical sampling depths that are only slightly more superficial than human brain tissue at low spatial frequencies. The ${\mu }_{\mathrm{s}}^{\prime }/{\mu }_{\mathrm{a}}$ values for human skin are not much larger than mouse skin but with much higher scattering properties resulting in the smallest $l^{*}$ values of the group and the most superficial optical sampling depths.

The spatial frequency dependence characteristics of these optical sampling depths are also of interest. For spatial frequencies larger than $f_x=0.1{\mathrm{mm}}^{-1}$, differences in the optical sampling in human breast, human brain, and mouse skin are barely distinguishable. For human brain, the median depth values for $f_x=0.1{\mathrm{mm}}^{-1}$ are roughly half the median depth values using $f_x=0{\mathrm{mm}}^{-1}$ for both wavelengths. For mouse skin, the median depths at $f_x=0.1{\mathrm{mm}}^{-1}$ are about two-thirds the median depth values for $f_x=0{\mathrm{mm}}^{-1}$ for both wavelengths. The lower ${\mu }_{\mathrm{s}}^{\prime }/{\mu }_{\mathrm{a}}$ properties of both human skin and mouse skin are indicative of a diminished effect of scattering on the light transport and result in far less spatial frequency variation in the optical sampling depth. For human skin, the median depth at $f_x=0.1{\mathrm{mm}}^{-1}$ is 94% of the median depth value for $f_x=0{\mathrm{mm}}^{-1}$ using $\lambda =731\mathrm{nm}$ and this factor is 90% at $\lambda =851\mathrm{nm}$. For mouse skin, the median depth values at $f_x=0.1{\mathrm{mm}}^{-1}$ are about 82% the median depth values for $f_x=0{\mathrm{mm}}^{-1}$ for both wavelengths. At $f_x=0.5{\mathrm{mm}}^{-1}$, the median depth for all real tissue types is within the range [0.15 to 0.26] mm.

3.5.

Determination of Optical Sampling Depths for Other Tissue Types

To enable the determination of optical sampling depths for any tissue type, we used our methodology to determine the variation of optical sampling depth with spatial frequency for a range of ${\mu }_{\mathrm{s}}^{\prime }/{\mu }_{\mathrm{a}}$ values between 1 and 1000 while keeping $l^{*}=1\mathrm{mm}$ fixed and assuming $g=0.8$ and $n=1.4$. The specific optical properties considered are listed in Table 2. Figure 7 shows the median depth of optical sampling as a function of ${\mu }_{\mathrm{s}}^{\prime }/{\mu }_{\mathrm{a}}$. In this figure, both the median sampling depth and the spatial frequency of illumination are scaled relative to $l^{*}$. The plots show that as the ${\mu }_{\mathrm{s}}^{\prime }/{\mu }_{\mathrm{a}}$ increases, so does the median depth. Moreover, at larger spatial frequencies, there is less sensitivity of the sampling depth to variations in ${\mu }_{\mathrm{s}}^{\prime }/{\mu }_{\mathrm{a}}$.

Table 2

General optical properties with l*=1  mm used for our sampling depth lookup table.

Fig. 7

Median sampling depth $d_{50}$ normalized by the transport mean-free path $l^{*}$ as a function of ${\mu }_{\mathrm{s}}^{\prime }/{\mu }_{\mathrm{a}}$ at spatial frequencies $f_x{l}^{*}=0$, 0.05, 0.1, 0.2, and 0.3. For clarity, results for only a subset of $f_x$ values are plotted here. Results for additional $f_x$ values are provided in the supplemental data.

Optical sampling depths determined using the general optical properties can be composed into a lookup table, then scaled, and interpolated to provide depth statistics for an arbitrary tissue. Only knowledge of the tissue absorption and reduced scattering properties is needed. The details of how this is performed are given in Appendix B. The advantage of this lookup table method is that it dispenses with the need to execute an MC simulation for the specific tissue optical properties in question and enables rapid estimation of SFD sampling depths. We performed the lookup table method using the ${\mu }_{\mathrm{a}}$ and ${\mu }_{\mathrm{s}}^{\prime }$ values of the real tissue optical properties listed in Table 1. Figure 8 shows the relative differences between the median depth $d_{\mathrm{table}}$ determined by the lookup table method and the median depth determined by the independent MC simulation at the real tissue optical properties $d_{\mathrm{MC}}$, $(d_{\text{table}}-d_{\mathrm{MC}})/d_{\mathrm{MC}}$. The lookup table estimates and the independent MC simulation results agree to within 7%, suggesting that the lookup table provides an accurate and convenient method for determining depth predictions.

Fig. 8

Relative difference between the median depth determined by the lookup table of general optical properties and the median depth determined by running an MC simulation using the real tissue optical properties.