Analysis of relative error in perturbation Monte Carlo simulations of radiative transport

Abstract Significance Perturbation and differential Monte Carlo (pMC/dMC) methods, used in conjunction with nonlinear optimization methods, have been successfully applied to solve inverse problems in diffuse optics. Application of pMC to systems over a large range of optical properties requires optimal “placement” of baseline conventional Monte Carlo (cMC) simulations to minimize the pMC variance. The inability to predict the growth in pMC solution uncertainty with perturbation size limits the application of pMC, especially for multispectral datasets where the variation of optical properties can be substantial. Aim We aim to predict the variation of pMC variance with perturbation size without explicit computation of perturbed photon weights. Our proposed method can be used to determine the range of optical properties over which pMC predictions provide sufficient accuracy. This method can be used to specify the optical properties for the reference cMC simulations that pMC utilizes to provide accurate predictions over a desired optical property range. Approach We utilize a conventional error propagation methodology to calculate changes in pMC relative error for Monte Carlo simulations. We demonstrate this methodology for spatially resolved diffuse reflectance measurements with ±20% scattering perturbations. We examine the performance of our method for reference simulations spanning a broad range of optical properties relevant for diffuse optical imaging of biological tissues. Our predictions are computed using the variance, covariance, and skewness of the photon weight, path length, and collision distributions generated by the reference simulation. Results We find that our methodology performs best when used in conjunction with reference cMC simulations that utilize Russian Roulette (RR) method. Specifically, we demonstrate that for a proximal detector placed immediately adjacent to the source, we can estimate the pMC relative error within 5% of the true value for scattering perturbations in the range of [−15%,+20%]. For a distal detector placed at ∼3 transport mean free paths relative to the source, our method provides relative error estimates within 20% for scattering perturbations in the range of [−8%,+15%]. Moreover, reference simulations performed at lower (μs′/μa) values showed better performance for both proximal and distal detectors. Conclusions These findings indicate that reference simulations utilizing continuous absorption weighting (CAW) with the Russian Roulette method and executed using optical properties with a low (μs′/μa) ratio spanning the desired range of μs values, are highly advantageous for the deployment of pMC to obtain radiative transport estimates over a wide range of optical properties.


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 1 ffiffiffi N p where N is the number of photons simulated. 1 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. 2 These methods can be generally categorized as follows: lookup table-based MC, 3 scaled or "white" MC, 4-8 perturbation MC, 9 parallel computing, cloud computing and/or graphic processor unit (GPU)-based methods [10][11][12] and variance reduction techniques. 13 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. 7 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][15][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. 1 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][18][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 21 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. 22 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. 16 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), 20 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.

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 postprocessing 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 P;i ) with perturbed optical properties (μ a;P ; μ s;P ) is computed from the i'th photon weight from the reference simulation (W R;i ), which utilized reference optical properties (μ a;R ; μ s;R ) as follows: 9 E Q -T A R G E T ; t e m p : i n t r a l i n k -; e 0 0 1 ; 1 1 7 ; 1 8 6 where μ s and μ t are the scattering and total interaction coefficients (μ t ¼ μ a þ μ s , μ a being the absorption coefficient), respectively, j i is the number of collisions the i'th tallied photon experiences in the medium, and L i is the path length taken by the i'th photon prior to detection. The subscripts 'P' and 'R' refer to the perturbed and reference cases, respectively. To generate a set of perturbed photon weights, Eq. (1) is applied to each photon that is tallied at the detector in the reference simulation. The variance of the entire population of photon weights σ 2 W P , for the perturbed case is given as where N is the total number of photons launched and W P is the mean weight of the population of photon weights for the perturbed case, W P;i . Our goal is to estimate the variance of the perturbed photon weights from the reference simulation alone, i.e., without explicitly calculating the perturbed photon weights needed to apply Eq. (2). 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, 23 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 σ 2 W P : 24 E Q -T A R G E T ; t e m p : i n t r a l i n k -; e 0 0 3 ; 1 1 4 ; 5 5 2 In this equation σ x and σ 2 x represent the standard deviation and variance of the random variable x, respectively, shown in Eq. (4a). σ x;y represents the covariance between the two random variables x and y shown in Eq. (4b). γ x provides a measure of skewness of the distribution of the random variable x, 25 as defined by Eq. (4c). The product γ x σ 3 x represents the third moment of the random variable x: E Q -T A R G E T ; t e m p : i n t r a l i n k -; e 0 0 4 a ; 1 1 4 ; 3 9 4 E Q -T A R G E T ; t e m p : i n t r a l i n k -; e 0 0 4 c ; 1 1 4 ; 3 1 6 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 T ) while the remaining photons (N U ) go untallied. We thus recast Eq. (2) defining the variance of the perturbed photon weights W 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: E Q -T A R G E T ; t e m p : i n t r a l i n k -; e 0 0 5 ; 1 1 4 ; 2 1 2 where σ 2 W P;T is the variance of the sub-population of photon weights that are tallied at the detector. Φ ¼ W P − W P;T represents the difference of the mean weight over the entire population of simulated photons (N) and only the population of photons that are tallied at the detector (N T ). In this restructured equation, the contribution of the untallied photon population to the overall variance is accounted for via the N U N W 2 P term and the impact of the tallied photons is expressed through the variance of the population of the tallied photon weights alone plus the correction factor N T N Φ 2 that accounts for the differing sizes of the tallied population and the entire population of simulated photons.
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, σ 2 lnðW P;T Þ as follows: E Q -T A R G E T ; t e m p : i n t r a l i n k -; e 0 0 6 ; 1 1 7 ; 5 0 4 Once σ 2 lnðW P;T Þ has been calculated, we estimate σ 2 W P;T using: ; t e m p : i n t r a l i n k -; e 0 0 7 ; 1 1 7 ; 4 5 0 E Q -T A R G E T ; t e m p : i n t r a l i n k -; e 0 0 8 ; 1 1 7 ; 4 0 1 Finally, throughout the estimation process, we replace W P with W 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 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).

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). 23 We consider a medium with fixed refractive index (n ¼ 1.4), single-scattering anisotropy (g ¼ 0.8), and transport mean free path (l Ã ¼ 1 mm). We consider five different media with different ratios of reduced scattering to absorption coefficients ðμ 0 s ∕μ a Þ. Collection of photons happens at a proximal detector positioned immediately adjacent to the source spanning radial locations ρ ∈ ½0 to 0.2 mm and a distal detector spanning radial locations ρ ∈ ½3 to 3.2 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. 14 As described below, we also performed cMC simulations utilizing the Russian Roulette method 26 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. 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. 7 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 Δ P called the "degradation degree" that quantifies the relative error of a pMC estimate relative to the reference simulation. Computation of the variation of Δ P with perturbation size (ϵ s ) enables the examination of the intrinsic accuracy of pMC. We also define δ 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: E Q -T A R G E T ; t e m p : i n t r a l i n k -; e 0 0 9 ; 1 1 4 ; 4 0 8 E Q -T A R G E T ; t e m p : i n t r a l i n k -; e 0 1 0 ; 1 1 4 ; 3 5 6 where W R and W P are the mean photon weight from the reference Monte Carlo simulation and the pMC simulation, respectively. We do see, however, a slight worsening of accuracy for higher values of ðμ 0 s ∕μ 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 ðμ 0 s ∕μ a Þ. Also, in the ϵ s range of ½−10%; 10%, we generally observe that for the reference simulations performed using larger ðμ 0 s ∕μ 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 ðμ 0 s ∕μ a Þ values, the pMC relative error degrades more sharply for negative perturbations.

Results and Discussions
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 ðμ 0 s ∕μ a Þ result in pMC predictions with much higher fidelity. Also, since the changes of Δ 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 ðμ 0 s ∕μ a Þ whereas for higher ðμ 0 s ∕μ 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 δ metric for all the cases examined.
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 ¼ AE10%. 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 R , number of collisions j, and photon  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 26 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 ðμ 0 s ∕μ a Þ of 5 and   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 ðμ 0 s ∕μ a Þ of 5 and 100, respectively, before and after applying RR.
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 ðμ 0 s ∕μ a Þ ¼ 100 compared with ðμ 0 s ∕μ a Þ ¼ 5. This is because for ðμ 0 s ∕μ 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. Figure 9 shows the relative difference δ between our estimate and the actual pMC relative error.
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 Table 2 Normalized variance ðσ 2 x ∕x 2 Þ, covariance ðσ x;y ∕σ x σ y Þ, and skewness γ x metrics of the logarithm of the photon weight [log (W R )], number of photon collisions (j), and photon path length (L) for reference MC simulations performed at ðμ 0 s ∕μ a Þ ¼ 5 before and after applying the RR method.
980 relative error estimates most notably for scattering perturbations beyond AE8%. 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 ðμ 0 s ∕μ 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 ðμ 0 s ∕μ 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 Table 3 Normalized variance ðσ 2 x ∕x 2 Þ, covariance ðσ x;y ∕σ x σ y Þ and skewness γ x metrics of the logarithm of the photon weight [log (W R )], number of photon collisions (j), and photon path length (L) for reference MC simulations performed at ðμ 0 s ∕μ a Þ ¼ 100 before and after applying the RR method.   Fig. 8 Variation of the degradation degree, Δ P , with scattering perturbation size ϵ s and optical properties ðμ 0 s ∕μ a Þ using RR in the reference simulations.
detector as compared with the proximal detector. Improvement in the prediction accuracy was observed to be greater for lower values of ðμ 0 s ∕μ 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, ðμ 0 s ∕μ 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 ðμ 0 s ∕μ a Þ ¼ 100. As a result, for the ðμ 0 s ∕μ a Þ ¼ 5 medium, 470 and 13,550 photons underwent RR reweighting for the proximal and distal detector, respectively. By contrast for the ðμ 0 s ∕μ 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 ðμ 0 s ∕μ a Þ since the intrinsic growth of the pMC variance, as reported by the degradation degree metric Δ 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 ðμ 0 s ∕μ 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.

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 ðμ 0 s ∕μ 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 ðμ 0 s ∕μ a Þ values over the desired range of μ 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 ðμ 0 s ∕μ 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.

Disclosures
The authors have no conflicts of interest.

Code, Data, and Materials Availability
The MCCL open-source Monte Carlo computation engine was used to generate the results of this study and can be accessed here: https://virtualphotonics.org/software-mccl.