Development of perturbation Monte Carlo methods for polarized light transport in a discrete particle scattering model

Jennifer Nguyen; Carole K. Hayakawa; Judith R. Mourant; Vasan Venugopalan; Jerome Spanier

doi:10.1364/BOE.7.002051

1. Introduction

The study of in vivo and in vitro cell and tissue microstructure, which yields information about underlying cellular processes and responses, has been historically addressed by optical biopsy techniques. Indeed, there appears to be growing interest in the use of polarization-sensitive optical probes as a tool to non-invasively obtain information regarding subcellular microstructure. Mourant and co-workers have integrated polarized wavelength-dependent measurement capabilities into reflectance based optical probes to increase sensitivity to scatterer size [1]. Backman and co-workers have employed polarization-gating techniques to characterize and control the depth of the probed tissue volume [2]. Sokolov and co-workers have used a probe with angled illumination collection geometry that detects light polarization to obtain depth-specific information of scatterer size and size distribution [3]. These methods, however, all represent instrumentation-based approaches to extract information about the underlying subcellular structure. Here we propose a versatile and rigorous computational method to analyze polarized signals in order to provide insight into subcellular morphology.

The perturbation Monte Carlo (pMC) method is well-established and has been applied to a wide range of biomedical optics problems including: determination of layered tissue optical properties [4, 5], tissue image reconstruction [6], optical tomography [7], and time-resolved functional imaging [8, 9]. Here we provide a novel extension to pMC that builds upon our previous work [10] to characterize the utility of pMC in providing estimates of polarization-sensitive reflectance measurements. Our exploration in this study also serves as a demonstration of one of many applications that can benefit from the use of pMC and as a mechanism for discovering the factors that affect pMC performance for a class of similar problems.

2. Model problem

2.1. Light scattering model

We develop a pMC computational framework within a tissue representation that models optical scattering using three log-normal distributions to represent the polydispersity of the scattering medium [1,11–13]. The pMC method introduced here can accommodate any number (m) of scatterer groups with arbitrary distributions of scatterer size. The scattering coefficient of the ith distribution, μ_s,i, is calculated using $μ_{s, i} = \int Q_{s c a t} (t) N_{a, i} (t) L_{i} ({\bar{r}}_{i}, t) d t$ where Q_scat is the scattering efficiency, N_a,i is the product of number density and the scattering cross-section of the particle, t is a particle radius, and L_i is the ith log-normal distribution characterized by its mean radius ${\bar{r}}_{i}$ and standard deviation, $σ_{i} : L_{i} ({\bar{r}}_{i}, t) = \exp {- {(\ln t - \ln {\bar{r}}_{i})}^{2} / (2 σ_{i}^{2})} / (t σ_{i} \sqrt{2 π})$ . The scattering coefficient of the medium ${\bar{μ}}_{s}$ can be calculated by summing the contribution of each distribution to the scattering coefficient ${\bar{μ}}_{s} = \sum_{i = 1}^{m} μ_{s, i}$ [14]. The phase function of each scatterer population is calculated by integrating across radii values and using the log-normal distribution to properly weight the phase functions that Mie Theory provides: f_i(θ) = 1/μ_s,i ∫ f_Mie(t, θ)Q_scat(t)N_a,i(t)dt where f_i is the phase function of the ith distribution of scatterers, f_Mie is the phase function of a mono-dispersed distribution of scatterers, and θ is the scattering angle. We approximate these integrals for μ_s,i and f_i using trapezoid rule integration with uniform grid spacing with the width of the bin set at 0.001 µm for all three distributions.

Table 1 provides scattering parameters that approximate the contributions of protein complexes, organelles, and nuclei present in epithelial cervical tissues based on previous experimental results [1,11]. These parameter values serve as baseline values in the pMC results presented this paper and are meant to model a class of problems (similar to the problem in [11]) in which pMC application could potentially benefit forward and inverse solutions.

Table 1. Distributions of scatter sizes in the tissue model

View Table | View all tables in this article

2.2. Measurement geometry

In this study we consider a probe with two detectors placed equidistant from a central source fiber which detect either unpolarized or polarized reflectance. A probe schematic is shown in Fig. 1(c). All detector fibers are angled at ω =20° relative to the outward pointing surface normal and towards the source fiber. The distance from the center of any detector fiber to the center of the source fiber is 550 µm and all the fibers are 480 µm in diameter. All fibers have a light cone detection half angle of 21.7°. The close source-detector separation and the angling of the detector fibers towards the source fibers promotes photon trajectories that interact with scatterers located at shallow depths.

Fig. 1 A conceptual diagram of the pMC application to the light scattering model used in this study. (a) Depiction of the light scattering model composed of three distinct distributions of scatterers. The dotted blue and red plots represent perturbations in mean radius and number density while solid green, blue and red plots represent “baseline” parameter values. (b) A schematic diagram of the source-detector configuration. (c) The top-down view of the probe. The yellow circle represents the source and the circles with numbers inside represent the two detectors. See Section 2.2 for probe details.

Download Full Size | PDF

3. Theory

The radiative transport equation (RTE) is an integro-differential equation describing the transport of photons in a turbid medium that can be transformed to an integral equation:

Ψ = K Ψ + Q

where

K Ψ (r, Ω) = \int_{S^{2}} \int_{R} K (r^{'}, Ω^{'} \to r, Ω) Ψ (r^{'}, Ω^{'}) d r^{'} d Ω^{'}

, Ψ is the radiance or transport solution, the function Q models some photon source, R is the set of all possible positions in space, and S² is the set of all solid angles. In this paper,

K

is an integral operator, K is the RTE kernel that characterizes light transport from a location r′ and unit direction Ω′ to another location r and unit direction Ω. To predict the detection of photonic signals, we are typically interested in estimating one or more integrals I with the general form

I = \int_{S^{2}} \int_{R} g (r, Ω) Ψ (r, Ω) d r d Ω

where g is a characteristic detector function and I is some value of interest, often reflectance or transmittance. A rigorous theory that exhibits the equivalence between an RTE-based analytic model, as described by Eqs. (1) and (2), and a Monte Carlo-based stochastic model is treated in detail in [10,15,16]. That equivalence produces the string of equations that relates the stochastic model with the physical model

E [ξ] = \int_{ℬ} ξ (b) d M (b) = \lim_{N \to \infty} \frac{1}{N} \sum_{n = 1}^{\infty} ξ (b_{n}) \int_{S^{2}} \int_{R} g (r, Ω) Ψ (r, Ω) d r d Ω

where

ξ : ℬ \to ℛ

is a random variable (estimator) on the sample space

ℬ

of biographies and M is a probability measure on

ℬ

. Roughly speaking, the measure M describes a probability density function defined on

ℬ

that can be interpreted as a “likelihood function” for each

b \in ℬ

. This likelihood function is uniquely induced by the specific functions used to generate the biographies in

ℬ

(details can be found in [15]). In Eq. (3)E is the expected value of ξ with respect to the measure induced on

ℬ

by the random walk process used to generate biographies, the second equality in Eq. (3) states that ξ is an unbiased estimator of E[ξ], and the third equality in Eq. (3) demonstrates the equivalence between the analytic model and the stochastic model. This equivalence is necessary to support the derivation of the pMC model.

3.1. Perturbation as change of measure

Suppose that $ℳ$ and $\hat{ℳ}$ are two different probability measures that are induced on $ℬ$ by two different methods of generating photon biographies $b \in ℬ$ . For example, $ℳ$ could represent the baseline measure (which for us is generated by the continuous ∈ absorption weighting (CAW) process [16]), and $\hat{ℳ}$ is a different measure representing any perturbation of $ℳ$ . If we denote $ℳ_{CAW}$ as the CAW-induced measure, the measure ${\hat{ℳ}}_{CAW}$ in our case would arise by replacing the baseline optical properties by any perturbation of them. The Radon-Nikodym Derivative, $\frac{d {\hat{ℳ}}_{CAW}}{d ℳ_{CAW}}$ , may be used to estimate the reflectance obtained when ${\hat{ℳ}}_{CAW}$ replaces $ℳ_{CAW}$ . The existence of the Radon-Nikodym Derivative can be guaranteed only if certain technical conditions are satisfied, namely, the sets of biographies that are “small” in the ${\hat{ℳ}}_{CAW}$ measure are also small in the $ℳ_{CAW}$ measure. A detected photon measure can also be expressed as an expected value with respect to the modified measure $\hat{ℳ}$ of a modified random variable $\hat{ξ}$ , where

\hat{ξ} = ξ \frac{d {\hat{ℳ}}_{CAW}}{d ℳ_{CAW}} so that \int \hat{ξ} d ℳ_{CAW} = \int ξ \frac{d {\hat{ℳ}}_{CAW}}{d ℳ_{CAW}} d ℳ_{CAW} = \int ξ d {\hat{ℳ}}_{CAW} .

We may interpret the application of the Radon-Nikodym derivative in the context of pMC as re-weighting the baseline collected photon weights to account for the change in the density functions used to generate the random walks.

3.2. Perturbation Monte Carlo

The pMC method is effective for estimating detected photon weights from small changes in parameters that characterize the forward RTE. When using pMC, a set of Monte Carlo biographies is generated to produce a baseline dataset consisting of the location, weight and incoming unit direction at every collision point for each collected photon. These biographies are then post-processed using pMC equations to determine the effect of perturbations in tissue composition on the detected photon weights. The advantage of this method is that it requires only a single set of Monte Carlo biographies, namely those in the baseline set, to perform simulations for tissues with perturbed properties using only a post-processing module. The computational expense of the post-processing module is 1-3 orders of magnitude smaller compared to that needed to generate a complete set of independent MC biographies.

The CAW random walk process is employed in all MC simulations in this study. Using CAW the photon is deweighted by the factor exp(−μ_as), where s is the pathlength between collisions and intercollision distances are sampled from ${\bar{μ}}_{s} \exp (- {\bar{μ}}_{s} s)$ . In addition to propagation and absorption, a photon may scatter through a given scattering angle as specified by the phase function. Our study uses a Distribution Selection Scattering Method (DSSM) where at each scattering interaction a photon has a probability of $μ_{s, i} / {\bar{μ}}_{s}$ to interact with particles of the ith scattering distribution. Our previous paper [10] uses the Composite Phase Function Scattering Method (CPFSM). A proof of the equivalence between CPFSM and DSSM is provided in the Appendix. Once the scattering distribution has been selected, the phase function corresponding to that distribution, f_i(θ), will be used to sample the scattering angle, θ. This results in the following RTE kernel [16]:

K = \frac{μ_{s, i}}{{\bar{μ}}_{s}} f_{i} (Ω^{'} \to Ω) \exp (- μ_{a} s) {\bar{μ}}_{s} \exp (- {\bar{μ}}_{s} s) .

In Eq. (5), $μ_{s, i} / {\bar{μ}}_{s}$ is the probability of interacting with the ith distribution through the DSSM, f_i(Ω′→ Ω) is the probability density that takes the photon from a direction of Ω′ to Ω for the ith distribution of scatterers, and $\exp (- μ_{a} s) {\bar{μ}}_{s} \exp (- {\bar{μ}}_{s} s)$ accounts for scattering probability density and the absorption attenuation associated with a photon pathlength s. This equation assumes a homogeneous absorption coefficient throughout the medium. Application of pMC through the Radon-Nikodym derivative using the kernel in Eq. (5) produces the perturbed estimator $\hat{ξ}$ :

\begin{matrix} \hat{ξ} = ξ {(\frac{{\hat{μ}}_{s, 1}}{μ_{s, 1}})}^{j_{1}} {(\frac{{\hat{μ}}_{s, 2}}{μ_{s, 2}})}^{j_{2}} \dots {(\frac{{\hat{μ}}_{s, m}}{μ_{s, m}})}^{j_{m}} \exp [- ({\hat{\bar{μ}}}_{s} - {\bar{μ}}_{s}) L] \exp [- ({\hat{μ}}_{a} - μ_{a}) L] \\ \prod_{l_{1} = 1}^{j_{1}} \frac{\hat{f} (θ_{l_{1}}, ϕ_{l_{1}})}{f (θ_{l_{1}}, ϕ_{l_{1}})} \prod_{l_{2} = 1}^{j_{2}} \frac{\hat{f} (θ_{l_{2}}, ϕ_{l_{2}})}{f (θ_{l_{2}}, ϕ_{l_{2}})} \dots \prod_{l_{m} = 1}^{j_{1}} \frac{\hat{f} (θ_{l_{m}}, ϕ_{l_{m}})}{f (θ_{l_{m}}, ϕ_{l_{m}})} \end{matrix}

where j₁ is the total number of collisions experienced by the photon with respect to the first scattering distribution, j₂ is the total number of collisions experienced by the photon with respect to the second scattering distribution, …, j_m is the total number of collisions experienced by the photon with respect to the last (the mth) scattering distribution of the light scattering model, L is the total pathlength traveled by the photon in the medium, and

θ_{l_{1}}

and

ϕ_{l_{1}}

are polar and azimuthal scattering angles, respectively, of the lth collision that result from interacting with the 1st distribution. The estimator ξ refers to the collected photon weight from the baseline cMC simulation, which can be polarized parallel or perpendicular to the source. All unhatted quantities refer to baseline parameter values for the quantity and hatted quantities refer to perturbed parameter values.

\hat{{\bar{μ}}_{s}}

is the total scattering coefficient of the perturbed medium,

{\hat{μ}}_{s}

is the absorption coefficient of the perturbed medium, and

{\hat{f}}_{i}

is the perturbed phase function of the ith distribution. The sample mean of

\hat{ξ}

based on N photon biographies can be used as an estimator for the reflectance and is denoted as

{\hat{\bar{ξ}}}_{N}

. The standard deviation of the sample mean is denoted by

σ [{\hat{\bar{ξ}}}_{N}]

.

4. Methods

To test pMC performance, pMC reflectance estimates were compared to those obtained using cMC. cMC reflectance estimates were obtained by launching 20 million photons in a tissue representation characterized by a combination of the parameter values in Table 2 and parameter baseline values in Table 1. The ranges in Table 2 were chosen with the intent of modeling the length scales over which the index of refraction changes in epithelial cells. pMC estimates were obtained by launching 20 million photons in a single Monte Carlo simulation at the baseline parameter values in Table 1 and calculating the sample averages of the estimator $\hat{ξ}$ to obtain perturbed reflectance values. Error bars shown on the plots represent the standard deviation of the sample means. All photon biographies were obtained from a modified version of the conventional Monte Carlo simulation code outlined in [17].

Table 2. Range of parameters examined in pMC calculations. Baseline ${\bar{μ}}_{s} = 125 {cm}^{- 1}$ and $\bar{g} = 0.9534$ .

View Table | View all tables in this article

cMC and pMC calculations were performed on the High Performance Cluster at University of California, Irvine and a single private computer, respectively. The cluster runs CentOS 6.6 and code on the cluster was compiled using gcc−4.8.2. The private computer has 2 × 1.33 GHz Quad-Core AMD processors and 8 GB 1333 MHz DDR2 RAM running Ubuntu 10.04.4 LTS and our code was compiled using gcc−4.1.3. Benchmark tests were performed to ensure similar performance across the two settings.

Both polarized and unpolarized simulations were performed using the same optical properties and probe geometry. Each cMC unpolarized light simulation launching 20 million photons took 22 – 45 minutes to run while the cMC simulations for polarized light transport took 6.7–18.9 hours. The post-processing module generating pMC reflectance estimates requires only ∼7 seconds and ∼13 seconds for unpolarized and polarized cases, respectively. Table 2 provides the perturbations used for the pMC tests we performed by applying perturbations in mean radii $\hat{{\bar{r}}_{i}}$ and number densities $\hat{N_{s, i}}$ . Our model has three distributions of scatterers and each distribution has an associated anisotropy factor. We refer to the anisotropy factor of the ith distribution as g_i and the ensemble anisotropy factor of the medium as $\bar{g}$ . Figure 2 shows how μ_s and g change as these parameters are changed with the parameter perturbation.

Fig. 2 (a) Changes in the mean radii and how it relates to g_i. (b) Changes in the ensemble scattering coefficient for all perturbations.

Download Full Size | PDF

In both cMC and pMC simulations, we treat the source fiber center as the origin of our coordinate system. The center of Detector 1 in Fig. 1(c) resides along the y-axis and the center of Detector 2 resides along the x-axis of the coordinate system. Our pMC implementation requires the biographies of detected photons generated by an initial baseline cMC simulation.

In the case of polarized light transport, polarization is tracked by rotating the Stokes vector as outlined in [18]. All photons initially enter the medium linearly polarized and with polarization parallel to the x-axis. Each scattering event is characterized by polar and azimuthal angles, θ and ϕ, which are sampled from the phase function through a rejection sampling method used by Bartel and Hielscher [19]. Once the photon exits the medium, the Stokes vector is used to calculate the components of the photon weight that are, respectively, parallel and perpendicular to the x-axis.

5. Results

5.1. Scattering model parameter perturbations

First, we investigate the ability of our pMC method to predict changes in reflectance produced by perturbations in the mean radii and/or number density of one of the three particle populations. Figure 3 shows a comparison between cMC and pMC reflectance estimates for Detector 1, which tallies photon weights with parallel or perpendicular polarization relative to the source. Figs. 3(a), 3(c), and 3(e) provide a comparison of pMC and cMC reflectance estimates resulting from perturbations in number density for each of the three scattering particle populations, whereas Figs. 3(b), 3(d), and 3(f) compare pMC and cMC reflectance estimates for perturbations in mean radii.

Fig. 3 pMC and cMC estimates of reflectance of parallel and perpendicular polarization in Detector 1 for perturbations in number density ((a), (c), and (e)) and mean radii ((b), (d), and (f)). Solid lines display the trend for parallel reflectance estimates and dashed lines display the trend for perpendicular reflectance estimates.

Download Full Size | PDF

The agreement between pMC and cMC estimates in Fig. 3 is linked to the range in which $\hat{{\bar{μ}}_{s}}$ is perturbed and the degree in which ${\hat{g}}_{i}$ is perturbed. For perturbations of the scattering model parameters, the two largest ranges in $\hat{{\bar{μ}}_{s}}$ occur when the second and third mean radii are perturbed, as seen in Table 2. Although perturbations in the second and third mean radii produce similar $\hat{{\bar{μ}}_{s}}$ ranges, the perturbations in the second mean radii coincides with a larger perturbation in anisotropy factor which causes a more pronounced disagreement between pMC and cMC. As shown in Table 2, we explored a perturbation range for mean radius that spans a larger range of ensemble $\hat{{\bar{μ}}_{s}}$ values as compared to perturbations in the number density. For example, perturbations we examined for number density of the first particle distribution ${\hat{N}}_{s, 1}$ result in a range of ensemble scattering coefficients $\hat{{\bar{μ}}_{s}} = 123 - 127 {cm}^{- 1}$ , whereas perturbations we examined for first mean radius in the range $\hat{{\bar{r}}_{1}}$ result in a broader range of $\hat{{\bar{μ}}_{s}} = 121 - 152 {cm}^{- 1}$ . Similar observations can be made for the second and third populations. Because number density perturbations do not result in perturbation of the phase function and because of the smaller changes in $\hat{{\bar{μ}}_{s}}$ , these pMC reflectance estimates have smaller variance than when perturbing mean radii. From this, we conclude that pMC performance is a function of the magnitude of $\hat{{\bar{μ}}_{s}}$ changes as well as the changes in the phase function which can be characterized by the change in g.

Interestingly, Fig. 3 reveals that reflectance is highly sensitive to changes in the first distribution since the slopes of the plots in Figs. 3(a) and 3(b) are much greater than those in Fig. 3(c)–(f). The sensitivity to the parameters of the first distribution can be explained by the low anisotropy factor of the associated phase function, which heavily promotes backscattering (See Table 2). This probe’s bias toward detecting photons that have undergone backscattering is connected to the geometry of the probe which features angled detector fibers and a small source-detector separation. Based on $μ_{s, i} / {\bar{μ}}_{s}$ ratios, we would expect that 3.0%, 49.8%, and 47.2% of all scattering interactions would be with the first, second, and third distributions, respectively. However, analysis of the collected photons reveal that 4.7%, 48.8%, 46.4% of the scattering interactions are with the particles from the first, second, and third distribution, respectively. This represents a 54% increase in the number of interactions that the detected photons have with particles from the first population over the expected number of interactions. This is evidence that collected photons have enriched sensitivity toward changes to parameters associated with the first particle distribution.

Another important observation from the results in Fig. 3 is that the polarization-sensitive measurements provide sensitivity to parameters of the second particle distribution. This result differs from findings in [10] where we saw sensitivity to only parameters of the first distribution with unpolarized reflectance estimates. Changes in the first and second scattering populations, which represent protein complexes and organelles, are associated with dysplasia [11] and sensitivity to parameters of both first and second scattering populations may be helpful. Figure 4 provides a comparison between pMC reflectance estimates for perturbations in the second mean radii for the cases of unpolarized and polarized light propagation and is intended to supplement the results in Fig. 3(d). In Fig. 4(a) and Fig. 4(b), the change in reflectance across $\hat{{\bar{r}}_{2}}$ values are observable where as in Fig. 4(c), the changes in reflectance are not discernible from the noise in pMC and cMC reflectance estimates.

Fig. 4 pMC estimates for perturbation in second mean radius in cases of (a) parallel polarization (b) perpendicular polarization and (c) unpolarized detection.

Download Full Size | PDF

Although pMC’s reflectance estimates are particularly poor for Detector 2 with perpendicular polarization, pMC accurately captures the 6% change in reflectance at $\hat{{\bar{r}}_{2}} = 0.41 - 0.47 μ m$ for Detector 2 with perpendicular polarization. The differences between pMC reflectance estimates in the two detectors and polarizations will be further addressed in §5.2.

In Fig. 5, we explore the ability of pMC to reproduce spectral reflectance measurements. The baseline wavelength is λ = 620 nm and the wavelength perturbations are performed in the range of λ = 500–720 nm. Alteration of the wavelength requires modification of both the scattering coefficients and phase functions for all three scatterer distributions. The spectral asymmetry in Fig. 5 of pMC’s performance relative to the baseline wavelength can be explained by the more rapid increase of the scattering coefficient as the wavelength decreases, as shown in Table 2. We also see that larger wavelength perturbations away from the baseline value of λ = 620nm result in larger standard deviations of the resulting pMC estimates. Perturbations from the baseline wavelength result in perturbations in both the ensemble scattering coefficient and the anisotropy factors, both of which drive increases in the variance of the estimator $\hat{ξ}$ .

Fig. 5 pMC and cMC estimates of reflectance across wavelength perturbations for polarized light propagation.

Download Full Size | PDF

To summarize, we find that pMC’s performance is negatively affected by: (1) increasing perturbations in the scattering coefficient and (2) extreme changes in the phase function which can be characterized by anisotropy factor perturbation. We have also demonstrated the utility of our method by showing that turning on polarization can add sensitivity to some parameters and by showing the accuracy of our method in the problem of perturbing across wavelength.

5.2. PMC performance and conventional Monte Carlo relative error

Second, we examine the relative performance of pMC reflectance estimates for parallel and perpendicular polarization as shown in Figs. 3 – 5. The top part of Table 3 shows the parameter ranges in which pMC estimates are within 5% of cMC reflectance estimates. The elements in the top part of Table 3 have a lower bound followed by an upper bound for the parameter values. These bounds are recorded as percentage changes from the baseline value for that parameter. For example, in the column Detector 1, ║ and in the row r₁, the parameter range is listed as −50%, +40%. This means that pMC reflectance estimates for Detector 1 with parallel polarization were within 5% of cMC reflectance estimates for the range r₁ = 0.015 – 0.042 µm since the baseline value for r₁ is 0.03 µm. In the bottom part of Table 3, we compare the relative errors, which is defined as the standard deviation of the estimate divided by the mean, of pMC estimates at baseline parameter values.

Table 3. Relative error and parameter ranges for polarized pMC reflectance estimates that are within 5% of cMC reflectance estimates for each detector.

View Table | View all tables in this article

The parameter ranges for accurate pMC estimates are governed by (a) the inherent noise of the baseline cMC simulation and (b) the size of the perturbation characterized by the changes in $\hat{{\bar{μ}}_{s}}$ , ${\hat{g}}_{i}$ , and ${\hat{μ}}_{s, i}$ . The effect of the inherent noise of the baseline cMC simulation can be determined through the comparison of the pMC reflectance estimates between Detector 1, ║ and Detector 1, ⊥. This is because reflectance estimates use the exact same photon biographies and the differences between pMC performance can be attributed to the final weights. The fact that Detector 1, ⊥ provides accurate pMC reflectance estimates over smaller parameter ranges than Detector 1, ║ can be explained by the larger relative error in the baseline cMC simulation for Detector 1, ⊥. Similar arguments can be made for the range of accuracy of the pMC estimates for Detector 2, ║ and Detector 2, ⊥. The accuracy range reported for N_s,₁ in Table 3 is the full range of parameter values for N_s,₁ in this study; the true limit of pMC’s performance lies beyond the range listed for N_s,₁.

The observed disparity between the relative errors of reflectance estimates for Detectors 1 and 2 is consistent with prior experimental work [11, 20, 21] and is consistent with the underlying physics of our model problem, which causes Detector 1 to detect a larger reflectance polarized parallel to the x-axis than Detector 2. This asymmetry in photon collection can be explained by a simple model outlined in [21] for a similar optical probe geometry. This model uses three assumptions: (1) photons undergo only 2 scattering interactions before entering a detector, (2) photons that are collected by Detector 1 have scattering events in the Y-Z plane and photons that are collected by Detector 2 have scattering events in the X-Z plane and (3) photons enter the medium with a trajectory parallel to the z-axis. Figure 6 shows plots of |s₁(θ)|² and (s₂θ)|², where s₁ and s₂ are elements of the scattering amplitude matrix yielded by Mie Theory. The plot of |s₁(θ)|² (blue) is the scattered irradiance per unit incident irradiance given that the incident light is polarized perpendicular to the scattering plane and the plot of |s₂(θ)|² (red) is the scattered irradiance per unit incident irradiance given that the incident light is polarized parallel to the scattering plane [22]. Since each photon initially starts with polarization parallel to the x-axis, |s₁(θ)|² can be interpreted as the phase function of scattering events in the Y-Z plane and |s₂(θ)|² can be interpreted as the phase function of scattering events in the X-Z plane in the context of this three assumption model. The greatest difference between |s₁(θ)|² and |s₂(θ)|² occurs at a scattering angle of around 90°. In the event that a photon scatters at two 90° angles before entering the detector, then the implication of Fig. 6 is that a photon will end up in the Y-Z plane more frequently.

Fig. 6 A plot of the s₁ and s₂ components produced by the Mie Scattering Method.

Download Full Size | PDF

In this section, we have shown that there is a link between the noise of the baseline cMC simulation and the subsequent pMC performance. This link is helpful in understanding the difference in performance between the pMC reflectance estimates for parallel and perpendicular polarizations. Furthermore, to ensure better pMC performance, it is critical to keep the relative error of the baseline cMC simulation as low as possible.

5.3. PMC and cMC computational efficiency comparisons

To examine the computational advantage provided by pMC over cMC, we examine their relative computational efficiencies, η, defined as

η = \frac{{\bar{R}}^{2}}{σ^{2} T}

where

\bar{R}

is the mean of the estimate, σ² is the variance, and T is the time required for the calculation [23]. This figure of merit allows for both characterization of the performance of the estimator and the amount of computational resources used by combining the relative error of the estimator and the computational time. Furthermore, this figure of merit is dimensionless since relative error is proportional to

1 / \sqrt{N}

where N is the number of photons launched and T is proportional to N; this lack of dimensionality then allows for comparison across different algorithms or estimators. Figs. 7(a) and 7(b) plots the computational efficiency for perturbations in wavelength, mean radii, and number densities. Both plots feature a relatively constant computational efficiency value for cMC estimates, while pMC estimates have the highest computational efficiency when the change in the ensemble scattering coefficient is zero. The gain in computational efficiency for pMC estimates degrades as the perturbation in the ensemble scattering coefficient increases, which results from an increase in standard deviation of the estimates. In Fig. 7(a), plots of the computational efficiency are shown for perturbations in the mean radii whereas in Fig. 7(b), plots of the computational efficiency are shown for perturbations in the number density. A comparison of the pMC computational efficiency plots in Fig. 7(a) and 7(b), reveals that changes in

\hat{{\bar{r}}_{1}}

degrade the computational efficiency even though the range in

\hat{{\bar{μ}}_{s}}

is much smaller than perturbations for

\hat{{\bar{r}}_{2}}

and

\hat{{\bar{r}}_{3}}

. This degradation is the result of increases in the size of the standard deviation of the estimate, which is caused by the perturbation in anisotropy factor of the phase function of the first population.

Fig. 7 A comparison of the computational efficiency of pMC and cMC estimates for (a) perturbations in $\hat{{\bar{r}}_{i}}$ and (b) perturbations in ${\hat{N}}_{s, i}$ . The pink and red symbols show the computational efficiency of pMC estimates. The blue and light blue symbols in the plot above are for the computational efficiencies of cMC estimates.

Download Full Size | PDF

Recall that pMC is a true perturbation method; that is it works best if the perturbations are “small”. The pMC estimator is unbiased for ALL perturbations, but the larger the perturbation, the larger the statistical error it incurs. (When we say that the pMC estimator is unbiased we mean that on average the estimator will not tend to over-estimate or under-estimate the true value of the parameter.) Thus, in the limit of very large perturbations, it is more efficient computationally to examine independent Monte Carlo simulations of the baseline pMC and at the perturbed values of parameters of interest. As the size of the perturbation shrinks to zero, the advantage of pMC over independent conventional simulations tends to increase.

6. Summary and conclusion

We have developed a new pMC algorithm that provides estimates of linearly polarized reflectance resulting from linearly polarized incident light. Unlike most other pMC algorithms, this algorithm takes into account phase function perturbations due to varying scatterer parameters values or incident wavelength.

The performance of this implementation of the pMC method depends primarily on the magnitude of the perturbation of the anisotropy factor and the ensemble scattering coefficient of the medium. We have applied the algorithm to the problem of accurately solving a forward problem for a scattering model that simulates realistic tissue optical properties and have shown that the pMC method is able to produce estimates that agree in a probabilistic sense with cMC estimates for a limited range of s and g perturbations.

The ability to model polarized light transport is important because polarized light transport measurements can sometimes provide information not available from unpolarized measurements. We have previously shown that ratios of measurements made by this type of probe can differentiate changes in size from changes in concentration in monodisperse populations of scatterers [17] using single wavelength measurements. The scattering media modeled in this paper are more complex. We demonstrate that polarized measurements provide sensitivity to specific morphological parameters (r₂, N_s,₂), which is not possible with unpolarized measurements. The sensitivity provides the opportunities to resolve inverse problems that may be out of reach when acquiring measurements of unpolarized light. Consequently, light transport parameters may be obtainable from more complicated systems.

7. Appendix

7.1. Scattering algorithms

Two distinct scattering algorithms were considered for this study: the Composite Phase Function Scattering Method and the Distribution Selection Scattering Method. We will briefly explain the algorithm behind each of these methods and show that these algorithms produce equivalent distributions of scattering angles. We chose DSSM in this paper because the derivative of the resulting estimator is computationally simple to obtain.

The Composite Phase Function Scattering Method used in [10] constructs a composite phase function by calculating the weighted average of the phase functions of each distribution of scatterers, where the weights are based on each distribution’s contribution to the scattering coefficient

\bar{f} = \sum_{i = 1}^{m} \frac{μ_{s, i}}{{\bar{μ}}_{s}} f_{i}

where

\bar{f}

is the composite phase function, f_i is the phase function of the ith distribution, μ_s,i is ith scatterer distribution’s contribution to the scattering coefficient and

{\bar{μ}}_{s}

is the scattering coefficient of the entire medium and is the sum of the scattering coefficient contributions of all distributions.

The Distribution Selection Scattering Method used in this study, calculates the probability of interacting separately with each of the scatterer distributions. The probability of interacting with the ith scattering distribution is proportional to that distribution’s contribution to the scattering coefficient, i.e., $P (Y = i) = μ_{s, i} / {\bar{μ}}_{s}$ , where Y is a random variable that selects the population for the scattering event. Given that the photon must interact with the ith scattering distribution, the phase function of the ith scattering distribution, f_i, is then sampled for a scattering angle. The composite scattering phase function is not constructed in this method.

7.2. Proof of scattering algorithm equivalence

We will demonstrate equivalence between the Composite Phase Function Scattering Method and the Distribution Selection Scattering method by showing that both methods ultimately produce the same cumulative distribution function (CDF).

Consider a medium that has m groups of scatterers that each contribute to the medium’s scattering coefficient. Recall that in the Distribution Selection Scattering Method, a distribution is randomly selected according to the percentage of the population of that scatterer to the total scatterer population. The variable Y is defined as a random discrete variable that selects the population for some scattering event and can take on integer values between 1 and m. The CDF of Y is

F_{Y} (y) = {\begin{matrix} 0, & y < 1 \\ \sum_{k = 1}^{f l o o r (y)} \frac{μ_{s, k}}{{\bar{μ}}_{s}}, & 1 \leq y \leq m \end{matrix} .

Given that Y selects the ith distribution, its phase function will be used to select a polar scattering angle θ. The conditional CDF of θ, given some value of Y is the CDF of that distribution’s phase function

F_{Θ Y} (θ | Y = i) = \int f_{i} \sin θ d θ .

Since the events that Y = 1, Y = 2, …, Y = m are disjoint, we can obtain the CDF for θ for DSSM through an application of Bayes’ Theorem using the probabilities in Eq. (9) and the conditional CDF in Eq. (10)

F_{n e w} (θ) = \sum_{i = 1}^{m} F_{Θ Y} (θ | Y = i) P (Y = i) = \sum_{i = 1}^{m} \int f_{i} \sin θ d θ {\frac{μ_{s, i}}{{\bar{μ}}_{s}}} = \int \sum_{i = 1}^{m} f_{i} \sin θ d θ {\frac{μ_{s, i}}{{\bar{μ}}_{s}}} .

In the last step above, we reverse the order of the summation and the integral so that Eq. (11) will more closely resemble Eq. (12). We invoke Tonelli’s theorem, which states that if f_n(x) ≥ 0 then ∑ ∫ f_n(x)dx = ∫ ∑ f_n(x)dx [24]. In our case, the expressions inside of the summation and the integral in Eq. (11) involve a probability density function and the sine of the scattering angle, both of which may take on values greater than or equal to zero, so reversal of integration and summation operations are valid. In the Composite Phase Function Scattering Method, each distribution’s phase function is calculated, and then the phase functions are weighted according to that distribution’s scattering coefficient and a composite phase function is then created. This was mentioned previously in Eq. (8). Next, a CDF is then constructed from the composite phase function. Substituting in Eq. (8), yields:

F (θ) = \int \bar{f} \sin θ d θ = \int \sum_{i = 1}^{m} \frac{μ_{s, i}}{{\bar{μ}}_{s}} f_{i} \sin θ d θ .

This demonstrates that the Composite Phase Function Scattering Method and the Distribution Selection Scattering Method both produce the same CDF since Eq. (11) and Eq. (12) are equal to one another.

Acknowledgments

The authors gratefully acknowledge support from the NIH-funded Laser Microbeam and Medical Program (P41 EB015890) and NIH CA71898.

References and links

1. J. R. Mourant, T. M. Johnson, S. Carpenter, A. Guerra, T. Aida, and J. P. Freyer, “Polarized angular dependent spectroscopy of epithelial cells and epithelial cell nuclei to determine the size scale of scattering structures,” J. Biomed. Opt. 7(3), 378–387 (2002). [CrossRef] [PubMed]

2. A. J. Gomes, V. Turzhitsky, S. Ruderman, and V. Backman, “Monte Carlo model of the penetration depth for polarization gating spectroscopy: Influence of illumination-collection geometry and sample optical properties,” Appl. Opt. 51(20), 4627–4637 (2012). [CrossRef] [PubMed]

3. L. Nieman, A. Myakov, J. Aaron, and K. Sokolov, “Optical sectioning using a fiber probe with an angled illumination-collection geometry: Evaluation in engineered tissue phantoms,” Appl. Opt. 43(6), 550–555 (2004). [CrossRef]

4. C. K. Hayakawa, J. Spanier, F. Bevilacqua, A. K. Dunn, J. S. You, B. J. Tromberg, and V. Venugopalan, “Perturbation Monte Carlo methods to solve inverse photon migration problems in heterogeneous tissues,” Opt. Lett. 26(17), 1335–1337 (2001). [CrossRef]

5. I. S. Seo, J. S. You, C. K. Hayakawa, and V. Venugopalan, “Perturbation and differential Monte Carlo methods for measurement of optical properties in a layered epithelial tissue model,” J. Biomed. Opt. 12(1), 014030 (2007). [CrossRef] [PubMed]

6. P. K. Yalavarthy, K. Karlekar, H. S. Patel, R. M. Vasu, M. Pramanik, P. C. Mathias, B. Jain, and P. K. Gupta, “Experimental investigation of perturbation Monte-Carlo based derivative estimation for imaging low-scattering tissue,” Opt. Express 13(3), 985–997 (2005). [CrossRef] [PubMed]

7. J. Heiskala, K. Kotilahti, L. Lipiinen, P. Hiltunen, P. E. Grant, and I. Nissila, “Optical tomographic imaging of activation of the infant auditory cortex using perturbation Monte Carlo with anatomical a priori information,” in Diffuse Optical Imaging of Tissue, Vol. 6629 of 2007 Proceedings of SPIE-OSA Biomedical Optics (Optical Society of America, 2007), paper 6629-29. [CrossRef]

8. J. Chen and X. Intes, “Time-gated perturbation Monte Carlo for whole body functional imaging in small animals,” Opt. Express 17(22), 19566–19579 (2009). [CrossRef] [PubMed]

9. J. Chen and X. Intes, “Comparison of Monte Carlo methods for fluorescence molecular tomography - computational efficiency,” Med. Phys. 38(10), 5788–5798 (2011). [CrossRef] [PubMed]

10. J. Nguyen, C. K Hayakawa, J. R. Mourant, and J. Spanier, “Perturbation Monte Carlo methods for tissue structure alterations,” Biomed. Opt. Express 4(10), 1946–1963 (2013). [CrossRef] [PubMed]

11. J. Ramachandran, T. M. Powers, S. Carpenter, A. Garcia-Lopez, J. P. Freyer, and J. R. Mourant, “Light scattering and microarchitectural differences between tumorigenic and non-tumorigenic cell models of tissue,” Opt. Express 15(7), 4039–4053 (2007). [CrossRef] [PubMed]

12. M. Bartlett, G. Huang, L. Larcom, and H. Jiang, “Measurement of particle size distribution in mammalian cells in vitro by use of polarized light spectroscopy,” Appl. Opt. 43(6), 1296–1307 (2004). [CrossRef] [PubMed]

13. J. M. Schmitt and G. Kumar, “Turbulent nature of refractive-index variations in biological tissue,” Opt. Lett. 21(16), 1310–1312 (1996). [CrossRef] [PubMed]

14. J. M. Schmitt and G. Kumar, “Optical scattering properties of soft tissue: a discrete particle model,” Appl. Opt. 37(13), 2788–2797 (1998). [CrossRef]

15. J. Spanier and E. M. Gelbard, Monte Carlo Principles and Neutron Transport Problems (Addison-Wesley, 1969).

16. C. K. Hayakawa, J. Spanier, and V. Venugopalan, “Comparative analysis of discrete and continuous absorption weighting estimators used in Monte Carlo simulations of radiative transport in turbid media,” J. Opt. Soc. Am. A 31(2), 301–311 (2014). [CrossRef]

17. J. R. Mourant, T. M. Johnson, and J. P. Freyer, “Characterizing mammalian cells and cell phantoms by polarized backscattering fiber-optic measurements,” Appl. Opt. 40(28), 5114–5123 (2001). [CrossRef]

18. J. Ramella-Roman, S. Prahl, and S. Jacques, “Three Monte Carlo programs of polarized light transport into scattering media: part I,” Opt. Express 13(12), 4420–4438 (2005). [CrossRef] [PubMed]

19. S. Bartel and A. H. Hielscher, “Monte Carlo simulations of the diffuse backscattering Mueller matrix for highly scattering media,” Appl. Opt. 39(10), 1580–1589 (2000). [CrossRef]

20. J. R. Mourant, T. M. Powers, T. J. Bocklage, H. M. Greene, M. H. Dorin, A. G. Waxman, M. M. Zsemlye, and H. O. Smith, “In vivo light scattering for the detection of cancerous and precancerous lesions of the cervix,” Appl. Opt. 48(10), D26–D35 (2009). [CrossRef] [PubMed]

21. T. M. Johnson and J. R. Mourant, “Polarized wavelength-dependent measurements of turbid media,” Opt. Express 4(6), 200–216 (1999). [CrossRef] [PubMed]

22. C. F. Bohren and D. R. Huffman, Absorption and Scattering of Light by Small Particles (Wiley-Interscience, 1983), Chap. 4.

23. X-5 Monte Carlo Team, “MCNP, a general Monte Carlo N-particle transport code, version 5, Report LA-UR-03-1987,” Technical Report LA-UR-03-1987 (Los Alamos National Laboratory, 2003).

24. H. Friedman, “A consistent Fubini-Tonelli theorem for nonmeasurable functions,” Illinois J. Math , 24(3), 390–395 (1980).

i	Average Size, ${\bar{r}}_{i} \pm σ (μ m)$	Number Density, N_s,i (cm⁻³)	n_medium	n_scatterer
1	0.03 ± 0.4	4 × 10¹³	1.33	1.46
2	0.45 ± 0.3	5 × 10¹⁰	1.35	1.40
3	4.8 ± 0.1	5 × 10⁷	1.37	1.49

Baseline Value, α	Perturbed Parameter, $\hat{α}$	g ranges	$\hat{{\bar{μ}}_{s}}$ ranges (cm⁻¹)
λ = 620 nm	500 – 720 nm	$\bar{g} = 0.952 - 0.953$	94.5 – 183
N_s,₁= 4.00 × 10¹³ cm⁻³	2.00 – 6.00 × 10¹³cm⁻³	g₁ = 0.183	123 – 127
N_s,₂= 5.00 × 10¹⁰ cm⁻³	2.5 – 7.5 × 10¹⁰cm⁻³	g₂ = 0.958	93.9 – 156
N_s,₃ = 5.00 × 10⁷ cm⁻³	2.5 – 7.5 × 10⁷cm⁻³	g₃ = 0.999	95.5 – 154
${\bar{r}}_{1} = 0.030 μ m$	0.015–0.045 µm	${\hat{g}}_{1} = 0.0502 - 0.345$	121 – 152
${\bar{r}}_{2} = 0.45 μ m$	0.38 – 0.52 µm	${\hat{g}}_{2} = 0.945 - 0.966$	94.2 – 174
${\bar{r}}_{3} = 4.8 μ m$	3.9 – 5.7 µm	${\hat{g}}_{3} = 0.998 - 0.999$	93.8 – 173

Parameter	Detector 1, ║	Detector 1, ⊥	Detector 2, ║	Detector 2, ⊥
λ	−12.9%, +16.1%	−6.4%, +12.9%	−6.4%, +12.9%	−6.4%, +6.4%
r₁	−50%, +40%	−30%, +30%	30%, +10%	−20%, +10%
r₂	−15.6%, +8.9%	−15.6%, +8.9%	−11.1%, +6.7%	−13.3%, 6.7%
r₃	−14.6%, +16.7%	−6.25%, +4.2%	−14.6%, +10.4%	−8.3%, +8.3%
N_s,₁	−50%, +50%	−50%, +50%	−50%, +50%	−50%, +50%
N_s,₂	−40%, +50%	−30%, +40%	−50%, +20%	−50%, +30%
N_s,₃	−50%, +50%	−50%, +30%	−40%, +40%	−10%, +20%
	Detector 1, ║	Detector, 1 ⊥	Detector 2, ║	Detector 2, ⊥
rel. err.	6.74 × 10⁻³	9.39 × 10⁻³	8.53 × 10⁻³	1.11 × 10⁻²

i	Average Size, ${\bar{r}}_{i} \pm σ (μ m)$	Number Density, N_s,i (cm⁻³)	n_medium	n_scatterer
1	0.03 ± 0.4	4 × 10¹³	1.33	1.46
2	0.45 ± 0.3	5 × 10¹⁰	1.35	1.40
3	4.8 ± 0.1	5 × 10⁷	1.37	1.49

Baseline Value, α	Perturbed Parameter, $\hat{α}$	g ranges	$\hat{{\bar{μ}}_{s}}$ ranges (cm⁻¹)
λ = 620 nm	500 – 720 nm	$\bar{g} = 0.952 - 0.953$	94.5 – 183
N_s,₁= 4.00 × 10¹³ cm⁻³	2.00 – 6.00 × 10¹³cm⁻³	g₁ = 0.183	123 – 127
N_s,₂= 5.00 × 10¹⁰ cm⁻³	2.5 – 7.5 × 10¹⁰cm⁻³	g₂ = 0.958	93.9 – 156
N_s,₃ = 5.00 × 10⁷ cm⁻³	2.5 – 7.5 × 10⁷cm⁻³	g₃ = 0.999	95.5 – 154
${\bar{r}}_{1} = 0.030 μ m$	0.015–0.045 µm	${\hat{g}}_{1} = 0.0502 - 0.345$	121 – 152
${\bar{r}}_{2} = 0.45 μ m$	0.38 – 0.52 µm	${\hat{g}}_{2} = 0.945 - 0.966$	94.2 – 174
${\bar{r}}_{3} = 4.8 μ m$	3.9 – 5.7 µm	${\hat{g}}_{3} = 0.998 - 0.999$	93.8 – 173

Development of perturbation Monte Carlo methods for polarized light transport in a discrete particle scattering model

Abstract

1. Introduction

2. Model problem

2.1. Light scattering model

2.2. Measurement geometry

3. Theory

3.1. Perturbation as change of measure

3.2. Perturbation Monte Carlo

4. Methods

5. Results

5.1. Scattering model parameter perturbations

5.2. PMC performance and conventional Monte Carlo relative error

5.3. PMC and cMC computational efficiency comparisons

6. Summary and conclusion

7. Appendix

7.1. Scattering algorithms

7.2. Proof of scattering algorithm equivalence

Acknowledgments

References and links

Cited By

Figures (7)

Tables (3)

Equations (12)

Biomedical Optics Express