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.

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.

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.

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.

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.

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.

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

Perturbation and differential Monte Carlo methods for measurement of optical properties in a layered epithelial tissue model

The use of perturbation and differential Monte Carlo (pMC/dMC) methods in conjunction with nonlinear optimization algorithms were proposed recently as a means to solve inverse photon migration problems in regionwise heterogeneous turbid media. We demonstrate the application of pMC/dMC methods for the recovery of optical properties in a two-layer extended epithelial tissue model from experimental measurements of spatially resolved diffuse reflectance. The results demonstrate that pMC/dMC methods provide a rapid and accurate approach to solve two-region inverse photon migration problems in the transport regime, that is, on spatial scales smaller than a transport mean free path and in media where optical scattering need not dominate absorption. The pMC/dMC approach is found to be effective over a broad range of absorption (50 to 400%) and scattering (70 to 130%) perturbations. The recovery of optical properties from spatially resolved diffuse reflectance measurements is examined for different sets of source-detector separation. These results provide some guidance for the design of compact fiber-based probes to determine and isolate optical properties from both epithelial and stromal layers of superficial tissues.

Frequency domain photon migration in the delta- P1 approximation: analysis of ballistic, transport, and diffuse regimes

The standard diffusion approximation (SDA) to the Boltzmann transport equation (BTE) is commonly used to describe radiative transport for biomedical applications of frequency-domain diffuse optical imaging and spectroscopy. Unfortunately, the SDA is unable to provide accurate radiative transport predictions on spatial scales comparable to the transport mean free path and for media in which optical scattering is not dominant over absorption. Here, we develop and demonstrate the use of the delta- P1 approximation to provide improved radiative transport estimates in the frequency domain via the addition of a Dirac delta function to both radiance and phase function approximations. Specifically, we consider photon density wave propagation resulting from the illumination of an infinite turbid medium with an embedded, intensity-modulated, spherical light source. We examine the accuracy of the standard diffusion and delta- P1 approximations relative to Monte Carlo simulations that provide exact solutions to the BTE. This comparison establishes the superior accuracy of the delta- P1 approximation relative to the SDA that is most notable at distances less than 3 transport mean free paths from the source. In addition, we demonstrate that the differences in photon density wave propagation in a highly forward scattering medium (g1=0.95) vs an isotropically scattering medium (g1=0) provides a basis to define three spatial regimes where the light field is dominated by (a) unscattered/ballistic light, (b) minimally scattered light, and (c) diffusely scattered light. We examine the impact of optical properties, source modulation frequency, and numerical aperture of detection on the spatial extent and location of these regimes.

Use of the delta-P1 approximation for recovery of optical absorption, scattering, and asymmetry coefficients in turbid media

We introduce a robust method to recover optical absorption, reduced scattering, and single-scattering asymmetry coefficients (microa, micro's, g1) of infinite turbid media over a range of (micro's/microa) spanning 3 orders of magnitude. This is accomplished through the spatially resolved measurement of irradiance at source-detector separations spanning 0.25-8 transport mean free paths (l*). These measurements are rapidly processed by a multistaged nonlinear optimization algorithm in which the measured irradiances are compared with predictions given by the delta-P1 variant of the diffusion approximation to the Boltzmann transport equation. The ability of the delta-P1 model to accurately describe radiative transport within media of arbitrary albedo and on spatial scales comparable to l* is the key element enabling the separation of g1 from micro's.

Sampling tissue volumes using frequency-domain photon migration

We demonstrate the use of Monte Carlo simulations to generate photon scattering density functions (PSDFs) that represent the tissue volume sampled by steady-state and frequency-domain photon migration. We use these results to illustrate how scaling laws can be developed to determine the mean sampling depth of the multiply scattered photons detected by photon migration methods that remain valid outside the bounds of the standard diffusion approximation, i.e., at small source-detector separations and in media where the optical absorption is significant relative to scattering. Using both the PSDF computation and the newly formulated scaling laws, we focus on a comprehensive description of the effects of source modulation frequency, optical absorption, and source-detector separation on the depth of the sampled tissue volume as well as the sensitivity of frequency-domain photon migration measurements to the presence of a localized absorption heterogeneity.

Perturbation Monte Carlo methods to solve inverse photon migration problems in heterogeneous tissues

We introduce a novel and efficient method to provide solutions to inverse photon migration problems in heterogeneous turbid media. The method extracts derivative information from a single Monte Carlo simulation to permit the rapid determination of rates of change in the detected photon signal with respect to perturbations in background tissue optical properties. We then feed this derivative information to a nonlinear optimization algorithm to determine the optical properties of the tissue heterogeneity under examination. We demonstrate the use of this approach to solve rapidly a two-region inverse problem of photon migration in the transport regime, for which diffusion-approximation-based approaches are not applicable.