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

Modeling voxel-based Monte Carlo light transport with curved and oblique boundary surfaces

SIGNIFICANCE

Monte Carlo (MC) light transport simulations are most often performed in regularly spaced three-dimensional voxels, a type of data representation that naturally struggles to represent boundary surfaces with curvature and oblique angles. Not accounting properly for such boundaries with an index of refractivity, mismatches can lead to important inaccuracies, not only in the calculated angles of reflection and transmission but also in the amount of light that transmits through or reflects from these mismatched boundary surfaces.

AIM

A new MC light transport algorithm is introduced to deal with curvature and oblique angles of incidence when simulated photons encounter mismatched boundary surfaces.

APPROACH

The core of the proposed algorithm applies the efficient preprocessing step of calculating a gradient map of the mismatched boundaries, a smoothing step on this calculated 3D vector field to remove surface roughness due to discretization and an interpolation scheme to improve the handling of curvature.

RESULTS

Through simulations of light hitting the side of a sphere and going through a lens, the agreement of this approach with analytical solutions is shown to be strong.

CONCLUSIONS

The MC method introduced here has the advantage of requiring only slight implementation changes from the current state-of-the-art to accurately simulate mismatched boundaries and readily exploit the acceleration of general-purpose graphics processing units. A code implementation, mcxyzn, is made available and maintained at https://omlc.org/software/mc/mcxyzn/.

Scatter Correction Based on GPU-Accelerated Full Monte Carlo Simulation for Brain PET/MRI

Accurate scatter correction is essential for qualitative and quantitative PET imaging. Until now, scatter correction based on Monte Carlo simulation (MCS) has been recognized as the most accurate method of scatter correction for PET. However, the major disadvantage of MCS is its long computational time, which makes it unfeasible for clinical usage. Meanwhile, single scatter simulation (SSS) is the most widely used method for scatter correction. Nevertheless, SSS has the disadvantage of limited robustness for dynamic measurements and for the measurement of large objects. In this work, a newly developed implementation of MCS using graphics processing unit (GPU) acceleration is employed, allowing full MCS-based scatter correction in clinical 3D brain PET imaging. Starting from the generation of annihilation photons to their detection in the simulated PET scanner, all relevant physical interactions and transport phenomena of the photons were simulated on GPUs. This resulted in an expected distribution of scattered events, which was subsequently used to correct the measured emission data. The accuracy of the approach was validated with simulations using GATE (Geant4 Application for Tomography Emission), and its performance was compared to SSS. The comparison of the computation time between a GPU and a single-threaded CPU showed an acceleration factor of 776 for a voxelized brain phantom study. The speedup of the MCS implemented on the GPU represents a major step toward the application of the more accurate MCS-based scatter correction for PET imaging in clinical routine.

Broadband Time Domain Diffuse Optical Reflectance Spectroscopy: A Review of Systems, Methods, and Applications

This review presents recent developments and a wide overview of broadband time domain diffuse optical spectroscopy (TD-DOS). Various topics including physics of photon migration, advanced instrumentation, methods of analysis, applications covering multiple domains (tissue chromophore, in vivo studies, food, wood, pharmaceutical industry) are elaborated. The key role of standardization and recent studies in that direction are discussed. Towards the end, a brief outlook is presented on the current status and future trends in broadband TD-DOS.

Physiological model using diffuse reflectance spectroscopy for nonmelanoma skin cancer diagnosis

Diffuse reflectance spectroscopy (DRS) is a noninvasive, fast, and low-cost technology with potential to assist cancer diagnosis. The goal of this study was to test the capability of our physiological model, a computational Monte Carlo lookup table inverse model, for nonmelanoma skin cancer diagnosis. We applied this model on a clinical DRS dataset to extract scattering parameters, blood volume fraction, oxygen saturation and vessel radius. We found that the model was able to capture physiological information relevant to skin cancer. We used the extracted parameters to classify (basal cell carcinoma [BCC], squamous cell carcinoma [SCC]) vs actinic keratosis (AK) and (BCC, SCC, AK) vs normal. The area under the receiver operating characteristic curve achieved by the classifiers trained on the parameters extracted using the physiological model is comparable to that of classifiers trained on features extracted via Principal Component Analysis. Our findings suggest that DRS can reveal physiologic characteristics of skin and this physiologic model offers greater flexibility for diagnosing skin cancer than a pure statistical analysis. Physiological parameters extracted from diffuse reflectance spectra data for nonmelanoma skin cancer diagnosis.

FullMonteCUDA: a fast, flexible, and accurate GPU-accelerated Monte Carlo simulator for light propagation in turbid media

Optimizing light delivery for photodynamic therapy, quantifying tissue optical properties or reconstructing 3D distributions of sources in bioluminescence imaging and absorbers in diffuse optical imaging all involve solving an inverse problem. This can require thousands of forward light propagation simulations to determine the parameters to optimize treatment, image tissue or quantify tissue optical properties, which is time-consuming and computationally expensive. Addressing this problem requires a light propagation simulator that produces results quickly given modelling parameters. In previous work, we developed FullMonteSW: currently the fastest, tetrahedral-mesh, Monte Carlo light propagation simulator written in software. Additional software optimizations showed diminishing performance improvements, so we investigated hardware acceleration methods. This work focuses on FullMonteCUDA: a GPU-accelerated version of FullMonteSW which targets NVIDIA GPUs. FullMonteCUDA has been validated across several benchmark models and, through various GPU-specific optimizations, achieves a 288-936x speedup over the single-threaded, non-vectorized version of FullMonteSW and a 4-13x speedup over the highly optimized, hand-vectorized and multi-threaded version. The increase in performance allows inverse problems to be solved more efficiently and effectively.

Head model based on the shape of the subject's head for optical brain imaging

Optical imaging methods such as near-infrared spectroscopy and diffuse optical tomography rely on models to solve the inverse problem. Imaging an adult human head also requires a head model. Using a model, which makes describing the structure of the head better, leads to acquiring a more accurate absorption map. Here, by combining the key features of layered slab models and head atlases, we introduce a new two-layered head model that is based on the surface geometry of the subject's head with variable thickness of the superficial layer. Using the Monte Carlo approach, we assess the performance of our model for fitting the optical properties from simulated time-resolved data of the adult head in a null distance source-detector configuration. Using our model, we observed improved results at 70 percent of the locations on the head and an overall 20 percent reduction in relative error compared to layered slab model.

Validation of an Inverse Fitting Method of Diffuse Reflectance Spectroscopy to Quantify Multi-Layered Skin Optical Properties

Skin consists of epidermis and dermis layers that have distinct optical properties. The quantification of skin optical properties is commonly achieved by modeling photon propagation in tissue using Monte Carlo (MC) simulations and iteratively fitting experimentally measured diffuse reflectance spectra. In order to speed up the inverse fitting process, time-consuming MC simulations have been replaced by artificial neural networks to quickly calculate reflectance spectra given tissue geometric and optical parameters. In this study the skin was modeled to consist of three layers and different scattering properties of the layers were considered. A new inverse fitting procedure was proposed to improve the extraction of chromophore-related information in the skin, including the hemoglobin concentration, oxygen saturation and melanin absorption. The performance of the new inverse fitting procedure was evaluated on 40 sets of simulated spectra. The results showed that the fitting procedure without knowing the epidermis thickness extracted chromophore information with accuracy similar to or better than fitting with known epidermis thickness, which is advantageous for practical applications due to simpler and more cost-effective instruments. In addition, the melanin volume fraction multiplied by the thickness of the melanin-containing epidermis layer was estimated more accurately than the melanin volume fraction itself. This product has the potential to provide a quantitative indicator of melanin absorption in the skin. In-vivo cuff occlusion experiments were conducted and skin optical properties extracted from the experiments were comparable to the results of previously reported in vivo studies. The results of the current study demonstrated the applicability of the proposed method to quantify the optical properties related to major chromophores in the skin, as well as scattering coefficients of the dermis. Therefore, it has the potential to be a useful tool for quantifying skin optical properties in vivo.

Contrast improvement in indocyanine green fluorescence sensing in thick tissue using a time-gating method

Indocyanine green-based fluorescence imaging techniques are very powerful in clinical applications, but the imaging is restricted to the signal from the near-surface region of tissue. Here, we focus on the method to discriminate the fluorescence signal from the background using a time-domain gating technique. The contrast of the fluorescence image from a fluorescence object at more than 1 cm depth in a meat phantom could be enhanced about 4-5 times relative to the continuous wave method if the time-gate range was properly selected. Further, a Monte Carlo simulation with a simple background model indicates that a shorter source and detector distance is more effective to improve the contrast. The simple time-gating method will enable a highly sensitive fluorescence detection in thick tissue.

Diffuse reflectance spectroscopy-enhanced drill for bone boundary detection

Intramedullary nailing is a routine orthopedic procedure used for treating fractures of femoral or tibial shafts. A critical part of this procedure involves the drilling of pilot holes in both ends of the bone for the placement of the screws that will secure the IM rod to sections of the fractured bone. This step introduces a risk of soft tissue damage because the drill bit, if not stopped in time, can transverse the bone-tissue boundary into the overlying muscle, causing unnecessary injury and prolonging healing time due to periosteum damage. In this respect, detecting the bone-tissue boundary before break-through can reduce the risks and complications associated with intramedullary nailing. Hence, in the present study, a two-wavelength diffuse reflectance spectroscopy technique was integrated into a surgical drill to optically detect bone-tissue boundary and automatically trigger the drill to stop. Furthermore, Monte-Carlo simulations were used to estimate the maximum distance from within the bone at which the bone-tissue boundary could be detected using DRS. The simulation results estimated that the detection distance, termed the "look-ahead-distance" was ∼1.5 mm for 1.3 mm source-detector fiber separation. Experimental measurements with 1.3 mm source-detector fiber separation showed that the look-ahead-distance was in the order of 250 µm in experiments with set drill rate and in the range of 1 mm in experiments where the holes were drilled by hand. Despite this difference, the automated DRS enhanced drill successfully detected the approaching bone tissue boundary when tested on samples of bovine femur and muscle tissue.

Fast and precise image generation of blood vessels embedded in skin

A software for fast rendering the visual appearance of a blood vessel located in human skin was developed based on a numerical solution of the radiative transfer equation. The user can specify geometrical properties, such as the depth and the diameter of the vessel, and physiological properties, such as the oxygen saturation of the vessel or the blood concentration in the skin. From these data, the spatially and spectrally resolved reflectance from the skin containing the blood vessel is calculated via Monte Carlo simulations, by which a two-dimensional image is generated. The short calculation time of about a second is achieved by precalculating and storing the spatially resolved reflectance for a variety of combinations of the optical and geometrical properties. This concept gives the user the opportunity to rapidly explore the influence of the physiological and geometrical properties of the investigated blood vessel on its visual appearance. The correctness of the lookup table was validated by comparison with independent Monte Carlo simulations. Rendering examples of different blood vessels in human skins are given. The current version of the software can be downloaded at https://www.ilm-ulm.de/software.

MCmatlab: an open-source, user-friendly, MATLAB-integrated three-dimensional Monte Carlo light transport solver with heat diffusion and tissue damage

While there exist many Monte Carlo (MC) programs for solving the radiative transfer equation (RTE) in biological tissues, we have identified a need for an open-source MC program that is sufficiently user-friendly for use in an education environment, in which detailed knowledge of compiling or UNIX command-line cannot be assumed. Therefore, we introduce MCmatlab, an open-source codebase thus far consisting of (a) a fast three-dimensional MC RTE solver and (b) a finite-element heat diffusion and Arrhenius-based thermal tissue damage simulator, both run in MATLAB. The kernel for both of these solvers is written in parallelized C and implemented as MATLAB MEX functions, combining the speed of C with the familiarity and versatility of MATLAB. We compare the RTE solver to Steven Jacques' mcxyz, which it is inspired by, and present example results generated by the thermal model. MCmatlab is easy to install and use and can be used by students and experienced researchers alike for simulating tissue light propagation and, optionally, thermal damage.

Graphics processing units-accelerated adaptive nonlocal means filter for denoising three-dimensional Monte Carlo photon transport simulations

The Monte Carlo (MC) method is widely recognized as the gold standard for modeling light propagation inside turbid media. Due to the stochastic nature of this method, MC simulations suffer from inherent stochastic noise. Launching large numbers of photons can reduce noise but results in significantly greater computation times, even with graphics processing units (GPU)-based acceleration. We develop a GPU-accelerated adaptive nonlocal means (ANLM) filter to denoise MC simulation outputs. This filter can effectively suppress the spatially varying stochastic noise present in low-photon MC simulations and improve the image signal-to-noise ratio (SNR) by over 5 dB. This is equivalent to the SNR improvement of running nearly 3.5  ×   more photons. We validate this denoising approach using both homogeneous and heterogeneous domains at various photon counts. The ability to preserve rapid optical fluence changes is also demonstrated using domains with inclusions. We demonstrate that this GPU-ANLM filter can shorten simulation runtimes in most photon counts and domain settings even combined with our highly accelerated GPU MC simulations. We also compare this GPU-ANLM filter with the CPU version and report a threefold to fourfold speedup. The developed GPU-ANLM filter not only can enhance three-dimensional MC photon simulation results but also be a valuable tool for noise reduction in other volumetric images such as MRI and CT scans.

Perspective: Prospects of non-invasive sensing of the human brain with diffuse optical imaging

Since the initial demonstration of near-infrared spectroscopy (NIRS) for noninvasive measurements of brain perfusion and metabolism in the 1970s, and its application to functional brain studies (fNIRS) in the 1990s, the field of noninvasive optical studies of the brain has been continuously growing. Technological developments, data analysis advances, and novel areas of application keep advancing the field. In this article, we provide a view of the state of the field of cerebral NIRS, starting with a brief historical introduction and a description of the information content of the NIRS signal. We argue that NIRS and fNIRS studies should always report data of both oxy- and deoxyhemoglobin concentrations in brain tissue, as they complement each other to provide more complete functional and physiological information, and may help identify different types of confounds. One significant challenge is the assessment of absolute tissue properties, be them optical or physiological, so that relative measurements account for the vast majority of NIRS and fNIRS applications. However, even relative measurements of hemodynamics or metabolic changes face the major problem of a potential contamination from extracerebral tissue layers. Accounting for extracerebral contributions to fNIRS signals is one of the most critical barriers in the field. We present some of the approaches that were proposed to tackle this challenge in the study of cerebral hemodynamics and functional connectivity. Finally, we critically compare fNIRS and functional magnetic resonance imaging (fMRI) by relating their measurements in terms of signal and noise, and by commenting on their complementarity.

Recent methodology advances in fluorescence molecular tomography

Molecular imaging (MI) is a novel imaging discipline that has been continuously developed in recent years. It combines biochemistry, multimodal imaging, biomathematics, bioinformatics, cell & molecular physiology, biophysics, and pharmacology, and it provides a new technology platform for the early diagnosis and quantitative analysis of diseases, treatment monitoring and evaluation, and the development of comprehensive physiology. Fluorescence Molecular Tomography (FMT) is a type of optical imaging modality in MI that captures the three-dimensional distribution of fluorescence within a biological tissue generated by a specific molecule of fluorescent material within a biological tissue. Compared with other optical molecular imaging methods, FMT has the characteristics of high sensitivity, low cost, and safety and reliability. It has become the research frontier and research hotspot of optical molecular imaging technology. This paper took an overview of the recent methodology advances in FMT, mainly focused on the photon propagation model of FMT based on the radiative transfer equation (RTE), and the reconstruction problem solution consist of forward problem and inverse problem. We introduce the detailed technologies utilized in reconstruction of FMT. Finally, the challenges in FMT were discussed. This survey aims at summarizing current research hotspots in methodology of FMT, from which future research may benefit.

Modeling photoacoustic imaging with a scanning focused detector using Monte Carlo simulation of energy deposition

Photoacoustic imaging using a focused, scanning detector in combination with a pulsed light source is a common technique to visualize light-absorbing structures in biological tissue. In the acoustic resolution mode, where the imaging resolution is given by the properties of the transducer, there are various challenges related to the choice of sensors and the optimization of the illumination. These are addressed by linking a Monte Carlo simulation of energy deposition to a time-domain model of acoustic propagation and detection. In this model, the spatial and electrical impulse responses of the focused transducer are combined with a model of acoustic attenuation in a single response matrix, which is used to calculate detector signals from a volumetric distribution of absorbed energy density. Using the radial symmetry of the detector, the calculation yields a single signal in less than a second on a standard personal computer. Various simulation results are shown, comparing different illumination geometries and demonstrating spectral imaging. Finally, simulation results and experimental images of an optically characterized phantom are compared, validating the accuracy of the model. The proposed method will facilitate the design of photoacoustic imaging devices and will be used as an accurate forward model for iterative reconstruction techniques.

Simultaneous retrieval of optical and geometrical parameters of multilayered turbid media via state-estimation algorithms

In the present paper we propose an implementation of the Kalman filter algorithm, which allows simultaneous recovery of the absorption coefficient, the reduced scattering coefficient and the thicknesses of multi-layered turbid media, with the deepest layer taken as semi-infinite. The approach is validated by both Monte Carlo simulations and experiments, showing good results in structures made up of four layers. As it is a Bayesian algorithm, prior knowledge can be included to improve the accuracy of the retrieved unknowns. One of the most promising applications of this approach is the capability of real-time monitoring of living organs by near infrared spectroscopy. In particular, determination of blood perfusion in the adult head is one of the desired goals, allowing continuous control of stroke patients. This demands accurate measurement of the optical properties, especially absorption, of the head layers, from scalp to the cortex.

High-performance, robustly verified Monte Carlo simulation with FullMonte

We introduce the FullMonte tetrahedral 3-D Monte Carlo (MC) software package for simulation, visualization, and analysis of light propagation in heterogeneous turbid media including tissue. It provides the highest computational performance and richest set of input, output, and analysis facilities of any open-source tetrahedral-mesh MC light simulator. It also provides a robust framework for statistical verification. A scripting interface makes set-up of simulation runs simple, including parameter sweeps, while simultaneously providing customization options. Data formats shared with class-leading visualization tools, VTK and Paraview, facilitate interactive generation of publication-quality fluence and irradiance maps. The simulator can read and write file formats supported by other similar simulators, such as TIM-OS, MMC, COMSOL (finite-element simulations), and MCML to support comparison. Where simulator features permit, FullMonte can take a single test case, run it in multiple software packages, and load the results together for comparison. Example meshes, optical properties, set-up scripts, and output files are provided for user convenience. We demonstrate its use in several test cases, including photodynamic therapy of the brain, bioluminescence imaging (BLI) in a mouse phantom, and a comparison against MCML for layered geometries. Application domains that can benefit from use of FullMonte include photodynamic, photothermal, and photobiomodulation therapies, BLI, diffuse optical tomography, MC software development, and biophotonics education. Since MC results may be used for preclinical or even clinical experiments, a robust and rigorous verification process is essential. Being a stochastic numerical method, MC simulation has unique challenges associated with verification of output results since observed differences may be due simply to output variance or actual differences in expected output. We describe and have implemented a rigorous and statistically justified framework for comparing between simulators of the same class and for performing regression testing.

Efficient estimation of subdiffusive optical parameters in real time from spatially resolved reflectance by artificial neural networks

Subdiffusive reflectance captured at short source-detector separations provides increased sensitivity to the scattering phase function and hence allows superficial probing of the tissue ultrastructure. Consequently, estimation of subdiffusive optical parameters has been the subject of many recent studies focusing on lookup-table-based (LUT) inverse models. Since an adequate description of the subdiffusive reflectance requires additional scattering phase function related optical parameters, the LUT inverse models, which grow exponentially with the number of estimated parameters, become excessively large and computationally inefficient. Herein, we propose, to the best of our knowledge, the first artificial-neural-network-based inverse Monte Carlo model that overcomes the limitations of the LUT inverse models and thus allows efficient real-time estimation of optical parameters from subdiffusive spatially resolved reflectance. The proposed inverse model retains the accuracy, is about four orders of magnitude faster than the LUT inverse models, grows only linearly with the number of estimated optical parameters, and can be easily extended to estimate additional optical parameters.

Parallelized Monte Carlo software to efficiently simulate the light propagation in arbitrarily shaped objects and aligned scattering media

A GPU-based Monte Carlo software (MCtet) was developed to calculate the light propagation in arbitrarily shaped objects, like a human tooth, represented by a tetrahedral mesh. A unique feature of MCtet is a concept to realize different kinds of light-sources illuminating the complex-shaped surface of an object, for which no preprocessing step is needed. With this concept, it is also possible to consider photons leaving a turbid media and reentering again in case of a concave object. The correct implementation was shown by comparison with five other Monte Carlo software packages. A hundredfold acceleration compared with central processing units-based programs was found. MCtet can simulate anisotropic light propagation, e.g., by accounting for scattering at cylindrical structures. The important influence of the anisotropic light propagation, caused, e.g., by the tubules in human dentin, is shown for the transmission spectrum through a tooth. It was found that the sensitivity to a change in the oxygen saturation inside the pulp for transmission spectra is much larger if the tubules are considered. Another "light guiding" effect based on a combination of a low scattering and a high refractive index in enamel is described.

Modelling spatially-resolved diffuse reflectance spectra of a multi-layered skin model by artificial neural networks trained with Monte Carlo simulations

A robust modelling method was proposed to extract chromophore information in multi-layered skin tissue with spatially-resolved diffuse reflectance spectroscopy. Artificial neural network models trained with a pre-simulated database were first built to map geometric and optical parameters into diffuse reflectance spectra. Nine fitting parameters including chromophore concentrations and oxygen saturation were then determined by solving the inverse problem of fitting spectral measurements from three different parts of the skin. Compared to the Monte Carlo simulation accelerated by a graphics processing unit, the proposed modelling method not only reduced the computation time, but also achieved a better fitting performance.

Scalable and massively parallel Monte Carlo photon transport simulations for heterogeneous computing platforms

We present a highly scalable Monte Carlo (MC) three-dimensional photon transport simulation platform designed for heterogeneous computing systems. Through the development of a massively parallel MC algorithm using the Open Computing Language framework, this research extends our existing graphics processing unit (GPU)-accelerated MC technique to a highly scalable vendor-independent heterogeneous computing environment, achieving significantly improved performance and software portability. A number of parallel computing techniques are investigated to achieve portable performance over a wide range of computing hardware. Furthermore, multiple thread-level and device-level load-balancing strategies are developed to obtain efficient simulations using multiple central processing units and GPUs.

Toward real-time diffuse optical tomography: accelerating light propagation modeling employing parallel computing on GPU and CPU

Parameter recovery in diffuse optical tomography is a computationally expensive algorithm, especially when used for large and complex volumes, as in the case of human brain functional imaging. The modeling of light propagation, also known as the forward problem, is the computational bottleneck of the recovery algorithm, whereby the lack of a real-time solution is impeding practical and clinical applications. The objective of this work is the acceleration of the forward model, within a diffusion approximation-based finite-element modeling framework, employing parallelization to expedite the calculation of light propagation in realistic adult head models. The proposed methodology is applicable for modeling both continuous wave and frequency-domain systems with the results demonstrating a 10-fold speed increase when GPU architectures are available, while maintaining high accuracy. It is shown that, for a very high-resolution finite-element model of the adult human head with ∼600,000 nodes, consisting of heterogeneous layers, light propagation can be calculated at ∼0.25  s/excitation source.

Analytical models for time-domain diffuse correlation spectroscopy for multi-layer and heterogeneous turbid media

A novel approach for time-domain diffuse correlation spectroscopy (TD-DCS) has been recently proposed, which has the unique advantage by simultaneous measurements of optical and dynamical properties in a scattering medium. In this study, analytical models for calculating the time-resolved electric-field autocorrelation function is presented for a multi-layer turbid sample, as well as a semi-infinite medium embedded with a small dynamic heterogeneity. To verify the analytical models, we used Monte Carlo simulations, which demonstrated that the theoretical prediction for the time-resolved autocorrelation function was highly consistent with the Monte Carlo simulation, validating the proposed analytical models. Using these analytical models, we also showed that TD-DCS has a higher sensitivity compared to conventional continuous-wave (CW) DCS for detecting the deeper dynamics. The presented analytical models and simulations can be utilized for quantification of optical and dynamical properties from future TD-DCS experimental data as well as for optimization of the experimental design to achieve maximum contrast for deep tissue dynamics.

Virtually increased acceptance angle for efficient estimation of spatially resolved reflectance in the subdiffusive regime: a Monte Carlo study

Light propagation in biological tissues is frequently modeled by the Monte Carlo (MC) method, which requires processing of many photon packets to obtain adequate quality of the observed backscattered signal. The computation times further increase for detection schemes with small acceptance angles and hence small fraction of the collected backscattered photon packets. In this paper, we investigate the use of a virtually increased acceptance angle for efficient MC simulation of spatially resolved reflectance and estimation of optical properties by an inverse model. We devise a robust criterion for approximation of the maximum virtual acceptance angle and evaluate the proposed methodology for a wide range of tissue-like optical properties and various source configurations.

Advances in Monte Carlo Simulation for Light Propagation in Tissue

Monte Carlo (MC) simulation for light propagation in tissue is the gold standard for studying the light propagation in biological tissue and has been used for years. Interaction of photons with a medium is simulated based on its optical properties. New simulation geometries, tissue-light interaction methods, and recording techniques recently have been designed. Applications, such as whole mouse body simulations for fluorescence imaging, eye modeling for blood vessel imaging, skin modeling for terahertz imaging, and human head modeling for sinus imaging, have emerged. Here, we review the technical advances and recent applications of MC simulation.

Finite element simulation of light transfer in turbid media under structured illumination

The spatial-frequency domain (SFD) imaging technique allows us to estimate the optical properties of biological tissues in a wide field of view. The technique is, however, prone to error in measurement because the two crucial assumptions used for deriving the analytical solution to the diffusion approximation cannot be met perfectly in practical applications. This research mainly focused on modeling light transfer in turbid media under the normal incidence of structured illumination using the finite element method (FEM). Finite element simulations were performed for 50 simulation samples with different combinations of optical absorption and scattering coefficients for varying spatial frequencies, and the results were then compared with the analytical method and Monte Carlo simulation. Relationships between diffuse reflectance and dimensionless absorption and dimensionless scattering coefficients were investigated. The results indicated that the FEM provided reasonable results for diffuse reflectance, compared with the analytical method. Both the FEM and the analytical method overestimated the reflectance for μ_tr/f_x values of greater than 2 and underestimated the reflectance for μ_tr/f_x values of smaller than 2. Larger values of μs'/μ_a yielded better diffuse reflectance estimations than did those of smaller than 10. The reflectance increased nonlinearly with the dimensionless scattering, whereas the reflectance decreased linearly with the dimensionless absorption. It was also observed that diffuse reflectance was relatively stable and insensitive to μs' when the dimensionless scattering was larger than 50. Overall results demonstrate that the FEM is effective for modeling light transfer in turbid media and can be used to explore the effects of crucial parameters for the SFD imaging technique.

Adopting higher-order similarity relations for improved estimation of optical properties from subdiffusive reflectance

Estimation of optical properties from subdiffusive reflectance acquired at short source-detector separations is challenging due to the sensitivity to the underlying scattering phase function. In recent studies, a second-order similarity parameter γ has been increasingly used alongside the absorption and reduced scattering coefficients to account for some of the phase function variability. By using Monte Carlo simulations, we show that the influence of the scattering phase function on the subdiffusive reflectance for the biologically relevant variations can be captured sufficiently well by considering γ and a third-order similarity parameter δ. Utilizing this knowledge, we construct an inverse model that estimates the absorption and reduced scattering coefficients, γ and δ, from spatially resolved reflectance. Nearly an order of magnitude smaller errors of the estimated optical properties are obtained in comparison to the inverse model that only composes γ.

Experimental assessment of the safety and potential efficacy of high irradiance photostimulation of brain tissues

Optogenetics is widely used in fundamental neuroscience. Its potential clinical translation for brain neuromodulation requires a careful assessment of the safety and efficacy of repeated, sustained optical stimulation of large volumes of brain tissues. This study was performed in rats and not in non-human primates for ethical reasons. We studied the spatial distribution of light, potential damage, and non-physiological effects in vivo, in anesthetized rat brains, on large brain volumes, following repeated high irradiance photo-stimulation. We generated 2D irradiance and temperature increase surface maps based on recordings taken during optical stimulation using irradiance and temporal parameters representative of common optogenetics experiments. Irradiances of 100 to 600 mW/mm2 with 5 ms pulses at 20, 40, and 60 Hz were applied during 90 s. In vivo electrophysiological recordings and post-mortem histological analyses showed that high power light stimulation had no obvious phototoxic effects and did not trigger non-physiological functional activation. This study demonstrates the ability to illuminate cortical layers to a depth of several millimeters using pulsed red light without detrimental thermal damages.

Lookup table-based sampling of the phase function for Monte Carlo simulations of light propagation in turbid media

Analytical expressions for sampling the scattering angle from a phase function in Monte Carlo simulations of light propagation are available only for a limited number of phase functions. Consequently, numerical sampling methods based on tabulated values are often required instead. By using Monte Carlo simulated reflectance, we compare two existing and propose an improved numerical sampling method and show that both the number of the tabulated values and the numerical sampling method significantly influence the accuracy of the simulated reflectance. The provided results and guidelines should serve as a good starting point for conducting computationally efficient Monte Carlo simulations with numerical phase function sampling.

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.

Time-domain Raman analytical forward solvers

A set of time-domain analytical forward solvers for Raman signals detected from homogeneous diffusive media is presented. The time-domain solvers have been developed for two geometries: the parallelepiped and the finite cylinder. The potential presence of a background fluorescence emission, contaminating the Raman signal, has also been taken into account. All the solvers have been obtained as solutions of the time dependent diffusion equation. The validation of the solvers has been performed by means of comparisons with the results of "gold standard" Monte Carlo simulations. These forward solvers provide an accurate tool to explore the information content encoded in the time-resolved Raman measurements.

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

Fluorescence lifetime measurements in heterogeneous scattering medium

Fluorescence lifetime in heterogeneous multiple light scattering systems is analyzed by an algorithm without solving the diffusion or radiative transfer equations. The algorithm assumes that the optical properties of medium are constant in the excitation and emission wavelength regions. If the assumption is correct and the fluorophore is a single species, the fluorescence lifetime can be determined by a set of measurements of temporal point-spread function of the excitation light and fluorescence at two different concentrations of the fluorophore. This method is not dependent on the heterogeneity of the optical properties of the medium as well as the geometry of the excitation–detection on an arbitrary shape of the sample. The algorithm was validated by an indocyanine green fluorescence in phantom measurements and demonstrated by an in vivo measurement.

Spatially and spectrally resolved particle swarm optimization for precise optical property estimation using diffuse-reflectance spectroscopy

This paper presents a new approach to estimate optical properties (absorption and scattering coefficients µa and µs) of biological tissues from spatially-resolved spectroscopy measurements. A Particle Swarm Optimization (PSO)-based algorithm was implemented and firstly modified to deal with spatial and spectral resolutions of the data, and to solve the corresponding inverse problem. Secondly, the optimization was improved by fitting exponential decays to the two best points among all clusters of the "particles" randomly distributed all over the parameter space (µs, µa) of possible solutions. The consequent acceleration of all the groups of particles to the "best" curve leads to significant error decrease in the optical property estimation. The study analyzes the estimated optical property error as a function of the various PSO parameter combinations, and several performance criteria such as the cost-function error and the number of iterations in the algorithms proposed. The final one led to error values between ground truth and estimated values of µs and µa less than 6%.

Radiative transfer equation modeling by streamline diffusion modified continuous Galerkin method

Optical tomography has a wide range of biomedical applications. Accurate prediction of photon transport in media is critical, as it directly affects the accuracy of the reconstructions. The radiative transfer equation (RTE) is the most accurate deterministic forward model, yet it has not been widely employed in practice due to the challenges in robust and efficient numerical implementations in high dimensions. Herein, we propose a method that combines the discrete ordinate method (DOM) with a streamline diffusion modified continuous Galerkin method to numerically solve RTE. Additionally, a phase function normalization technique was employed to dramatically reduce the instability of the DOM with fewer discrete angular points. To illustrate the accuracy and robustness of our method, the computed solutions to RTE were compared with Monte Carlo (MC) simulations when two types of sources (ideal pencil beam and Gaussian beam) and multiple optical properties were tested. Results show that with standard optical properties of human tissue, photon densities obtained using RTE are, on average, around 5% of those predicted by MC simulations in the entire/deeper region. These results suggest that this implementation of the finite element method-RTE is an accurate forward model for optical tomography in human tissues.

Generalized mesh-based Monte Carlo for wide-field illumination and detection via mesh retessellation

Monte Carlo methods are commonly used as the gold standard in modeling photon transport through turbid media. With the rapid development of structured light applications, an accurate and efficient method capable of simulating arbitrary illumination patterns and complex detection schemes over large surface area is in great need. Here we report a generalized mesh-based Monte Carlo algorithm to support a variety of wide-field illumination methods, including spatial-frequency-domain imaging (SFDI) patterns and arbitrary 2-D patterns. The extended algorithm can also model wide-field detectors such as a free-space CCD camera. The significantly enhanced flexibility of source and detector modeling is achieved via a fast mesh retessellation process that combines the target domain and the source/detector space in a single tetrahedral mesh. Both simulations of complex domains and comparisons with phantom measurements are included to demonstrate the flexibility, efficiency and accuracy of the extended algorithm. Our updated open-source software is provided at http://mcx.space/mmc.

Mesoscopic Fluorescence Molecular Tomography for Evaluating Engineered Tissues

Optimization of regenerative medicine strategies includes the design of biomaterials, development of cell-seeding methods, and control of cell-biomaterial interactions within the engineered tissues. Among these steps, one paramount challenge is to non-destructively image the engineered tissues in their entirety to assess structure, function, and molecular expression. It is especially important to be able to enable cell phenotyping and monitor the distribution and migration of cells throughout the bulk scaffold. Advanced fluorescence microscopic techniques are commonly employed to perform such tasks; however, they are limited to superficial examination of tissue constructs. Therefore, the field of tissue engineering and regenerative medicine would greatly benefit from the development of molecular imaging techniques which are capable of non-destructive imaging of three-dimensional cellular distribution and maturation within a tissue-engineered scaffold beyond the limited depth of current microscopic techniques. In this review, we focus on an emerging depth-resolved optical mesoscopic imaging technique, termed laminar optical tomography (LOT) or mesoscopic fluorescence molecular tomography (MFMT), which enables longitudinal imaging of cellular distribution in thick tissue engineering constructs at depths of a few millimeters and with relatively high resolution. The physical principle, image formation, and instrumentation of LOT/MFMT systems are introduced. Representative applications in tissue engineering include imaging the distribution of human mesenchymal stem cells embedded in hydrogels, imaging of bio-printed tissues, and in vivo applications.

Fast and efficient image reconstruction for high density diffuse optical imaging of the human brain

Real-time imaging of human brain has become an important technique within neuroimaging. In this study, a fast and efficient sensitivity map generation based on Finite Element Models (FEM) is developed which utilises a reduced sensitivitys matrix taking advantage of sparsity and parallelisation processes. Time and memory efficiency of these processes are evaluated and compared with conventional method showing that for a range of mesh densities from 50000 to 320000 nodes, the required memory is reduced over tenfold and computational time fourfold allowing for near real-time image recovery.

Accelerating fDOT image reconstruction based on path-history fluorescence Monte Carlo model by using three-level parallel architecture

The excessive time required by fluorescence diffuse optical tomography (fDOT) image reconstruction based on path-history fluorescence Monte Carlo model is its primary limiting factor. Herein, we present a method that accelerates fDOT image reconstruction. We employ three-level parallel architecture including multiple nodes in cluster, multiple cores in central processing unit (CPU), and multiple streaming multiprocessors in graphics processing unit (GPU). Different GPU memories are selectively used, the data-writing time is effectively eliminated, and the data transport per iteration is minimized. Simulation experiments demonstrated that this method can utilize general-purpose computing platforms to efficiently implement and accelerate fDOT image reconstruction, thus providing a practical means of using path-history-based fluorescence Monte Carlo model for fDOT imaging.

Limitations of the commonly used simplified laterally uniform optical fiber probe-tissue interface in Monte Carlo simulations of diffuse reflectance

Light propagation models often simplify the interface between the optical fiber probe tip and tissue to a laterally uniform boundary with mismatched refractive indices. Such simplification neglects the precise optical properties of the commonly used probe tip materials, e.g. stainless steel or black epoxy. In this paper, we investigate the limitations of the laterally uniform probe-tissue interface in Monte Carlo simulations of diffuse reflectance. In comparison to a realistic probe-tissue interface that accounts for the layout and properties of the probe tip materials, the simplified laterally uniform interface is shown to introduce significant errors into the simulated diffuse reflectance.

Near-infrared spectroscopy integrated catheter for characterization of myocardial tissues: preliminary demonstrations to radiofrequency ablation therapy for atrial fibrillation

Effects of radiofrequency ablation (RFA) treatment of atrial fibrillation can be limited by the ability to characterize the tissue in contact. Parameters obtained by conventional catheters, such as impedance and temperature can be insufficient in providing physiological information pertaining to effective treatment. In this report, we present a near-infrared spectroscopy (NIRS)-integrated catheter capable of extracting tissue optical properties. Validation experiments were first performed in tissue phantoms with known optical properties. We then apply the technique for characterization of myocardial tissues in swine and human hearts, ex vivo. Additionally, we demonstrate the recovery of critical parameters relevant to RFA therapy including contact verification, and lesion transmurality. These findings support the application of NIRS for improved guidance in RFA therapeutic interventions.

Dynamic light scattering Monte Carlo: a method for simulating time-varying dynamics for ordered motion in heterogeneous media

Few methods exist that can accurately handle dynamic light scattering in the regime between single and highly multiple scattering. We demonstrate dynamic light scattering Monte Carlo (DLS-MC), a numerical method by which the electric field autocorrelation function may be calculated for arbitrary geometries if the optical properties and particle motion are known or assumed. DLS-MC requires no assumptions regarding the number of scattering events, the final form of the autocorrelation function, or the degree of correlation between scattering events. Furthermore, the method is capable of rapidly determining the effect of particle motion changes on the autocorrelation function in heterogeneous samples. We experimentally validated the method and demonstrated that the simulations match both the expected form and the experimental results. We also demonstrate the perturbation capabilities of the method by calculating the autocorrelation function of flow in a representation of mouse microvasculature and determining the sensitivity to flow changes as a function of depth.

Mesh Optimization for Monte Carlo-Based Optical Tomography

Mesh-based Monte Carlo techniques for optical imaging allow for accurate modeling of light propagation in complex biological tissues. Recently, they have been developed within an efficient computational framework to be used as a forward model in optical tomography. However, commonly employed adaptive mesh discretization techniques have not yet been implemented for Monte Carlo based tomography. Herein, we propose a methodology to optimize the mesh discretization and analytically rescale the associated Jacobian based on the characteristics of the forward model. We demonstrate that this method maintains the accuracy of the forward model even in the case of temporal data sets while allowing for significant coarsening or refinement of the mesh.