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

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.

MCCL: an open-source software application for Monte Carlo simulations of radiative transport.

The Monte Carlo Command Line application (MCCL) is an open-source software package that provides Monte Carlo simulations of radiative transport through heterogeneous turbid media. MCCL is available on GitHub through our virtualphotonics.org website, is actively supported, and carries extensive documentation. Here, we describe the main technical capabilities, the overall software architecture, and the operational details of MCCL.

Comparative analysis of discrete and continuous absorption weighting estimators used in Monte Carlo simulations of radiative transport in turbid media: erratum.

This erratum corrects the relative error plots and references in our paper [J. Opt. Soc. Am. A31, 301 (2014)JOAOD60740-323210.1364/JOSAA.31.000301].

Optical sampling depth in the spatial frequency domain.

We present a Monte Carlo (MC) method to determine depth-dependent probability distributions of photon visitation and detection for optical reflectance measurements performed in the spatial frequency domain (SFD). These distributions are formed using an MC simulation for radiative transport that utilizes a photon packet weighting procedure consistent with the two-dimensional spatial Fourier transform of the radiative transport equation. This method enables the development of quantitative metrics for SFD optical sampling depth in layered tissue and its dependence on both tissue optical properties and spatial frequency. We validate the computed depth-dependent probability distributions using SFD measurements in a layered phantom system with a highly scattering top layer of variable thickness supported by a highly absorbing base layer. We utilize our method to establish the spatial frequency-dependent optical sampling depth for a number of tissue types and also provide a general tool to determine such depths for tissues of arbitrary optical properties.

Momentum transfer Monte Carlo for the simulation of laser speckle imaging and its application in the skin.

Due to its simplicity and low cost, laser speckle imaging (LSI) has achieved widespread use in biomedical applications. However, interpretation of the blood-flow maps remains ambiguous, as LSI enables only limited visualization of vasculature below scattering layers such as the epidermis and skull. Here, we describe a computational model that enables flexible in-silico study of the impact of these factors on LSI measurements. The model uses Monte Carlo methods to simulate light and momentum transport in a heterogeneous tissue geometry. The virtual detectors of the model track several important characteristics of light. This model enables study of LSI aspects that may be difficult or unwieldy to address in an experimental setting, and enables detailed study of the fundamental origins of speckle contrast modulation in tissue-specific geometries. We applied the model to an in-depth exploration of the spectral dependence of speckle contrast signal in the skin, the effects of epidermal melanin content on LSI, and the depth-dependent origins of our signal. We found that LSI of transmitted light allows for a more homogeneous integration of the signal from the entire bulk of the tissue, whereas epi-illumination measurements of contrast are limited to a fraction of the light penetration depth. We quantified the spectral depth dependence of our contrast signal in the skin, and did not observe a statistically significant effect of epidermal melanin on speckle contrast. Finally, we corroborated these simulated results with experimental LSI measurements of flow beneath a thin absorbing layer. The results of this study suggest the use of LSI in the clinic to monitor perfusion in patients with different skin types, or inhomogeneous epidermal melanin distributions.

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

We present a polarization-sensitive, transport-rigorous perturbation Monte Carlo (pMC) method to model the impact of optical property changes on reflectance measurements within a discrete particle scattering model. The model consists of three log-normally distributed populations of Mie scatterers that approximate biologically relevant cervical tissue properties. Our method provides reflectance estimates for perturbations across wavelength and/or scattering model parameters. We test our pMC model performance by perturbing across number densities and mean particle radii, and compare pMC reflectance estimates with those obtained from conventional Monte Carlo simulations. These tests allow us to explore different factors that control pMC performance and to evaluate the gains in computational efficiency that our pMC method provides.

Targeted narrowband intense pulsed light on cutaneous vasculature.

BACKGROUND AND OBJECTIVES: Laser based therapies are the standard treatment protocol for port wine stain in the United States, but complete removal is infrequently achieved. Intense pulsed light (IPL) offers a broadband light spectrum approach as a viable treatment alternative. Previous studies suggest that IPL can be more effective in treatment of port wine stain by utilizing multiple wavelengths to selectively target different peaks in oxy- and deoxy-hemoglobin. Our study objectives were to (i) determine a characteristic radiant exposure able to achieve persistent vascular shutdown with narrowband IPL irradiation, (ii) determine the degree to which narrowband IPL irradiation can achieve persistent vascular shutdown, and (iii) compare the effectiveness of narrowband IPL radiation to single wavelength pulsed dye laser (PDL) irradiation in achieving persistent vascular shutdown. STUDY DESIGN/MATERIALS AND METHODS: We utlized either single pulse or double, stacked pulses in narrowband IPL experiments, with the IPL operating over a 500-600 nm wavelength range on the rodent dorsal window chamber model. We compared the results from our narrowband IPL experiments to acquired PDL data from a previous study and determined that narrowband IPL treatments can also produce persistent vascular shutdown. We ran Monte Carlo simulations to investigate the relationship between absorbed energy, wavelength, and penetration depth. RESULTS: For single and double pulse narrowband IPL irradiation we observed (i) little to no change in blood flow, resulting in no persistent vascular shutdown, (ii) marked acute disruption in blood flow and vascular structure, followed by partial to full recovery of blood flow, also resulting in no persistent vascular shutdown, and (iii) immediate changes in blood flow and vascular structure, resulting in prolonged and complete vascular shutdown. Monte Carlo modeling resulted in a 53.2% and 69.0% higher absorbed energy distribution in the top half and the total simulated vessel when comparing the composite narrowband IPL to the 595 nm (PDL), respectively. CONCLUSIONS: Our data collectively demonstrate the potential to achieve removal of vascular lesions using a 500-600 nm range. Additionally, the narrowband IPL was tuned to optimize a specific wavelength range that can be used to treat PWS, whereas the PDL can only operate at one discrete wavelength.

Rapid computation of the amplitude and phase of tightly focused optical fields distorted by scattering particles.

We develop an efficient method for accurately calculating the electric field of tightly focused laser beams in the presence of specific configurations of microscopic scatterers. This Huygens-Fresnel wave-based electric field superposition (HF-WEFS) method computes the amplitude and phase of the scattered electric field in excellent agreement with finite difference time-domain (FDTD) solutions of Maxwell's equations. Our HF-WEFS implementation is 2-4 orders of magnitude faster than the FDTD method and enables systematic investigations of the effects of scatterer size and configuration on the focal field. We demonstrate the power of the new HF-WEFS approach by mapping several metrics of focal field distortion as a function of scatterer position. This analysis shows that the maximum focal field distortion occurs for single scatterers placed below the focal plane with an offset from the optical axis. The HF-WEFS method represents an important first step toward the development of a computational model of laser-scanning microscopy of thick cellular/tissue specimens.

Coupled forward-adjoint Monte Carlo simulation of spatial-angular light fields to determine optical sensitivity in turbid media.

We present a coupled forward-adjoint Monte Carlo (cFAMC) method to determine the spatially resolved sensitivity distributions produced by optical interrogation of three-dimensional (3-D) tissue volumes. We develop a general computational framework that computes the spatial and angular distributions of the forward-adjoint light fields to provide accurate computations in mesoscopic tissue volumes. We provide full computational details of the cFAMC method and provide results for low- and high-scattering tissues probed using a single pair of optical fibers. We examine the effects of source-detector separation and orientation on the sensitivity distributions and consider how the degree of angular discretization used in the 3-D tissue model impacts the accuracy of the resulting absorption sensitivity profiles. We discuss the value of such computations for optical imaging and the design of optical measurements.

Comparative analysis of discrete and continuous absorption weighting estimators used in Monte Carlo simulations of radiative transport in turbid media.

We examine the relative error of Monte Carlo simulations of radiative transport that employ two commonly used estimators that account for absorption differently, either discretely, at interaction points, or continuously, between interaction points. We provide a rigorous derivation of these discrete and continuous absorption weighting estimators within a stochastic model that we show to be equivalent to an analytic model, based on the radiative transport equation (RTE). We establish that both absorption weighting estimators are unbiased and, therefore, converge to the solution of the RTE. An analysis of spatially resolved reflectance predictions provided by these two estimators reveals no advantage to either in cases of highly scattering and highly anisotropic media. However, for moderate to highly absorbing media or isotropically scattering media, the discrete estimator provides smaller errors at proximal source locations while the continuous estimator provides smaller errors at distal locations. The origin of these differing variance characteristics can be understood through examination of the distribution of exiting photon weights.

Perturbation Monte Carlo methods for tissue structure alterations.

This paper describes an extension of the perturbation Monte Carlo method to model light transport when the phase function is arbitrarily perturbed. Current perturbation Monte Carlo methods allow perturbation of both the scattering and absorption coefficients, however, the phase function can not be varied. The more complex method we develop and test here is not limited in this way. We derive a rigorous perturbation Monte Carlo extension that can be applied to a large family of important biomedical light transport problems and demonstrate its greater computational efficiency compared with using conventional Monte Carlo simulations to produce forward transport problem solutions. The gains of the perturbation method occur because only a single baseline Monte Carlo simulation is needed to obtain forward solutions to other closely related problems whose input is described by perturbing one or more parameters from the input of the baseline problem. The new perturbation Monte Carlo methods are tested using tissue light scattering parameters relevant to epithelia where many tumors originate. The tissue model has parameters for the number density and average size of three classes of scatterers; whole nuclei, organelles such as lysosomes and mitochondria, and small particles such as ribosomes or large protein complexes. When these parameters or the wavelength is varied the scattering coefficient and the phase function vary. Perturbation calculations give accurate results over variations of ∼15-25% of the scattering parameters.

In vivo multiphoton NADH fluorescence reveals depth-dependent keratinocyte metabolism in human skin.

We employ a clinical multiphoton microscope to monitor in vivo and noninvasively the changes in reduced nicotinamide adenine dinucleotide (NADH) fluorescence of human epidermal cells during arterial occlusion. We correlate these results with measurements of tissue oxy- and deoxyhemoglobin concentration during oxygen deprivation using spatial frequency domain imaging. During arterial occlusion, a decrease in oxyhemoglobin corresponds to an increase in NADH fluorescence in the basal epidermal cells, implying a reduction in basal cell oxidative phosphorylation. The ischemia-induced oxygen deprivation is associated with a strong increase in NADH fluorescence of keratinocytes in layers close to the stratum basale, whereas keratinocytes from epidermal layers closer to the skin surface are not affected. Spatial frequency domain imaging optical property measurements, combined with a multilayer Monte Carlo-based radiative transport model of multiphoton microscopy signal collection in skin, establish that localized tissue optical property changes during occlusion do not impact the observed NADH signal increase. This outcome supports the hypothesis that the vascular contribution to the basal layer oxygen supply is significant and these cells engage in oxidative metabolism. Keratinocytes in the more superficial stratum granulosum are either supplied by atmospheric oxygen or are functionally anaerobic. Based on combined hemodynamic and two-photon excited fluorescence data, the oxygen consumption rate in the stratum basale is estimated to be ∼0.035 µmoles/10(6) cells/h.

Quantitative, depth-resolved determination of particle motion using multi-exposure, spatial frequency domain laser speckle imaging.

Laser Speckle Imaging (LSI) is a simple, noninvasive technique for rapid imaging of particle motion in scattering media such as biological tissue. LSI is generally used to derive a qualitative index of relative blood flow due to unknown impact from several variables that affect speckle contrast. These variables may include optical absorption and scattering coefficients, multi-layer dynamics including static, non-ergodic regions, and systematic effects such as laser coherence length. In order to account for these effects and move toward quantitative, depth-resolved LSI, we have developed a method that combines Monte Carlo modeling, multi-exposure speckle imaging (MESI), spatial frequency domain imaging (SFDI), and careful instrument calibration. Monte Carlo models were used to generate total and layer-specific fractional momentum transfer distributions. This information was used to predict speckle contrast as a function of exposure time, spatial frequency, layer thickness, and layer dynamics. To verify with experimental data, controlled phantom experiments with characteristic tissue optical properties were performed using a structured light speckle imaging system. Three main geometries were explored: 1) diffusive dynamic layer beneath a static layer, 2) static layer beneath a diffuse dynamic layer, and 3) directed flow (tube) submerged in a dynamic scattering layer. Data fits were performed using the Monte Carlo model, which accurately reconstructed the type of particle flow (diffusive or directed) in each layer, the layer thickness, and absolute flow speeds to within 15% or better.

Analysis of single Monte Carlo methods for prediction of reflectance from turbid media.

Starting from the radiative transport equation we derive the scaling relationships that enable a single Monte Carlo (MC) simulation to predict the spatially- and temporally-resolved reflectance from homogeneous semi-infinite media with arbitrary scattering and absorption coefficients. This derivation shows that a rigorous application of this single Monte Carlo (sMC) approach requires the rescaling to be done individually for each photon biography. We examine the accuracy of the sMC method when processing simulations on an individual photon basis and also demonstrate the use of adaptive binning and interpolation using non-uniform rational B-splines (NURBS) to achieve order of magnitude reductions in the relative error as compared to the use of uniform binning and linear interpolation. This improved implementation for sMC simulation serves as a fast and accurate solver to address both forward and inverse problems and is available for use at http://www.virtualphotonics.org/.

Electric field Monte Carlo simulations of focal field distributions produced by tightly focused laser beams in tissues.

The focal field distribution of tightly focused laser beams in turbid media is sensitive to optical scattering and therefore of direct relevance to image quality in confocal and nonlinear microscopy. A model that considers both the influence of scattering and diffraction on the amplitude and phase of the electric field in focused beam geometries is required to describe these distorted focal fields. We combine an electric field Monte Carlo approach that simulates the electric field propagation in turbid media with an angular-spectrum representation of diffraction theory to analyze the effect of tissue scattering properties on the focal field. In particular, we examine the impact of variations in the scattering coefficient (µ(s)), single-scattering anisotropy (g), of the turbid medium and the numerical aperture of the focusing lens on the focal volume at various depths. The model predicts a scattering-induced broadening, amplitude loss, and depolarization of the focal field that corroborates experimental results. We find that both the width and the amplitude of the focal field are dictated primarily by µ(s) with little influence from g. In addition, our model confirms that the depolarization rate is small compared to the amplitude loss of the tightly focused field.

Investigation of a probe design for facilitating the uses of the standard photon diffusion equation at short source-detector separations: Monte Carlo simulations.

We design a special diffusing probe to investigate the optical properties of human skin in vivo. The special geometry of the probe enables a modified two-layer (MTL) diffusion model to precisely describe the photon transport even when the source-detector separation is shorter than 3 mean free paths. We provide a frequency domain comparison between the Monte Carlo model and the diffusion model in both the MTL geometry and conventional semiinfinite geometry. We show that using the Monte Carlo model as a benchmark method, the MTL diffusion theory performs better than the diffusion theory in the semiinfinite geometry. In addition, we carry out Monte Carlo simulations with the goal of investigating the dependence of the interrogation depth of this probe on several parameters including source-detector separation, sample optical properties, and properties of the diffusing high-scattering layer. From the simulations, we find that the optical properties of samples modulate the interrogation volume greatly, and the source-detector separation and the thickness of the diffusing layer are the two dominant probe parameters that impact the interrogation volume. Our simulation results provide design guidelines for a MTL geometry probe.

Amplitude and phase of tightly focused laser beams in turbid media.

A framework is developed that combines electric field Monte Carlo simulations of random scattering with an angular-spectrum representation of diffraction theory to determine the amplitude and phase characteristics of tightly focused laser beams in turbid media. For planar sample geometries, the scattering-induced coherence loss of wave vectors at larger angles is shown to be the primary mechanism for broadening the focal volume. This approach for evaluating the formation of the focal volume in turbid media is of direct relevance to the imaging properties of nonlinear coherent microscopy, which rely on both the amplitude and phase of the focused fields.

Radiative transport in the delta-P1 approximation for semi-infinite turbid media.

We have developed an analytic solution for spatially resolved diffuse reflectance within the deltaP1 approximation to the radiative transport equation for a semi-infinite homogeneous turbid medium. We evaluate the performance of this solution by comparing its predictions with those provided by Monte Carlo simulations and the standard diffusion approximation. We demonstrate that the delta-P1 approximation provides accurate estimates for spatially resolved diffuse reflectance in both low and high scattering media. We also develop a multi-stage nonlinear optimization algorithm in which the radiative transport estimates provided by the delta-P1 approximation are used to recover the optical absorption (microa), reduced scattering (micros'), and single-scattering asymmetry coefficients (g1) of liquid and solid phantoms from experimental measurements of spatially resolved diffuse reflectance. Specifically, the delta-P1 approximation can be used to recover microa, micros', and g1 with errors within +/- 22%, +/- 18%, and +/- 17%, respectively, for both intralipid-based and siloxane-based tissue phantoms. These phantoms span the optical property range 4 < (micros' /microa) < 117. Using these same measurements, application of the standard diffusion approximation resulted in the recovery of microa and micros' with errors o f +/- 29% and +/- 25%, respectively. Collectively, these results demonstrate that the delta-P1 approximation provides accurate radiative transport estimates that can be used to determine accurately the optical properties of biological tissues, particularly in spectral regions where tissue may display moderate/low ratios of reduced scattering to absorption (micros'/microa).

Determination of optical properties of superficial volumes of layered tissue phantoms.

Previously, we reported the design of a new diffusing probe that employs a standard two-layer diffusion model to recover the optical properties of turbid samples. This particular probe had a source-detector separation of 2.5 mm and performance was validated with Monte Carlo simulations and homogeneous phantom experiments. The goal of the current study is to characterize the performance of this new method in the context of two-layer phantoms that mimic the optical properties of human skin. We analyze the accuracy of the recovered top layer optical properties and their dependences on the thickness of the top layer of two-layer phantoms. Our results demonstrate that the optical properties of the top layer can be accurately determined with a 1.6 mm source-detector separation diffusing probe when this layer thickness is as thin as 1 mm. Monte Carlo simulations illustrate that the interrogation depth can be further decreased by shortening the source-detector separation.

Method for recovering quantitative broadband diffuse optical spectra from layered media.

We report the recovery of broadband (650-1000 nm) diffuse optical absorption and reduced scattering spectra stratified by layer in a two-layer phantom. The broadband optical properties of the phantom featured top and bottom layers designed to simulate adipose and muscle, respectively. The absolute value and dynamic variation of chromophore concentrations in both layers (top layer thickness greater than 5 mm) were calculated with an average 10% error and 3% error, respectively. In addition to spectra, the algorithm recovers the top layer thickness up to 12 mm within 10% error. It is insensitive to initial guesses of both layers' optical properties as long as the layer thickness initial guess is within +/-2 mm.