1.

Introduction

Monte Carlo simulations have been broadly adopted by the biomedical optics community to simulate light propagation in scattering tissues on mesoscopic and macroscopic scales; i.e., on spatial scales comparable to and larger than the single scattering mean free path. Conventional Monte Carlo (cMC) simulations provide rigorous solutions to the radiative transport equation and can be configured for systems with complicated geometric and material features. However, cMC simulations can be computationally costly since this stochastic solver carries with it a solution uncertainty that scales as $\frac{1}{\sqrt{N}}$ where $N$ is the number of photons simulated.¹ Thus, any application that requires the execution of multiple cMC simulations; e.g., the resolution of an inverse problem, can easily have an associated computational cost that is impractical given the number of photons that must be simulated in order to obtain an estimate with sufficiently low uncertainty.

The computational cost associated with Monte Carlo simulation has motivated many groups to develop methods to improve its speed and efficiency for simulations of light transport in turbid media.² These methods can be generally categorized as follows: lookup table-based MC,³ scaled or “white” MC,^4–8 perturbation MC,⁹ parallel computing, cloud computing and/or graphic processor unit (GPU)-based methods^10–12 and variance reduction techniques.¹³ Amongst these, both lookup table and scaled methods suffer from restrictions to a fixed pre-defined geometry and prior binning of results that lead to reductions in accuracy.⁷ Parallel computing and GPU-based methods accelerate the speed of Monte Carlo simulations using innovations in compilers and hardware. However, the MC simulation engine remains unchanged, and often existing codes must be restructured to reap the benefits.

The perturbation Monte Carlo (pMC) method has been developed to rapidly obtain estimates for systems that are slightly modified, in terms of optical properties and/or geometry, relative to a reference cMC simulation. Moreover, the pMC framework facilitates implementation of a ‘sister’ method known as differential Monte Carlo (dMC) that enables the computation of sensitivities (Jacobian) for the resolution of inverse problems using gradient-based optimization methods.^14–16

The pMC method leverages correlated sampling by using a single set of random walks for simultaneous analysis of a ‘reference’ system together with any number of closely related systems which are modified in terms of optical properties and/or geometric characteristics.¹ pMC methods enable the rapid computation of RTE solutions for these closely related systems by post-processing the random walks from a database formed by characteristics of the reference MC simulation.^9,14,17–19 In this database, the weight, path length, and the number of collisions experienced by each detected photon are stored. pMC analysis modifies the weight of each tallied photon in the reference database based on its path length and number of collisions and change of the optical properties relative to the reference system.^1,9,20

Several studies have implemented pMC to improve computational efficiency and accuracy as compared to cMC simulations. Yamamoto and Sakamoto²¹ used pMC to reconstruct the optical characteristics of a heterogeneous, two-dimensional tissue model using temporal frequency domain data. This approach effectively reconstructed both scattering and absorption coefficients. For diffuse optical tomography applications, Yao et al.²² extended pMC to compute spatially and temporally resolved sensitivity profiles. A novel method for solving the inverse problem of quantitative photoacoustic tomography using pMC was provided in another paper by Leino et al.¹⁶ In this study, pMC was shown to be capable of estimating spatial distributions of both absorption and scattering parameters. These estimated distributions were quantitatively accurate over the full range of parameter values typical for biological tissues.

Despite these achievements, challenges remain to broadly apply pMC methods for the analysis of multispectral datasets. Prior studies have applied pMC to datasets acquired at a single or small number of wavelengths.^9,14,15,17 Given the broad range of optical properties spanned in multispectral datasets, reference MC simulations for multiple sets of optical properties are likely needed. Yet, the specific optical property values that should be chosen to minimize the number of reference simulations required are unknown. While it is known that pMC uncertainty, and therefore accuracy, degrades with increases in perturbation size (degree of the tissue optical property change),²⁰ a general method to quantify the growth in the pMC uncertainty with perturbation size has not been proposed. Currently, the accuracy of a pMC simulation can only be assessed after the pMC computation is performed. This is clearly undesirable since we would like to know a priori, how large a perturbation can be computed from a reference simulation before the accuracy of the resulting pMC estimate becomes unacceptable.

A priori quantitative prediction of the growth of pMC uncertainty with perturbation size using data from the reference simulation alone would facilitate the implementation of pMC/dMC methods. This is because once the reference simulation is performed the growth in pMC uncertainty with perturbation size could be quantified thereby identifying the range of optical properties for which the reference simulation could be utilized. This would facilitate the analysis of multispectral datasets where there are often large variations in optical absorption and scattering properties. Our first objective in this study is to identify the range of perturbation size that could be applied to a reference simulation while retaining the pMC estimate variance within a certain limit. To do so we determine the pMC estimate variance (which is directly related to the pMC estimate relative error) as a function of perturbation size and optical properties to identify an acceptable perturbation range for each reference simulation. We then develop a method for a priori prediction of pMC uncertainty using information from the reference simulation alone without explicitly computing the perturbed photon weights. This enables the prediction of the largest perturbation size for which pMC can still be used to provide an estimate with an acceptable relative error.

2.

Method

We start by describing the formulation of pMC and its use to obtain a mean detected photon weight and associated variance which captures the uncertainty of the mean estimate. Rigorous computation of the variance associated with a pMC estimate is best obtained by analyzing the population of perturbed photon weights. These photon weights are determined by post-processing the database that stores the characteristics of each detected photon from the reference simulation. Although the photons perturbed weight variance can be accurately calculated for each perturbation size, our objective is to avoid such rigorous calculations and determine an accurate variance estimate for a range of perturbation sizes using only data from the reference simulation. This includes various order moments of weight, the number of collisions, and path length distributions, along with their corresponding covariance values.

We consider a pMC simulation of a homogeneous semi-infinite tissue whereby the optical properties of the perturbed system are changed globally. In pMC, the weight of only those photons that are tallied at the detector under consideration is modified. The perturbed photon weight for the $i$’th photon ($W_{\mathrm{P},i}$) with perturbed optical properties (${\mu }_{a,\mathrm{P}},{\mu }_{s,\mathrm{P}}$) is computed from the $i$’th photon weight from the reference simulation ($W_{\mathrm{R},i}$), which utilized reference optical properties (${\mu }_{a,\mathrm{R}},{\mu }_{s,\mathrm{R}}$) as follows:⁹

The database obtained from the reference Monte Carlo simulation contains not only the weight of each photon but also the number of photon collisions and path length. Given that the photon weights, collisions, and path lengths in MC simulations of radiative transport are frequently not normally distributed,²³ we will examine the utility of applying a classical error propagation approach inclusive of the second-order (covariance) and third-order (skewness) terms to estimate the variance associated with the mean perturbed photon weight. Application of this approach to Eq. (1) provides the following equation to estimate the variance of the population of perturbed photon weights ${\sigma }_{W_{\mathrm{P}}}^2$:²⁴

In this equation ${\sigma }_x$ and ${\sigma }_x^2$ represent the standard deviation and variance of the random variable $x$, respectively, shown in Eq. (4a). ${\sigma }_{x,y}$ represents the covariance between the two random variables $x$ and $y$ shown in Eq. (4b). ${\gamma }_x$ provides a measure of skewness of the distribution of the random variable $x$,²⁵ as defined by Eq. (4c). The product ${\gamma }_x{\sigma }_x^3$ represents the third moment of the random variable $x$:

In general, for a specific detector collecting photons over a finite interval of space, time, and/or propagation direction, only a subset of the $N$ photons that are simulated are tallied at the detector ($N_{\mathrm{T}}$) while the remaining photons ($N_{\mathrm{U}}$) go untallied. We thus recast Eq. (2) defining the variance of the perturbed photon weights $W_{\mathrm{P}}$, in a form that separates contributions of the tallied and untallied photons to the variance of the perturbed photon weight. As detailed in Section A of the Supplemental Materials, this reformulation results in the following expression for the variance of the perturbed estimate:

The nonlinearities inherent in Eq. (1) can lead to a large dynamic range of the perturbed photon weights. For this reason, we choose to apply Eq. (3) on the linearized form of Eq. (1) and estimate the variance of the natural logarithm of the photon weights, ${\sigma }_{\ln (W_{\mathrm{P},\mathrm{T}})}^2$ as follows:

Once ${\sigma }_{\ln (W_{\mathrm{P},\mathrm{T}})}^2$ has been calculated, we estimate ${\sigma }_{W_{\mathrm{P},\mathrm{T}}}^2$ using:

Finally, throughout the estimation process, we replace ${\overline{W}}_{\mathrm{P}}$ with ${\overline{W}}_{\mathrm{R}}$ to eliminate the need to calculate the mean perturbed weight. This replacement is expected to provide a reasonable approximation for small perturbations.

Figure 1 summarizes the conventional process to calculate the pMC variance along with our proposed process. The conventional process requires performing the pMC computation (steps 2 to 4) for each perturbation size of interest. By contrast, our proposed strategy computes statistical metrics that characterize the distributions of $W_{\mathrm{R}}$, $j$, and $L$ from the reference simulation alone (step 2) from which the pMC variance can be estimated (steps 3 and 4) for any perturbation size of interest. The equations characterizing the variation of pMC variance with perturbation size and the accuracy of our pMC variance estimate are presented in Eqs. (9) and (10).

Fig. 1

Schematic representation of the conventional variance calculation and our proposed variance estimation methods.

3.

Model Problem

To examine the performance of our method, we considered the case of spatially resolved reflectance in a homogeneous, semi-infinite medium. We performed conventional Monte Carlo simulations in which 20 million photons are launched from a directional point source in the reference simulations which utilize continuous absorption weighting (CAW).²³ We consider a medium with fixed refractive index ($n=1.4$), single-scattering anisotropy ($g=0.8$), and transport mean free path ($l^{*}=1\mathrm{mm}$). We consider five different media with different ratios of reduced scattering to absorption coefficients $({\mu }_s^{\prime }/{\mu }_a)$. Collection of photons happens at a proximal detector positioned immediately adjacent to the source spanning radial locations $\rho \in [0\text{to}0.2]\mathrm{mm}$ and a distal detector spanning radial locations $\rho \in [3\text{to}3.2]\mathrm{mm}$. These two detectors are chosen to examine the characteristics of the pMC estimates as the collected signals at these two detectors have differing sensitivities to changes in optical absorption and scattering.¹⁴ As described below, we also performed cMC simulations utilizing the Russian Roulette method²⁶ with a weight threshold of ${10}^{-3}$ and a survival probability of 0.1. Table 1 provides the optical properties for the five reference cases investigated.

Table 1

Optical parameters used for reference MC simulations with detectors placed at ρ∈[0 to 0.2] mm (proximal) and ρ∈[3 to 3.2]  mm (distal), g=0.8 was used in all simulations.

The database resulting from each reference simulation was processed using pMC for scattering perturbations (${ϵ}_s$) over the range of $[-20\%,20\%]$ in 5% increments and all other optical properties were left unchanged. We restricted our analysis to the consideration of scattering perturbations, since for conventional MC simulations utilizing CAW, perturbations in absorption can be accommodated without the loss of accuracy regardless of perturbation size.⁷

We introduce two metrics to characterize pMC performance and our ability to accurately calculate the pMC relative error using information from the reference MC simulation alone. First, we define a metric ${\overline{\mathrm{\Delta }}}_{\mathrm{P}}$ called the “degradation degree” that quantifies the relative error of a pMC estimate relative to the reference simulation. Computation of the variation of ${\overline{\mathrm{\Delta }}}_{\mathrm{P}}$ with perturbation size (${ϵ}_s$) enables the examination of the intrinsic accuracy of pMC. We also define $\delta$ which quantifies the relative difference between our estimate for the pMC relative error using information from the reference MC simulation alone compared with the actual pMC relative error. These two metrics are defined as follows:

4.

Results and Discussions

Figure 2 shows the variation in the degradation degree ${\overline{\mathrm{\Delta }}}_{\mathrm{P}}$ for the five sets of optical properties based on rigorous application of pMC calculations for ${ϵ}_s$ range of $[-20\%,20\%]$. The mean and variance of all the pMC calculations are provided in Section B of the Supplemental Material. For the proximal detector, we see minimal changes in the pMC relative error as compared with the reference simulation, with only slight increases observed at scattering perturbations of $\pm 20\%$. We do see, however, a slight worsening of accuracy for higher values of $({\mu }_s^{\prime }/{\mu }_a)$. The situation is notably different for the distal detector with much more significant increases in relative error as compared with the reference simulations for both positive and negative scattering perturbations. Importantly, we see sharper escalations in pMC relative error for the simulations utilizing larger values of $({\mu }_s^{\prime }/{\mu }_a)$. Also, in the ${ϵ}_s$ range of $[-10\%,10\%]$, we generally observe that for the reference simulations performed using larger $({\mu }_s^{\prime }/{\mu }_a)$ values, the pMC relative error increases more sharply for positive scattering perturbations as compared with negative perturbations. Whereas for reference simulations performed at lower $({\mu }_s^{\prime }/{\mu }_a)$ values, the pMC relative error degrades more sharply for negative perturbations.

Fig. 2

Variation of the degradation degree, ${\overline{\mathrm{\Delta }}}_{\mathrm{P}}$, with scattering perturbation size ${ϵ}_s$ and optical properties $({\mu }_s^{\prime }/{\mu }_a)$.

These results indicate that the optical properties do not play as prominent role in pMC accuracy for the proximal detector as compared with the distal detector, where reference simulations performed at lower values of $({\mu }_s^{\prime }/{\mu }_a)$ result in pMC predictions with much higher fidelity. Also, since the changes of ${\overline{\mathrm{\Delta }}}_{\mathrm{P}}$ relative to ${ϵ}_s$ are asymmetric, it appears that positive scattering perturbations may be more robust in terms of reduced degradation in the pMC relative error for cases with lower $({\mu }_s^{\prime }/{\mu }_a)$ whereas for higher $({\mu }_s^{\prime }/{\mu }_a)$ increases in the degradation degree are smaller for negative scattering perturbations.

To evaluate the accuracy in estimating the pMC relative error from the reference simulation as compared with a rigorous analysis obtained through conventional process of variance calculation, in Fig. 3 we display values for the $\delta$ metric for all the cases examined.

Fig. 3

Variation of $\delta$ with scattering perturbation size ${ϵ}_s$ and optical properties $({\mu }_s^{\prime }/{\mu }_a)$.

Our estimation of the pMC relative error based on distribution metrics of the reference simulation data alone incurs a relative error that generally grows with perturbation size for all reference cases. Our estimation method shows high fidelity in predicting the pMC relative error for the proximal detector over a substantial range of ${ϵ}_s$ values. On the other hand, for the distal detector, the range of ${ϵ}_s$ values that yield accurate estimates are confined to a much narrower range of ${ϵ}_s$. Moreover, the growth in relative error is asymmetrical and our method provides better accuracy for positive scattering perturbations. Figure 3 also shows that our pMC error estimates are conservative and that we almost invariably provide overestimates of the relative error. In the case of the proximal detector, our method provides pMC relative error estimates within 10% of the true value for scattering perturbations in the range of ${ϵ}_s=\pm 10\%$. The accuracy of our method for the distal detector is notably worse with error estimates within 15% for scattering perturbations in the range of ${ϵ}_s=[-4\%,4\%]$ with much poorer performance outside this range.

To gain insight into the performance characteristics of our method to estimate pMC relative error, we examined the distributions of the photon weight $W_{\mathrm{R}}$, number of collisions $j$, and photon path length $L$ for both proximal and distal detector locations. The distributions for the reference cases with $({\mu }_s^{\prime }/{\mu }_a)$ of 5 and 100 are shown in Figs. 4 and 5, respectively. Distributions for the other three cases can be found in Section C of the Supplemental Material.

Fig. 4

Histograms for the detected photons tallying $W_{\mathrm{R}}$, $j$, and $L$ for the reference simulation performed at $({\mu }_s^{\prime }/{\mu }_a)=5$. $N_{\mathrm{T}}$ = 680,088 and 411,540 for the (a) proximal and (b) distal detector, respectively.

Fig. 5

Histograms for the detected photons tallying $W_{\mathrm{R}}$, $j$, and $L$ for the reference simulation performed at $({\mu }_s^{\prime }/{\mu }_a)=100$. ${\mathit{N}}_{\mathrm{T}}=\mathrm{817,515}$ and 405,562 for the (a) proximal and (b) distal detector, respectively.

For both cases shown in Figs. 4 and 5, we observe strongly skewed distributions with long, sparsely populated tails. This feature is stronger in case of the distal detector for all reference simulations. Given that the original error propagation formula is based on a Taylor series expansion and assumes that the random variables are normally distributed, the performance of the estimation provided by Eq. (3) degrades when applied to random variables that deviate from these conditions.

To improve the accuracy of our method, we attempt to reduce the magnitude of the moments of the photon weight, path length, and collision distributions by utilizing RR²⁶ in our reference simulations. RR is an unbiased method for terminating simulated photons once their weight falls below a predefined threshold with a fair game probability. Once a photon’s weight drops below a specified weight threshold during its propagation, at the next collision its weight is either amplified by a factor of $1/p$ with a survival probability $p$ or terminated with probability ($1-p$).

Using this general approach, we selected a weight threshold of ${10}^{-3}$ and a survival probability of 0.1. Figures 6 and 7 show the distribution of the same random variables shown in Figs. 4 and 5 after implementation of RR for the reference simulations with $({\mu }_s^{\prime }/{\mu }_a)$ of 5 and 100, respectively. Distributions for the other three cases can be found in Section D of the Supplemental Material. Appearance of discontinuities in the $j$ and $L$ histograms is reflective of the tallied photon subpopulation that has undergone the RR. Comparison of the histograms utilizing the RR method with their counterparts in Figs. 4 and 5 illustrates a reduced range of values for each random variable and more compact and densely populated distributions.

Fig. 6

Histograms for the detected photons tallying $W_{\mathrm{R}}$, $j$, and $L$ for the reference simulation performed at $({\mu }_s^{\prime }/{\mu }_a)=5$ using RR. ${\mathit{N}}_{\mathrm{T}}=\mathrm{681,805}$ and 382,485 for the (a) proximal and (b) distal detector, respectively. $L_{\mathrm{thr}}$ indicates the path length corresponding to the RR weight threshold.

Fig. 7

Histograms for the detected photons tallying $W_{\mathrm{R}}$, $j$, and $L$ for the reference simulation performed at $({\mu }_s^{\prime }/{\mu }_a)=100$ using RR. ${\mathit{N}}_{\mathrm{T}}=\mathrm{816,995}$ and 404,955 for the (a) proximal and (b) distal detector, respectively. $L_{\mathrm{thr}}$ indicates the path length corresponding to the RR weight threshold.

We computed the various moments for the reference simulations both before and after implementation of RR and analyzed them after normalization to eliminate any dependence on the actual values. Tables 2 and 3 show the normalized characteristics of the distribution of random variables corresponding to the reference MC simulations for $({\mu }_s^{\prime }/{\mu }_a)$ of 5 and 100, respectively, before and after applying RR.

Table 2

Normalized variance (σx2/x¯2), covariance (σx,y/σxσy), and skewness γx metrics of the logarithm of the photon weight [log (WR)], number of photon collisions (j), and photon path length (L) for reference MC simulations performed at (μs′/μa)=5 before and after applying the RR method.

Table 3

Normalized variance (σx2/x¯2), covariance (σx,y/σxσy) and skewness γx metrics of the logarithm of the photon weight [log (WR)], number of photon collisions (j), and photon path length (L) for reference MC simulations performed at (μs′/μa)=100 before and after applying the RR method.

Based on these results, the normalized variance for all random variables decreases after applying RR technique. Similarly, a reduction in normalized covariance (correlation) is observed. This makes sense since the photon reweighting accomplished by RR weakens the strict dependence between the photon weight and both path length and number of collisions. The skewness of the random variable distributions also reduced dramatically after using RR. It is also clear that the reductions are less dramatic for the reference with $({\mu }_s^{\prime }/{\mu }_a)=100$ compared with $({\mu }_s^{\prime }/{\mu }_a)=5$. This is because for $({\mu }_s^{\prime }/{\mu }_a)=100$, fewer photons undergo RR because of the larger path length is necessary for the threshold photon weight to be reached.

Using the photon databases generated from these new reference simulations that utilized RR, we again performed a rigorous calculation of the pMC relative error as well as estimated the relative error utilizing only the distribution characteristics of the reference simulation data and our error propagation method. Figure 8 shows the degradation degree metric indicating how the pMC relative error varies with perturbation size using the RR reference simulations. These results are comparable to those shown in Fig. 2, which reports the same metric for the reference simulations that were performed without use of RR.

Fig. 8

Variation of the degradation degree, ${\overline{\mathrm{\Delta }}}_{\mathrm{P}}$, with scattering perturbation size ${ϵ}_s$ and optical properties $({\mu }_s^{\prime }/{\mu }_a)$ using RR in the reference simulations.

Figure 9 shows the relative difference $\delta$ between our estimate and the actual pMC relative error.

Fig. 9

Variation of $\delta$ with scattering perturbation size ${ϵ}_s$ and optical properties $({\mu }_s^{\prime }/{\mu }_a)$ using RR in the reference simulations.

The utilization of RR clearly improves our ability to estimate pMC relative error using data from the reference simulation alone. For the proximal detector, the use of RR improves our pMC relative error estimates most notably for scattering perturbations beyond $\pm 8\%$. With the usage of RR, we can estimate the pMC relative error within 5% of the true value for scattering perturbations in the range of $[-15\%,+20\%]$. Our predictions for the distal detector are also notably improved and the usage of RR provides relative error estimates within 20% for scattering perturbations in the range of $[-8\%,+15\%]$. For the distal detector, reference simulation with $({\mu }_s^{\prime }/{\mu }_a)=5$ seems to provide the largest range of perturbation for which our approximation method provides the highest accuracy. While the accuracy of our method is very strong overall for the proximal detector, we continue to observe poorer performance for higher values of $({\mu }_s^{\prime }/{\mu }_a)$, which is also the case for the distal detector.

To explain the underlying reasons for this, we should first note that the use of reference simulations utilizing RR improved the overall pMC error predictions more so for the distal detector as compared with the proximal detector. Improvement in the prediction accuracy was observed to be greater for lower values of $({\mu }_s^{\prime }/{\mu }_a)$. These characteristics are expected as photons collected at the distal detector typically have a larger path length as compared with those collected at the proximal detector. Moreover, photons propagating in the more highly absorbing medium, $({\mu }_s^{\prime }/{\mu }_a)=5$, need only travel a path length of 41 mm before RR is invoked as opposed to 698 mm for the highly scattering medium of $({\mu }_s^{\prime }/{\mu }_a)=100$. As a result, for the $({\mu }_s^{\prime }/{\mu }_a)=5$ medium, 470 and 13,550 photons underwent RR reweighting for the proximal and distal detector, respectively. By contrast for the $({\mu }_s^{\prime }/{\mu }_a)=100$ medium, only 40 and 140 of the detected photons underwent RR reweighting for the proximal and distal detector, respectively.

Taken collectively, our results suggest that reference simulations should be run using the lowest possible value of $({\mu }_s^{\prime }/{\mu }_a)$ since the intrinsic growth of the pMC variance, as reported by the degradation degree metric ${\overline{\mathrm{\Delta }}}_{\mathrm{P}}$, is minimized. This result is consistent with the theoretical analysis of Ref. 20 that shows the perturbation range over which pMC estimates can be obtained, grows as the probability of scattering is decreased. Moreover, it is also reassuring to observe that our ability to estimate the pMC relative error using reference simulation data alone performs best for lower $({\mu }_s^{\prime }/{\mu }_a)$ values. The results also show that the use of RR reduces the intrinsic variance of the reference simulations while also improving our ability to accurately estimate the pMC variance. Moreover, the threshold weight should be chosen to adequately limit excessively long path lengths that result in distributions of photon weight, path length, and collisions with extended, sparsely populated tails.

5.

Conclusions

In conclusion, we have presented a method to estimate the relative error associated with the use of perturbation Monte Carlo estimates using distribution metrics of the reference simulation data alone. This ability reduces pMC computational cost and provides specific guidance for the selection of optical properties for the placement of reference simulations. Moreover, we have shown that the use of RR is advantageous in reducing the intrinsic relative error characteristics of reference simulations used for deriving pMC estimates as well as providing a large improvement in the perturbation range over which we can predict the relative error. Our results show conclusively that the range of scattering perturbation while minimizing the growth in the relative error of the resulting pMC estimates is best accomplished when $({\mu }_s^{\prime }/{\mu }_a)$ is low. This result is consistent with the analysis of Ref. 20 who showed that the allowable perturbation range of pMC grows as the probability of scattering is decreased. Our results suggest that to utilize pMC for predictions over a wide range of optical properties, reference simulations should utilize CAW with optical properties corresponding to low $({\mu }_s^{\prime }/{\mu }_a)$ values over the desired range of ${\mu }_s$ values. This is because absorption perturbations can be computed exactly when using CAW and can then be employed to compute pMC estimates for cases where $({\mu }_s^{\prime }/{\mu }_a)$ is large. Finally, we note that the results provided in this paper represent a “worst case” for the application of pMC in that we perturbed the properties of the entire medium. However, in most applications, the perturbation will be applied to only a subdomain of the entire volume being considered. In these cases, we expect that our methodology will provide accurate results for a larger range of scattering perturbations.