Condensed Monte Carlo simulations for the description of light transport

Optical Goniometer Paired with Digital Monte Carlo Twin to Determine the Optical Properties of Turbid Media

We present a goniometer designed for capturing spectral and angular-resolved data from scattering and absorbing media. The experimental apparatus is complemented by a comprehensive Monte Carlo simulation, meticulously replicating the radiative transport processes within the instrument's optical components and simulating scattering and absorption across arbitrary volumes. Consequently, we were able to construct a precise digital replica, or "twin", of the experimental setup. This digital counterpart enabled us to tackle the inverse problem of deducing optical parameters such as absorption and scattering coefficients, along with the scattering anisotropy factor from measurements. We achieved this by fitting Monte Carlo simulations to our goniometric measurements using a Levenberg-Marquardt algorithm. Validation of our approach was performed using polystyrene particles, characterized by Mie scattering, supplemented by a theoretical analysis of algorithmic convergence. Ultimately, we demonstrate strong agreement between optical parameters derived using our novel methodology and those obtained via established measurement protocols.

Efficient algorithm to calculate the optical properties of breast tumors by high-order perturbation theory

An efficient algorithm to obtain the solutions for n-th order terms of perturbation expansions in absorption, scattering, and cross-coupling for light propagating in human tissue is presented. The proposed solution is free of any approximations and makes possible fast and efficient estimates of mammographic, optical tomographic, and fluorescent images, applying a perturbation order of 30 and more. The presented analysis sets the general limits for the applicability of the perturbation approach as a function of tumor size and optical properties of the human tissue. The convergence tests of the efficient calculations for large absorbing objects show excellent agreement with the reference data from finite element method calculations. The applicability of the theory is demonstrated in experiments on breast-like phantoms with high absorbing and low-scattering lesions.

Improved topographic reconstruction of turbid media in the spatial frequency domain including the determination of the reduced scattering and absorption coefficients

The separation of scattering and absorption is of great importance for studying the radiative transfer in turbid media. Obtaining the corresponding coefficients for non-flat objects is difficult and needs special consideration. Building on our previous work [J. Opt. Soc. Am. A39, 1823 (2022)JOAOD60740-323210.1364/JOSAA.464007], we present an approach that takes the changing incident and detection angles relative to the surface normal of curved surfaces into account to improve the determination of the reduced scattering and absorption coefficients with measurements in the spatial frequency domain (SFD). The optical coefficients are reconstructed using a pre-calculated lookup table generated with Monte Carlo simulations on graphical processing units. With the obtained values, the error in the captured surface geometry of the object, which is due to the volume scattering, is compensated and reduced by 1 order of magnitude for measurements in the SFD. Considering the approximate surface geometry, the absorption and reduced scattering are accurately resolved for moderate object curvatures, with very low dependence on the tilt angle. In contrast to models that only correct the amplitudes of the SFD signal, our approach, in addition to the optical properties, predicts the phase values correctly, which is the reason why it can be used to correct the surface geometry.

Compact setup to determine size and concentration of spherical particles in a turbid medium

We propose a compact setup to determine the size and concentration of spherical particles in a turbid medium. A pair of plane mirrors is used to multifold the undeviated laser beam, and measure it at a detector placed close to the sample, to determine the interaction coefficient. The size of particles is uniquely determined by comparison of the scattered light from the medium, measured at two separate detectors placed at two different angular positions, with that from Monte Carlo simulations. The methodology is verified using measurements with turbid samples comprising polystyrene spheres.

In-silico investigation towards the non-invasive optical detection of blood lactate

This paper uses Monte Carlo simulations to investigate the interaction of short-wave infrared (SWIR) light with vascular tissue as a step toward the development of a non-invasive optical sensor for measuring blood lactate in humans. The primary focus of this work was to determine the optimal source-detector separation, penetration depth of light at SWIR wavelengths in tissue, and the optimal light power required for reliable detection of lactate. The investigation also focused on determining the non-linear variations in absorbance of lactate at a few select SWIR wavelengths. SWIR photons only penetrated 1.3 mm and did not travel beyond the hypodermal fat layer. The maximum output power was only 2.51% of the input power, demonstrating the need for a highly sensitive detection system. Simulations optimized a source-detector separation of 1 mm at 1684 nm for accurate measurement of lactate in blood.

Evaluation of Fluence Correction Algorithms in Multispectral Photoacoustic Imaging

Multispectral photoacoustic imaging (MPAI) is a promising emerging diagnostic technology, but fluence artifacts can degrade device performance. Our goal was to develop well-validated phantom-based test methods for evaluating and comparing MPAI fluence correction algorithms, including a heuristic diffusion approximation, Monte Carlo simulations, and an algorithm we developed based on novel application of the diffusion dipole model (DDM). Phantoms simulated a range of breast-mimicking optical properties and contained channels filled with chromophore solutions (ink, hemoglobin, or copper sulfate) or connected to a previously developed blood flow circuit providing tunable oxygen saturation (SO₂). The DDM algorithm achieved similar spectral recovery and SO₂ measurement accuracy to Monte Carlo-based corrections with lower computational cost, potentially providing an accurate, real-time correction approach. Algorithms were sensitive to optical property uncertainty, but error was minimized by matching phantom albedo. The developed test methods may provide a foundation for standardized assessment of MPAI fluence correction algorithm performance.

Efficient inversion strategies for estimating optical properties with Monte Carlo radiative transport models

SIGNIFICANCE

Indirect imaging problems in biomedical optics generally require repeated evaluation of forward models of radiative transport, for which Monte Carlo is accurate yet computationally costly. We develop an approach to reduce this bottleneck, which has significant implications for quantitative tomographic imaging in a variety of medical and industrial applications.

AIM

Our aim is to enable computationally efficient image reconstruction in (hybrid) diffuse optical modalities using stochastic forward models.

APPROACH

Using Monte Carlo, we compute a fully stochastic gradient of an objective function for a given imaging problem. Leveraging techniques from the machine learning community, we then adaptively control the accuracy of this gradient throughout the iterative inversion scheme to substantially reduce computational resources at each step.

RESULTS

For example problems of quantitative photoacoustic tomography and ultrasound-modulated optical tomography, we demonstrate that solutions are attainable using a total computational expense that is comparable to (or less than) that which is required for a single high-accuracy forward run of the same Monte Carlo model.

CONCLUSIONS

This approach demonstrates significant computational savings when approaching the full nonlinear inverse problem of optical property estimation using stochastic methods.

Methods of extraction of optical properties from diffuse reflectance measurements of ex-vivo human colon tissue using thin film silicon photodetector arrays

Spatially resolved diffuse reflectance spectroscopy (SRDRS) is a promising technique for characterization of colon tissue. Herein, two methods for extracting the reduced scattering and absorption coefficients ( μ s ' ( λ ) and μ a ( λ ) ) from SRDRS data using lookup tables of simulated diffuse reflectance are reported. Experimental measurements of liquid tissue phantoms performed with a custom multi-pixel silicon SRDRS sensor spanning the 450 - 750 nm wavelength range were used to evaluate the extraction methods, demonstrating that the combined use of spatial and spectral data reduces extraction error compared to use of spectral data alone. Additionally, SRDRS measurements of normal and tumor ex-vivo human colon tissue are presented along with μ s ' ( λ ) and μ a ( λ ) extracted from these measurements.

High-density diffuse optical tomography for imaging human brain function

This review describes the unique opportunities and challenges for noninvasive optical mapping of human brain function. Diffuse optical methods offer safe, portable, and radiation free alternatives to traditional technologies like positron emission tomography or functional magnetic resonance imaging (fMRI). Recent developments in high-density diffuse optical tomography (HD-DOT) have demonstrated capabilities for mapping human cortical brain function over an extended field of view with image quality approaching that of fMRI. In this review, we cover fundamental principles of the diffusion of near infrared light in biological tissue. We discuss the challenges involved in the HD-DOT system design and implementation that must be overcome to acquire the signal-to-noise necessary to measure and locate brain function at the depth of the cortex. We discuss strategies for validation of the sensitivity, specificity, and reliability of HD-DOT acquired maps of cortical brain function. We then provide a brief overview of some clinical applications of HD-DOT. Though diffuse optical measurements of neurophysiology have existed for several decades, tremendous opportunity remains to advance optical imaging of brain function to address a crucial niche in basic and clinical neuroscience: that of bedside and minimally constrained high fidelity imaging of brain function.

Hybrid method to estimate two-layered superficial tissue optical properties from simulated data of diffuse reflectance spectroscopy

An iterative curve fitting method has been applied in both simulation [J. Biomed. Opt.17, 107003 (2012)JBOPFO1083-366810.1117/1.JBO.17.10.107003] and phantom [J. Biomed. Opt.19, 077002 (2014)JBOPFO1083-366810.1117/1.JBO.19.7.077002] studies to accurately extract optical properties and the top layer thickness of a two-layered superficial tissue model from diffuse reflectance spectroscopy (DRS) data. This paper describes a hybrid two-step parameter estimation procedure to address two main issues of the previous method, including (1) high computational intensity and (2) converging to local minima. The parameter estimation procedure contained a novel initial estimation step to obtain an initial guess, which was used by a subsequent iterative fitting step to optimize the parameter estimation. A lookup table was used in both steps to quickly obtain reflectance spectra and reduce computational intensity. On simulated DRS data, the proposed parameter estimation procedure achieved high estimation accuracy and a 95% reduction of computational time compared to previous studies. Furthermore, the proposed initial estimation step led to better convergence of the following fitting step. Strategies used in the proposed procedure could benefit both the modeling and experimental data processing of not only DRS but also related approaches such as near-infrared spectroscopy.

Review of mesoscopic optical tomography for depth-resolved imaging of hemodynamic changes and neural activities

Understanding the functional wiring of neural circuits and their patterns of activation following sensory stimulations is a fundamental task in the field of neuroscience. Furthermore, charting the activity patterns is undoubtedly important to elucidate how neural networks operate in the living brain. However, optical imaging must overcome the effects of light scattering in the tissue, which limit the light penetration depth and affect both the imaging quantitation and sensitivity. Laminar optical tomography (LOT) is a three-dimensional (3-D) in-vivo optical imaging technique that can be used for functional imaging. LOT can achieve both a resolution of 100 to [Formula: see text] and a penetration depth of 2 to 3 mm based either on absorption or fluorescence contrast, as well as large field-of-view and high acquisition speed. These advantages make LOT suitable for 3-D depth-resolved functional imaging of the neural functions in the brain and spinal cords. We review the basic principles and instrumentations of representative LOT systems, followed by recent applications of LOT on 3-D imaging of neural activities in the rat forepaw stimulation model and mouse whisker-barrel system.

Estimation of optical properties by spatially resolved reflectance spectroscopy in the subdiffusive regime

We propose and objectively evaluate an inverse Monte Carlo model for estimation of absorption and reduced scattering coefficients and similarity parameter ? from spatially resolved reflectance (SRR) profiles in the subdiffusive regime. The similarity parameter ? carries additional information on the phase function that governs the angular properties of scattering in turbid media. The SRR profiles at five source-detector separations were acquired with an optical fiber probe. The inverse Monte Carlo model was based on a cost function that enabled robust estimation of optical properties from a few SRR measurements without a priori knowledge about spectral dependencies of the optical properties. Validation of the inverse Monte Carlo model was performed on synthetic datasets and measured SRR profiles of turbid phantoms comprising molecular dye and polystyrene microspheres. We observed that the additional similarity parameter ? substantially reduced the reflectance variability arising from the phase function properties and significantly improved the accuracy of the inverse Monte Carlo model. However, the observed improvement of the extended inverse Monte Carlo model was limited to reduced scattering coefficients exceeding ?15??cm?1, where the relative root-mean-square errors of the estimated optical properties were well within 10%.

Estimation of bulk optical properties of turbid media from hyperspectral scatter imaging measurements: metamodeling approach

In many research areas and application domains, the bulk optical properties of biological materials are of great interest. Unfortunately, these properties cannot be obtained easily for complex turbid media. In this study, a metamodeling approach has been proposed and applied for the fast and accurate estimation of the bulk optical properties from contactless and non-destructive hyperspectral scatter imaging (HSI) measurements. A set of liquid optical phantoms, based on intralipid, methylene blue and water, were prepared and the Vis/NIR bulk optical properties were characterized with a double integrating sphere and unscattered transmittance setup. Accordingly, the phantoms were measured with the HSI technique and metamodels were constructed, relating the Vis/NIR reflectance images to the reference bulk optical properties of the samples. The independent inverse validation showed good prediction performance for the absorption coefficient and the reduced scattering coefficient, with R(2)(p) values of 0.980 and 0.998, and RMSE(P) values of 0.032 cm(-1) and 0.197 cm(-1) respectively. The results clearly support the potential of this approach for fast and accurate estimation of the bulk optical properties of turbid media from contactless HSI measurements.

Polarization influence on reflectance measurements in the spatial frequency domain

In this work, we quantify the influence of crossed polarizers on reflectance measurements in the spatial frequency domain. The use of crossed polarizers is a very common approach for suppression of specular surface reflections. However, measurements are typically evaluated using a non-polarized scalar theory. The consequences of this discrepancy are the focus of our study, and we also quantify the related errors of the derived optical properties. We used polarized Monte Carlo simulations for forward calculation of the reflectance from different samples. The samples' scatterers are assumed to be spherical, allowing for the calculation of the scattering functions by Mie theory. From the forward calculations, the reduced scattering coefficient [Formula: see text] and the absorption coefficient μa were derived by means of a scalar theory, as commonly used. Here, we use the analytical solution of the scalar radiative transfer equation. With this evaluation approach, which does not consider polarization, we found large errors in [Formula: see text] and μa in the range of 25% and above. Furthermore, we investigated the applicability of the use of a reference measurement to reduce these errors as suggested in literature. We found that this method is not able to generally improve the accuracy of measurements in the spatial frequency domain. Our general recommendation is to apply a polarized theory when using crossed polarizers.

Optical microsystem for analysis of diffuse reflectance and fluorescence signals applied to early gastrointestinal cancer detection

The detection of cancer at its earliest stage is crucial in order to increase the probability of a successful treatment. Optical techniques, specifically diffuse reflectance and fluorescence, may considerably improve the ability to detect pre-cancerous lesions. These techniques have high sensitivity to some biomarkers present on the tissues, providing morphological and biochemical information of normal and diseased tissue. The development of a chip sized spectroscopy microsystem, based on these techniques, will greatly improve the early diagnosis of gastrointestinal cancers. The main innovation is the detection of the spectroscopic signals using only few, but representative, spectral bands allowing for miniaturization. This paper presents the mathematical models, its validation and analysis for retrieving data of the measured spectroscopic signals. These models were applied to a set of phantoms clearly representative of gastrointestinal tissues, leading to a more accurate diagnostic by a pathologist. Moreover, it was demonstrated that the models can use the reconstructed spectroscopic signals based only on its extraction on those specific spectral bands. As a result, the viability of the spectroscopy microsystem implementation was proved.

Two electric field Monte Carlo models of coherent backscattering of polarized light

Modeling of coherent polarized light propagation in turbid scattering medium by the Monte Carlo method provides an ultimate understanding of coherent effects of multiple scattering, such as enhancement of coherent backscattering and peculiarities of laser speckle formation in dynamic light scattering (DLS) and optical coherence tomography (OCT) diagnostic modalities. In this report, we consider two major ways of modeling the coherent polarized light propagation in scattering tissue-like turbid media. The first approach is based on tracking transformations of the electric field along the ray propagation. The second one is developed in analogy to the iterative procedure of the solution of the Bethe-Salpeter equation. To achieve a higher accuracy in the results and to speed up the modeling, both codes utilize the implementation of parallel computing on NVIDIA Graphics Processing Units (GPUs) with Compute Unified Device Architecture (CUDA). We compare these two approaches through simulations of the enhancement of coherent backscattering of polarized light and evaluate the accuracy of each technique with the results of a known analytical solution. The advantages and disadvantages of each computational approach and their further developments are discussed. Both codes are available online and are ready for immediate use or download.

Design and characterization of a novel multimodal fiber-optic probe and spectroscopy system for skin cancer applications

The design and characterization of an instrument combining Raman, fluorescence, and reflectance spectroscopic modalities is presented. Instrument development has targeted skin cancer applications as a novel fiber-optic probe has been specially designed to interrogate cutaneous lesions. The instrument is modular and both its software and hardware components are described in depth. Characterization of the fiber-optic probe is also presented, which details the probe's ability to measure diagnostically important parameters such as intrinsic fluorescence and absorption and reduced scattering coefficients along with critical performance metrics such as high Raman signal-to-noise ratios at clinically practical exposure times. Validation results using liquid phantoms show that the probe and system can extract absorption and scattering coefficients with less than 10% error. As the goal is to use the instrument for the clinical early detection of skin cancer, preliminary clinical data are also presented, which indicates our system's ability to measure physiological quantities such as relative collagen and nicotinamide adenine dinucleotide concentration, oxygen saturation, blood volume fraction, and mean vessel diameter.

Robust estimation of cerebral hemodynamics in neonates using multilayered diffusion model for normal and oblique incidences

The diffusion approximation is useful for many optical diagnostics modalities, such as near-infrared spectroscopy. However, the simple normal incidence, semi-infinite layer model may prove lacking in estimation of deep-tissue optical properties such as required for monitoring cerebral hemodynamics, especially in neonates. To answer this need, we present an analytical multilayered, oblique incidence diffusion model. Initially, the model equations are derived in vector-matrix form to facilitate fast and simple computation. Then, the spatiotemporal reflectance predicted by the model for a complex neonate head is compared with time-resolved Monte Carlo (TRMC) simulations under a wide range of physiologically feasible parameters. The high accuracy of the multilayer model is demonstrated in that the deviation from TRMC simulations is only a few percent even under the toughest conditions. We then turn to solve the inverse problem and estimate the oxygen saturation of deep brain tissues based on the temporal and spatial behaviors of the reflectance. Results indicate that temporal features of the reflectance are more sensitive to deep-layer optical parameters. The accuracy of estimation is shown to be more accurate and robust than the commonly used single-layer diffusion model. Finally, the limitations of such approaches are discussed thoroughly.

Verification of a two-layer inverse Monte Carlo absorption model using multiple source-detector separation diffuse reflectance spectroscopy

A two-layer Monte Carlo lookup table-based inverse model is validated with two-layered phantoms across physiologically relevant optical property ranges. Reflectance data for source-detector separations of 370 μm and 740 μm were collected from these two-layered phantoms and top layer thickness, reduced scattering coefficient and the top and bottom layer absorption coefficients were extracted using the inverse model and compared to the known values. The results of the phantom verification show that this method is able to accurately extract top layer thickness and scattering when the top layer thickness ranges from 0 to 550 μm. In this range, top layer thicknesses were measured with an average error of 10% and the reduced scattering coefficient was measured with an average error of 15%. The accuracy of top and bottom layer absorption coefficient measurements was found to be highly dependent on top layer thickness, which agrees with physical expectation; however, within appropriate thickness ranges, the error for absorption properties varies from 12-25%.

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.

Monte Carlo lookup table-based inverse model for extracting optical properties from tissue-simulating phantoms using diffuse reflectance spectroscopy

We present a Monte Carlo lookup table (MCLUT)-based inverse model for extracting optical properties from tissue-simulating phantoms. This model is valid for close source-detector separation and highly absorbing tissues. The MCLUT is based entirely on Monte Carlo simulation, which was implemented using a graphics processing unit. We used tissue-simulating phantoms to determine the accuracy of the MCLUT inverse model. Our results show strong agreement between extracted and expected optical properties, with errors rate of 1.74% for extracted reduced scattering values, 0.74% for extracted absorption values, and 2.42% for extracted hemoglobin concentration values.

System analysis of spatial frequency domain imaging for quantitative mapping of surgically resected breast tissues

The feasibility of spatial frequency domain imaging (SFDI) for breast surgical margin assessment was evaluated in tissue-simulating phantoms and in fully intact lumpectomy specimens at the time of surgery. Phantom data was evaluated according to contrast-detail resolution, quantitative accuracy and model-data goodness of fit, where optical parameters were estimated by minimizing the residual sum of squares between the measured modulation amplitude and its solutions, modeled according to diffusion and scaled-Monte Carlo simulations. In contrast-detail phantoms, a 1.25-mm-diameter surface inclusion was detectable for scattering contrast >28%; a fraction of this scattering contrast (7%) was detectable for a 10 mm surface inclusion and at least 33% scattering contrast was detected up to 1.5 mm below the phantom surface, a probing depth relevant to breast surgical margin assessment. Recovered hemoglobin concentrations were insensitive to changes in scattering, except for overestimation at visible wavelengths for total hemoglobin concentrations <15  μM. The scattering amplitude increased linearly with scattering concentration, but the scattering slope depended on both the particle size and number density. Goodness of fit was comparable for the diffusion and scaled-Monte Carlo models of transport in spatially modulated, near-infrared reflectance acquired from 47 lumpectomy tissues, but recovered absorption parameters varied more linearly with expected hemoglobin concentration in liquid phantoms for the scaled-Monte Carlo forward model. SFDI could potentially reduce the high secondary excision rate associated with breast conserving surgery; its clinical translation further requires reduced image reconstruction time and smart inking strategies.

Equivalence of four Monte Carlo methods for photon migration in turbid media

In the field of photon migration in turbid media, different Monte Carlo methods are usually employed to solve the radiative transfer equation. We consider four different Monte Carlo methods, widely used in the field of tissue optics, that are based on four different ways to build photons' trajectories. We provide both theoretical arguments and numerical results showing the statistical equivalence of the four methods. In the numerical results we compare the temporal point spread functions calculated by the four methods for a wide range of the optical properties in the slab and semi-infinite medium geometry. The convergence of the methods is also briefly discussed.

Penetration of light into multiple scattering media: model calculations and reflectance experiments. Part I: the axial transfer

The article presents two general equations of radiation penetration into layers of diffuse reflectors. One of the equations describes the depth origins of reflection, the other the depth profiles of absorption. The equations are evaluated within the theory of radiative transfer applying various degrees of analytical approximations and Monte Carlo simulations. The data are presented for different scattering and absorption coefficients, arbitrary layer thicknesses, collimated and diffused irradiation, and anisotropic forward scattering. The calculated mean depths of reflection are always lower than the mean depths of absorption. For nearly non-absorbing layers, the mean depths of absorption are about one third of the physical layer thickness. In contrast, penetration saturates for strong absorbers at very low depth levels. From the simulated data, methods are derived for the determination of the penetration depth from reflectance and transmittance data of thin layers or from radially diffused reflectance profiles upon spot irradiation. The methods are experimentally verified for a series of metal oxide powders with particle sizes ranging from much smaller to much larger than the wavelength of irradiation and for microcrystalline cellulose stained with different concentrations of an organic dye.

Monte Carlo model of the penetration depth for polarization gating spectroscopy: influence of illumination-collection geometry and sample optical properties

Polarization-gating has been widely used to probe superficial tissue structures, but the penetration depth properties of this method have not been completely elucidated. This study employs a polarization-sensitive Monte Carlo method to characterize the penetration depth statistics of polarization-gating. The analysis demonstrates that the penetration depth depends on both the illumination-collection geometry [illumination-collection area (R) and collection angle (θ(c))] and on the optical properties of the sample, which include the scattering coefficient (μ(s)), absorption coefficient (μ(a)), anisotropy factor (g), and the type of the phase function. We develop a mathematical expression relating the average penetration depth to the illumination-collection beam properties and optical properties of the medium. Finally, we quantify the sensitivity of the average penetration depth to changes in optical properties for different geometries of illumination and collection. The penetration depth model derived in this study can be applied to optimizing application-specific fiber-optic probes to target a sampling depth of interest with minimal sensitivity to the optical properties of the sample.

Modeling diffuse reflectance from homogeneous semi-infinite turbid media for biological tissue applications: a Monte Carlo study

Diffuse reflectance spectroscopy is one of the simplest and widely used techniques for the non-invasive study of biological tissues but no exact analytical solution exists for the problem of diffuse reflectance from turbid media such as biological tissues. In this work, a general treatment of the problem of diffuse reflectance from a homogeneous semi-infinite turbid medium is presented using Monte Carlo simulations. Based on the results of the Monte Carlo method, simple semi-empirical analytical solutions are developed valid for a wide range of collection geometries corresponding to various optical detector diameters. This approach may be useful for the quick and accurate modeling of diffuse reflectance from tissues.

Evaluation of algorithms for microperfusion assessment by fast simulations of laser Doppler power spectral density

In classical laser Doppler (LD) perfusion measurements, zeroth- and first-order moments of the power spectral density of the LD signal are utilized for the calculation of a signal corresponding to the concentration, speed and flow of red blood cells (RBCs). We have analysed the nonlinearities of the moments in relation to RBC speed distributions, parameters of filters utilized in LD instruments and the signal-to-noise ratio. We have developed a new method for fast simulation of the spectrum of the LD signal. The method is based on a superposition of analytically calculated Doppler shift probability distributions derived for the assumed light scattering phase function. We have validated the method by a comparison of the analytically calculated spectra with results of Monte Carlo (MC) simulations. For the semi-infinite, homogeneous medium and the single Doppler scattering regime, the analytical calculation describes LD spectra with the same accuracy as the MC simulation. The method allows for simulating the LD signal in time domain and furthermore analysing the index of perfusion for the assumed wavelength of the light, optical properties of the tissue and concentration of RBCs. Fast simulations of the LD signal in time domain and its frequency spectrum can be utilized in applications where knowledge of the LD photocurrent is required, e.g. in the development of detectors for tissue microperfusion monitoring or in measurements of the LD autocorrelation function for perfusion measurements. The presented fast method for LD spectra calculation can be used as a tool for evaluation of signal processing algorithms used in the LD method and/or for the development of new algorithms of the LD flowmetry and imaging. We analysed LD spectra obtained by analytical calculations using a classical algorithm applied in classical LD perfusion measurements. We observed nonlinearity of the first moment M₁ for low and high speeds of particles (v < 2 mm s⁻¹, v > 10 mm s⁻¹). It was also noted that the first moment M(1) is less sensitive to the change of the mean RBC speed for flat speed distributions. The low-pass filter frequency f₂ implemented in the LD instrument has a significant influence on the first moment of the spectrum. In particular, for a cut-off frequency lower than 10 kHz the M₁ value is strongly underestimated.

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

Multiple scattering model for the penetration depth of low-coherence enhanced backscattering

Low-coherence enhanced backscattering (LEBS) is a depth-selective self-interference phenomenon that originates from light traveling time-reversed paths in a scattering medium. The depth selectivity of LEBS and its sensitivity to optical properties of the scattering medium has made it a promising technique for probing the structure of biological tissue with applications to disease diagnosis and, in particular, precancerous conditions. The ability to accurately predict the penetration depth of the LEBS signal is important in targeting an optimal tissue depth for detecting precancerous cells. This prediction is further complicated by the variation in optical properties of different tissue types. In this paper, the effects of the reduced scattering coefficient (μ(s)'), the phase function and the instrument spatial coherence length (L(sc)) on the LEBS penetration depth are quantified. It is determined that the LEBS penetration depth is primarily dependent on L(sc), μ(s)', and the anisotropy factor (g), but has minimal dependence on higher moments of the phase function. An empirical expression, having a similar form as the double scattering approximation for LEBS, is found to accurately predict the average penetration depth in the multiple scattering regime. The expression is shown to be accurate for a broad range of experimentally relevant optical properties and spatial coherence lengths.

Impact of model parameters on Monte Carlo simulations of backscattering Mueller matrix images of colon tissue

Polarimetric imaging is emerging as a viable technique for tumor detection and staging. As a preliminary step towards a thorough understanding of the observed contrasts, we present a set of numerical Monte Carlo simulations of the polarimetric response of multilayer structures representing colon samples in the backscattering geometry. In a first instance, a typical colon sample was modeled as one or two scattering "slabs" with monodisperse non absorbing scatterers representing the most superficial tissue layers (the mucosa and submucosa), above a totally depolarizing Lambertian lumping the contributions of the deeper layers (muscularis and pericolic tissue). The model parameters were the number of layers, their thicknesses and morphology, the sizes and concentrations of the scatterers, the optical index contrast between the scatterers and the surrounding medium, and the Lambertian albedo. With quite similar results for single and double layer structures, this model does not reproduce the experimentally observed stability of the relative magnitudes of the depolarizing powers for incident linear and circular polarizations. This issue was solved by considering bimodal populations including large and small scatterers in a single layer above the Lambertian, a result which shows the importance of taking into account the various types of scatterers (nuclei, collagen fibers and organelles) in the same model.

Models of light propagation in human tissue applied to cancer diagnostics

Optical methods such as reflectance and fluorescence spectroscopy are being investigated for their potential to aid cancer detection in a quantitative, minimally invasive manner. Mathematical models of reflectance and fluorescence provide an important link between measured optical data and biomedically-relevant tissue parameters that can be extracted from these data to characterize the presence or absence of disease. The most commonly-used mathematical models in biomedical optics are the diffusion approximation (DA) to the radiative transfer equation, Monte Carlo (MC) computational models of light transport, and semi-empirical models. This paper presents a review of the applications of these models to reflectance and endogenous fluorescence sensing for cancer diagnostics in human tissues. Specific examples are given for cervical, breast, and pancreatic tissues. A comparison of the DA and MC methods in two biologically-relevant regimes of optical parameter space will also be discussed.

Fast perturbation Monte Carlo method for photon migration in heterogeneous turbid media

We present a two-step Monte Carlo (MC) method that is used to solve the radiative transfer equation in heterogeneous turbid media. The method exploits the one-to-one correspondence between the seed value of a random number generator and the sequence of random numbers. In the first step, a full MC simulation is run for the initial distribution of the optical properties and the "good" seeds (the ones leading to detected photons) are stored in an array. In the second step, we run a new MC simulation with only the good seeds stored in the first step, i.e., we propagate only detected photons. The effect of a change in the optical properties is calculated in a short time by using two scaling relationships. By this method we can increase the speed of a simulation up to a factor of 1300 in typical situations found in near-IR tissue spectroscopy and diffuse optical tomography, with a minimal requirement for hard disk space. Potential applications of this method for imaging of turbid media and the inverse problem are discussed.

Measurement of optical scattering properties with low-coherence enhanced backscattering spectroscopy

Low-coherence enhanced backscattering (LEBS) is a depth selective technique that allows noninvasive characterization of turbid media such as biological tissue. LEBS provides a spectral measurement of the tissue reflectance distribution as a function of distance between incident and reflected ray pairs through the use of partial spatial coherence broadband illumination. We present LEBS as a new depth-selective technique to measure optical properties of tissue in situ. Because LEBS enables measurements of reflectance due to initial scattering events, LEBS is sensitive to the shape of the phase function in addition to the reduced scattering coefficient (μ(s) (*)). We introduce a simulation of LEBS that implements a two parameter phase function based on the Whittle-Matérn refractive index correlation function model. We show that the LEBS enhancement factor (E) primarily depends on μ(s) (*), the normalized spectral dependence of E (S(n)) depends on one of the two parameters of the phase function that also defines the functional type of the refractive index correlation function (m), and the LEBS peak width depends on both the anisotropy factor (g) and m. Three inverse models for calculating these optical properties are described and the calculations are validated with an experimental measurement from a tissue phantom.

Quantification of the optical properties of two-layered turbid media by simultaneously analyzing the spectral and spatial information of steady-state diffuse reflectance spectroscopy

We applied hyperspectral imaging to measure spatially-resolved diffuse reflectance spectra in the visible range and an iterative inversion method based on forward Monte Carlo modeling to quantify optical properties of two-layered tissue models. We validated the inversion method using spectra experimentally measured from liquid tissue mimicking phantoms with known optical properties. Results of fitting simulated data showed that simultaneously considering the spatial and spectral information in the inversion process improves the accuracies of estimating the optical properties and the top layer thickness in comparison to methods fitting reflectance spectra measured with a single source-detector separation or fitting spatially-resolved reflectance at a single wavelength. Further development of the method could improve noninvasive assessment of physiological status and pathological conditions of stratified squamous epithelium and superficial stroma.

Quantification of the optical properties of two-layered turbid media by simultaneously analyzing the spectral and spatial information of steady-state diffuse reflectance spectroscopy

We applied hyperspectral imaging to measure spatially-resolved diffuse reflectance spectra in the visible range and an iterative inversion method based on forward Monte Carlo modeling to quantify optical properties of two-layered tissue models. We validated the inversion method using spectra experimentally measured from liquid tissue mimicking phantoms with known optical properties. Results of fitting simulated data showed that simultaneously considering the spatial and spectral information in the inversion process improves the accuracies of estimating the optical properties and the top layer thickness in comparison to methods fitting reflectance spectra measured with a single source-detector separation or fitting spatially-resolved reflectance at a single wavelength. Further development of the method could improve noninvasive assessment of physiological status and pathological conditions of stratified squamous epithelium and superficial stroma.

A robust Monte Carlo model for the extraction of biological absorption and scattering in vivo

We have a toolbox to quantify tissue optical properties that is composed of specialized fiberoptic probes for UV-visible diffuse reflectance spectroscopy and a fast, scalable inverse Monte Carlo (MC) model. In this paper, we assess the robustness of the toolbox for quantifying physiologically relevant parameters from turbid tissue-like media. In particular, we consider the effects of using different instruments, fiberoptic probes, and instrument-specific settings for a wide range of optical properties. Additionally, we test the quantitative accuracy of the inverse MC model for extracting the biologically relevant parameters of hemoglobin saturation and total hemoglobin concentration. We also test the effect of double-absorber phantoms (hemoglobin and crocin to model the absorption of hemoglobin and beta carotene, respectively, in the breast) for a range of absorption and scattering properties. We include an assessment on which reference phantom serves as the best calibration standard to enable accurate extraction of the absorption and scattering properties of the target sample. We found the best reference-target phantom combinations to be ones with similar scattering levels. The results from these phantom studies provide a set of guidelines for extracting optical parameters from clinical studies.

Higher-order perturbation theory for the diffusion equation in heterogeneous media: application to layered and slab geometries

We apply a previously proposed perturbation theory of the diffusion equation for studying light propagation through heterogeneous media in the presence of absorbing defects. The theory is based on the knowledge of (a) the geometric characteristics of a focal inclusion, (b) the mean optical path length inside the inclusion, and (c) the optical properties of the inclusion. The potential of this method is shown in the layered and slab geometries, where calculations are carried out up to the fourth order. The relative changes of intensity with respect to the unperturbed (heterogeneous) medium are predicted by the theory to within 10% for a wide range of contrasts dDeltamu(a) (up to dDeltamu(a) approximately 0.4-0.8), where d is the effective diameter of the defect and Deltamu(a) the absorption contrast between defect and local background. We also show how the method of Padé approximants can be used to extend the validity of the theory for a larger range of absorption contrasts. Finally, we study the possibility of using the proposed method for calculating the effect of a colocalized scattering and absorbing perturbation.

Monte Carlo modeling of multilayer phantoms with multiple fluorophores: simulation algorithm and experimental validation

This work is first a description of a statistical simulation algorithm developed for simulating the spectral absorption and emission of several fluorophores in an absorbing and diffusing multilayer model. Second, a detailed experimental validation of the simulation program is conducted on two sets of liquid and solid multilayer phantoms, containing one, two, or three fluorophores, within absorbing and scattering media. Experimental spatially resolved reflectance spectra are acquired in the wavelength band 400 to 800 nm and compared to corresponding simulated spectra. The degree of similarity between experimentation and simulation data is quantified. The results obtained underline good correlations with mean errors varying from 2 to 10%, depending on the number of layers and on the complexity of the phantom's composition.

Parallel computing with graphics processing units for high-speed Monte Carlo simulation of photon migration

General-purpose computing on graphics processing units (GPGPU) is shown to dramatically increase the speed of Monte Carlo simulations of photon migration. In a standard simulation of time-resolved photon migration in a semi-infinite geometry, the proposed methodology executed on a low-cost graphics processing unit (GPU) is a factor 1000 faster than simulation performed on a single standard processor. In addition, we address important technical aspects of GPU-based simulations of photon migration. The technique is expected to become a standard method in Monte Carlo simulations of photon migration.

Investigation of diffuse correlation spectroscopy in multi-layered media including the human head

In this work, diffuse correlation spectroscopy (DCS) is explored in multi-layered geometries. A quantitative comparison of an homogeneous versus a two-layered model efficiencies to recover flow changes is presented. By simulating a realistic human head with MRI anatomical data, we show that the two-layered model allows distinction between changes in superficial layers and brain. We also show that the two-layered model provides a better estimate of the flow change than the homogeneous one. Experimental measurements with a two-layered dynamical phantom confirm the ability of the two-layered analytical model to distinguish flow increase in each layer.

White Monte Carlo for time-resolved photon migration

A novel scheme for fully scalable White Monte Carlo (WMC) has been developed and is used as a forward solver in the evaluation of experimental time-resolved spectroscopy. Previously reported scaling problems are avoided by storing detection events individually, turning spatial and temporal binning into post-simulation activities. The approach is suitable for modeling of both interstitial and noninvasive settings (i.e., infinite and semi-infinite geometries). Motivated by an interest in in vivo optical properties of human prostate tissue, we utilize WMC to explore the low albedo regime of time-domain photon migration--a regime where the diffusion approximation of radiative transport theory breaks down, leading to the risk of overestimating both reduced scattering (mu(s)') and absorption (mu(a)). Experimental work supports our findings and establishes the advantages of Monte Carlo-based evaluation.

Diagnosis of breast cancer using fluorescence and diffuse reflectance spectroscopy: a Monte-Carlo-model-based approach

We explore the use of Monte-Carlo-model-based approaches for the analysis of fluorescence and diffuse reflectance spectra measured ex vivo from breast tissues. These models are used to extract the absorption, scattering, and fluorescence properties of malignant and nonmalignant tissues and to diagnose breast cancer based on these intrinsic tissue properties. Absorption and scattering properties, including beta-carotene concentration, total hemoglobin concentration, hemoglobin saturation, and the mean reduced scattering coefficient are derived from diffuse reflectance spectra using a previously developed Monte Carlo model of diffuse reflectance. A Monte Carlo model of fluorescence described in an earlier manuscript was employed to retrieve the intrinsic fluorescence spectra. The intrinsic fluorescence spectra were decomposed into several contributing components, which we attribute to endogenous fluorophores that may present in breast tissues including collagen, NADH, and retinol/vitamin A. The model-based approaches removes any dependency on the instrument and probe geometry. The relative fluorescence contributions of individual fluorescing components, as well as beta-carotene concentration, hemoglobin saturation, and the mean reduced scattering coefficient display statistically significant differences between malignant and adipose breast tissues. The hemoglobin saturation and the reduced scattering coefficient display statistically significant differences between malignant and fibrous/benign breast tissues. A linear support vector machine classification using (1) fluorescence properties alone, (2) absorption and scattering properties alone, and (3) the combination of all tissue properties achieves comparable classification accuracies of 81 to 84% in sensitivity and 75 to 89% in specificity for discriminating malignant from nonmalignant breast tissues, suggesting each set of tissue properties are diagnostically useful for the discrimination of breast malignancy.

Monte-Carlo-based model for the extraction of intrinsic fluorescence from turbid media

A Monte-Carlo-based model of fluorescence is developed that is capable of extracting the intrinsic fluorescence properties of tissue, which are independent of the absorption and scattering properties of tissue. This model is flexible in its applicability to different illumination-collection geometries and is also valid for a wide range of optical properties, representative of tissue in the UV-visible spectrum. This is potentially useful in a variety of biomedical applications, including cancer diagnostics and monitoring the physiological response to therapy. The model is validated using phantoms composed of hemoglobin (absorber), polystyrene spheres (scatterer), and furan-2 (fluorophore). It is found that this model is able to retrieve the intrinsic fluorescence spectra of the tissue phantoms and recover the intrinsic fluorescence intensity of furan within the phantoms to within a mean error of less than 10%.

Non-invasive diagnosing malignant melanoma by multi-spectral optical Nevoscope

Malignant melanoma is one of the most fatal forms of skin cancer. There are two significant signs indicating malignancy of a melanoma, that is, abnormal melanin distribution in dermis layer and a peripheral blood net. We developed a multi-spectral optical Nevoscope aimed to diagnosing malignant melanoma non-invasively. An algorithm is proposed to reconstruct the melanoma in terms of Nevoscope geometry. The algorithm has been verified on an optical tumor model at 580 nm and 800 nm. The reconstructed melanoma is consistent with the tumor model which suggests a great potential of using Nevoscope to investigate malignant melanoma.

Scaling method for fast Monte Carlo simulation of diffuse reflectance spectra from multilayered turbid media

A scaling Monte Carlo method has been developed to calculate diffuse reflectance from multilayered media with a wide range of optical properties in the ultraviolet-visible wavelength range. This multilayered scaling method employs the photon trajectory information generated from a single baseline Monte Carlo simulation of a homogeneous medium to scale the exit distance and exit weight of photons for a new set of optical properties in the multilayered medium. The scaling method is particularly suited to simulating diffuse reflectance spectra or creating a Monte Carlo database to extract optical properties of layered media, both of which are demonstrated in this paper. Particularly, it was found that the root-mean-square error (RMSE) between scaled diffuse reflectance, for which the anisotropy factor and refractive index in the baseline simulation were, respectively, 0.9 and 1.338, and independently simulated diffuse reflectance was less than or equal to 5% for source-detector separations from 200 to 1500 microm when the anisotropy factor of the top layer in a two-layered epithelial tissue model was varied from 0.8 to 0.99; in contrast, the RMSE was always less than 5% for all separations (from 0 to 1500 microm) when the anisotropy factor of the bottom layer was varied from 0.7 to 0.99. When the refractive index of either layer in the two-layered tissue model was varied from 1.3 to 1.4, the RMSE was less than 10%. The scaling method can reduce computation time by more than 2 orders of magnitude compared with independent Monte Carlo simulations.

Modeling granular phosphor screens by Monte Carlo methods

The intrinsic phosphor properties are of significant importance for the performance of phosphor screens used in medical imaging systems. In previous analytical-theoretical and Monte Carlo studies on granular phosphor materials, values of optical properties, and light interaction cross sections were found by fitting to experimental data. These values were then employed for the assessment of phosphor screen imaging performance. However, it was found that, depending on the experimental technique and fitting methodology, the optical parameters of a specific phosphor material varied within a wide range of values, i.e., variations of light scattering with respect to light absorption coefficients were often observed for the same phosphor material. In this study, x-ray and light transport within granular phosphor materials was studied by developing a computational model using Monte Carlo methods. The model was based on the intrinsic physical characteristics of the phosphor. Input values required to feed the model can be easily obtained from tabulated data. The complex refractive index was introduced and microscopic probabilities for light interactions were produced, using Mie scattering theory. Model validation was carried out by comparing model results on x-ray and light parameters (x-ray absorption, statistical fluctuations in the x-ray to light conversion process, number of emitted light photons, output light spatial distribution) with previous published experimental data on Gd2O2S: Tb phosphor material (Kodak Min-R screen). Results showed the dependence of the modulation transfer function (MTF) on phosphor grain size and material packing density. It was predicted that granular Gd2O2S: Tb screens of high packing density and small grain size may exhibit considerably better resolution and light emission properties than the conventional Gd2O2S: Tb screens, under similar conditions (x-ray incident energy, screen thickness).

Optimization of the Monte Carlo code for modeling of photon migration in tissue

The Monte Carlo method is frequently used to simulate light transport in turbid media because of its simplicity and flexibility, allowing to analyze complicated geometrical structures. Monte Carlo simulations are, however, time consuming because of the necessity to track the paths of individual photons. The time consuming computation is mainly associated with the calculation of the logarithmic and trigonometric functions as well as the generation of pseudo-random numbers. In this paper, the Monte Carlo algorithm was developed and optimized, by approximation of the logarithmic and trigonometric functions. The approximations were based on polynomial and rational functions, and the errors of these approximations are less than 1% of the values of the original functions. The proposed algorithm was verified by simulations of the time-resolved reflectance at several source-detector separations. The results of the calculation using the approximated algorithm were compared with those of the Monte Carlo simulations obtained with an exact computation of the logarithm and trigonometric functions as well as with the solution of the diffusion equation. The errors of the moments of the simulated distributions of times of flight of photons (total number of photons, mean time of flight and variance) are less than 2% for a range of optical properties, typical of living tissues. The proposed approximated algorithm allows to speed up the Monte Carlo simulations by a factor of 4. The developed code can be used on parallel machines, allowing for further acceleration.

Perturbation theory for the diffusion equation by use of the moments of the generalized temporal point-spread function. I. Theory

We approach the perturbative solution to the diffusion equation for the case of absorbing inclusions embedded in a heterogeneous scattering medium by using general properties of the radiative transfer equation and the solution of the Fredholm equation of the second kind given by the Neumann series. The terms of the Neumann series are used to obtain the expression of the moments of the generalized temporal point-spread function derived in transport theory. The moments are calculated independently by using Monte Carlo simulations for validation of the theory. While the mixed moments are correctly derived from the theory by using the solution of the diffusion equation in the geometry of interest, in order to obtain the self moments we should reframe the problem in transport theory and use a suitable solution of the radiative transfer equation for the calculation of the multiple integrals of the corresponding Neumann series. Since the rigorous theory leads to impractical formulas, in order to simplify and speed up the calculation of the self moments, we propose a heuristic method based on the calculation of only a single integral and some scaling parameters. We also propose simple quadrature rules for the calculation of the mixed moments for speeding up the computation of perturbations due to multiple defects. The theory can be developed in the continuous-wave domain, the time domain, and the frequency domain. In a companion paper [J. Opt. Soc. Am. A23, 2119 (2006)] we discuss the conditions of applicability of the theory in practical cases found in diffuse optical imaging of biological tissues.