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

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].

Hybrid Monte Carlo estimators for multilayer transport problems

This paper introduces a new family of hybrid estimators aimed at controlling the efficiency of Monte Carlo computations in particle transport problems. In this context, efficiency is usually measured by the figure of merit (FOM) given by the inverse product of the estimator variance Var[ξ] and the run time T: FOM := (Var[ξ] T)^-1. Previously, we developed a new family of transport-constrained unbiased radiance estimators (T-CURE) that generalize the conventional collision and track length estimators [1] and provide 1-2 orders of magnitude additional variance reduction. However, these gains in variance reduction are partly offset by increases in overhead time [2], lowering their computational efficiency. Here we show that combining T-CURE estimation with conventional terminal estimation within each individual biography can moderate the efficiency of the resulting "hybrid" estimator without introducing bias in the computation. This is achieved by treating only the refractive interface crossings with the extended next event estimator, and all others by standard terminal estimators. This is because when there are index-mismatched interfaces between the collision location and the detector, the T-CURE computation rapidly becomes intractable due to the large number of refractions and reflections that can arise. We illustrate the gains in efficiency by comparing our hybrid strategy with more conventional estimation methods in a series of multi-layer numerical examples.

Perturbation Monte Carlo Method for Quantitative Photoacoustic Tomography

Quantitative photoacoustic tomography aims at estimating optical parameters from photoacoustic images that are formed utilizing the photoacoustic effect caused by the absorption of an externally introduced light pulse. This optical parameter estimation is an ill-posed inverse problem, and thus it is sensitive to measurement and modeling errors. In this work, we propose a novel way to solve the inverse problem of quantitative photoacoustic tomography based on the perturbation Monte Carlo method. Monte Carlo method for light propagation is a stochastic approach for simulating photon trajectories in a medium with scattering particles. It is widely accepted as an accurate method to simulate light propagation in tissues. Furthermore, it is numerically robust and easy to implement. Perturbation Monte Carlo maintains this robustness and enables forming gradients for the solution of the inverse problem. We validate the method and apply it in the framework of Bayesian inverse problems. The simulations show that the perturbation Monte Carlo method can be used to estimate spatial distributions of both absorption and scattering parameters simultaneously. These estimates are qualitatively good and quantitatively accurate also in parameter scales that are realistic for biological tissues.

Predictive model for the quantitative analysis of human skin using photothermal radiometry and diffuse reflectance spectroscopy

We have recently introduced a novel methodology for the noninvasive analysis of the structure and composition of human skin in vivo. The approach combines pulsed photothermal radiometry (PPTR), involving time-resolved measurements of mid-infrared emission after irradiation with a millisecond light pulse, and diffuse reflectance spectroscopy (DRS) in the visible part of the spectrum. Simultaneous fitting of both data sets with respective predictions from a numerical model of light transport in human skin enables the assessment of the contents of skin chromophores (melanin, oxy-, and deoxy-hemoglobin), as well as scattering properties and thicknesses of the epidermis and dermis. However, the involved iterative optimization of 14 skin model parameters using a numerical forward model (i.e., inverse Monte Carlo - IMC) is computationally very expensive. In order to overcome this drawback, we have constructed a very fast predictive model (PM) based on machine learning. The PM involves random forests, trained on ∼9,000 examples computed using our forward MC model. We show that the performance of such a PM is very satisfying, both in objective testing using cross-validation and in direct comparisons with the IMC procedure. We also present a hybrid approach (HA), which combines the speed of the PM with versatility of the IMC procedure. Compared with the latter, the HA improves both the accuracy and robustness of the inverse analysis, while significantly reducing the computation times.

Incoherent scattering from dielectric metasurfaces under the influence of electromagnetic eigenmodes

Anomalous optical properties of microscopically inhomogeneous dielectric films placed on a thick metal sublayer are investigated. We study the reflection, scattering, and absorption of the coherent electromagnetic radiation as a function of the incidence angle. Computer simulations show the existence of the incidence angle of the laser beam when the scattering and absorption increase simultaneously for the s-polarization so that almost 60% of the incident light goes in the scattering channel. The critical angle corresponds to the excitation of Fabry-Perot mode. The effect makes it possible to manipulate the reflection from the metafilms.

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.

Comparison of speckleplethysmographic (SPG) and photoplethysmographic (PPG) imaging by Monte Carlo simulations and in vivo measurements

Noncontact photoplethysmography (PPG) is limited by a poor signal-to-noise ratio (SNR). A solution to this limitation is the use of alternate sources of optical contrast to generate a complementary pulsatile waveform. One such source is laser speckle contrast, which is modulated in biological tissues by the flow rate of red blood cells. Averaging a region of interest from a speckle contrast image over time allows for the calculation of a speckleplethysmogram (SPG). Similar to PPG, SPG enables monitoring of heart rate and respiratory rate. A gap in the knowledge base exists as to the precise spatiotemporal relationship between PPG and SPG signals. We have developed an eight-layer tissue model to simulate both PPG and SPG signals in a reflectance geometry via Monte Carlo methods. We modeled PPG by compression of the upper and lower blood nets due to expansion of the larger arterial layer below. The in silico PPG peak-to-peak amplitude percent was greater at 532 nm than at 860 nm (5.6% vs. 3.0%, respectively), which matches trends from the literature. We modeled SPG by changing flow speeds of red blood cells in both the capillaries and arterioles over the cardiac cycle. The in silico SPG peak-to-peak amplitude percent was 24% at 532 nm and 40% at 860 nm. In silico results are similar to in vivo results measured with a two-camera set up for simultaneous imaging of PPG and SPG. Both in silico and in vivo data suggest SPG has a much larger SNR than PPG, which may prove beneficial for noncontact, wide-field optical monitoring of cardiovascular health.

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.

Hyperspectral optical tomography of intrinsic signals in the rat cortex

We introduce a tomographic approach for three-dimensional imaging of evoked hemodynamic activity, using broadband illumination and diffuse optical tomography (DOT) image reconstruction. Changes in diffuse reflectance in the rat somatosensory cortex due to stimulation of a single whisker were imaged at a frame rate of 5 Hz using a hyperspectral image mapping spectrometer. In each frame, images in 38 wavelength bands from 484 to 652 nm were acquired simultaneously. For data analysis, we developed a hyperspectral DOT algorithm that used the Rytov approximation to quantify changes in tissue concentration of oxyhemoglobin ([Formula: see text]) and deoxyhemoglobin (ctHb) in three dimensions. Using this algorithm, the maximum changes in [Formula: see text] and ctHb were found to occur at [Formula: see text] and [Formula: see text] beneath the surface of the cortex, respectively. Rytov tomographic reconstructions revealed maximal spatially localized increases and decreases in [Formula: see text] and ctHb of [Formula: see text] and [Formula: see text], respectively, with these maximum changes occurring at [Formula: see text] poststimulus. The localized optical signals from the Rytov approximation were greater than those from modified Beer-Lambert, likely due in part to the inability of planar reflectance to account for partial volume effects.

Extending generalized Kubelka-Munk to three-dimensional radiative transfer

The generalized Kubelka-Munk (gKM) approximation is a linear transformation of the double spherical harmonics of order one (DP₁) approximation of the radiative transfer equation. Here, we extend the gKM approximation to study problems in three-dimensional radiative transfer. In particular, we derive the gKM approximation for the problem of collimated beam propagation and scattering in a plane-parallel slab composed of a uniform absorbing and scattering medium. The result is an 8×8 system of partial differential equations that is much easier to solve than the radiative transfer equation. We compare the solutions of the gKM approximation with Monte Carlo simulations of the radiative transfer equation to identify the range of validity for this approximation. We find that the gKM approximation is accurate for isotropic scattering media that are sufficiently thick and much less accurate for anisotropic, forward-peaked scattering media.

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.

Design and validation of a diffuse reflectance and spectroscopic microendoscope with poly(dimethylsioxane)-based phantoms

Many cases of epithelial cancer originate in basal layers of tissue and are initially undetected by conventional microendoscopy techniques. We present a bench-top, fiber-bundle microendoscope capable of providing high resolution images of surface cell morphology. Additionally, the microendoscope has the capability to interrogate deeper into material by using diffuse reflectance and broadband diffuse reflectance spectroscopy. The purpose of this multimodal technique was to overcome the limitation of microendoscopy techniques that are limited to only visualizing morphology at the tissue or cellular level. Using a custom fiber optic probe, high resolution surface images were acquired using topical proflavine to fluorescently stain non-keratinized epithelia. A 635 nm laser coupled to a 200 μm multimode fiber delivers light to the sample and the diffuse reflectance signal was captured by a 1 mm image guide fiber. Finally, a tungsten-halogen lamp coupled to a 200 μm multimode fiber delivers broadband light to the sample to acquire spectra at source-detector separations of 374, 729, and 1051 μm. To test the instrumentation, a high resolution proflavine-induced fluorescent image of resected healthy mouse colon was acquired. Additionally, five monolayer poly(dimethylsiloxane)-based optical phantoms with varying absorption and scattering properties were created to acquire diffuse reflectance profiles and broadband spectra.

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.