Real-space Green's function calculation for the solution of the diffusion equation in stratified turbid media

Efficient computation of the steady-state and time-domain solutions of the photon diffusion equation in layered turbid media

Accurate and efficient forward models of photon migration in heterogeneous geometries are important for many applications of light in medicine because many biological tissues exhibit a layered structure of independent optical properties and thickness. However, closed form analytical solutions are not readily available for layered tissue-models, and often are modeled using computationally expensive numerical techniques or theoretical approximations that limit accuracy and real-time analysis. Here, we develop an open-source accurate, efficient, and stable numerical routine to solve the diffusion equation in the steady-state and time-domain for a layered cylinder tissue model with an arbitrary number of layers and specified thickness and optical coefficients. We show that the steady-state ([Formula: see text] ms) and time-domain ([Formula: see text] ms) fluence (for an 8-layer medium) can be calculated with absolute numerical errors approaching machine precision. The numerical implementation increased computation speed by 3 to 4 orders of magnitude compared to previously reported theoretical solutions in layered media. We verify our solutions asymptotically to homogeneous tissue geometries using closed form analytical solutions to assess convergence and numerical accuracy. Approximate solutions to compute the reflected intensity are presented which can decrease the computation time by an additional 2-3 orders of magnitude. We also compare our solutions for 2, 3, and 5 layered media to gold-standard Monte Carlo simulations in layered tissue models of high interest in biomedical optics (e.g. skin/fat/muscle and brain). The presented routine could enable more robust real-time data analysis tools in heterogeneous tissues that are important in many clinical applications such as functional brain imaging and diffuse optical spectroscopy.

Time-Domain Near-Infrared Spectroscopy and Imaging: A Review

This article reviews the past and current statuses of time-domain near-infrared spectroscopy (TD-NIRS) and imaging. Although time-domain technology is not yet widely employed due to its drawbacks of being cumbersome, bulky, and very expensive compared to commercial continuous wave (CW) and frequency-domain (FD) fNIRS systems, TD-NIRS has great advantages over CW and FD systems because time-resolved data measured by TD systems contain the richest information about optical properties inside measured objects. This article focuses on reviewing the theoretical background, advanced theories and methods, instruments, and studies on clinical applications for TD-NIRS including some clinical studies which used TD-NIRS systems. Major events in the development of TD-NIRS and imaging are identified and summarized in chronological tables and figures. Finally, prospects for TD-NIRS in the near future are briefly described.

Analytical reconstruction of the bioluminescent source with priors

Bioluminescence imaging has been a popular tool in small animal imaging. During the last decade, the efforts have focused on the development of tomographic systems. However, due to the difficulties in the nature of inverse source problem, multi-modal systems have been the center of attention for the last couple of years. These systems provide complementary information such that the difficulties of the inverse source problem could be overcome using the a priori information obtained. Motivated by these advances in multi-modal systems, this work presents a novel analytical reconstruction of the bioluminescent source. It is shown that if source strength is known a priori then source position could be calculated or vice versa, if source location is known a priori, source strength could be calculated as well as the photon fluence rate. The determination of the source location can be achieved by another imaging system such as X-ray computed tomography. Therefore, in bioluminescence tomography together with an imaging system can be utilized as a multi-modal system. In this work, conventional finite element based simulations are also performed and the numerical results are compared with the analytical ones. It turns out to be that the analytical results are in a good accordance with the numerical results.

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.

Application of T-matrix method in solving mixed boundary separable obstacle problem

The practical problem of imaging scatterers enclosed by separable obstacles with mixed boundary is addressed. Both the unknown scatterers and the known obstacle media can be mixture of dielectric and perfect electric conducting (PEC) materials. The scattering phenomenon of such problem is well modeled by T-matrix method. By usage of separable prior information, the obstacle media are treated as known scatterers rather than part of the background. The number of unknowns is thus reduced greatly. After recovering the profiles of scatterers by T-matrix method, a criterion is further provided to classify the PEC and dielectric scatterers. Various numerical examples are presented to show the effectiveness and good performance of the method.

Optical characterization of two-layered turbid media for non-invasive, absolute oximetry in cerebral and extracerebral tissue

We introduce a multi-distance, frequency-domain, near-infrared spectroscopy (NIRS) method to measure the optical coefficients of two-layered media and the thickness of the top layer from diffuse reflectance measurements. This method features a direct solution based on diffusion theory and an inversion procedure based on the Levenberg-Marquardt algorithm. We have validated our method through Monte Carlo simulations, experiments on tissue-like phantoms, and measurements on the forehead of three human subjects. The Monte Carlo simulations and phantom measurements have shown that, in ideal two-layered samples, our method accurately recovers the top layer thickness (L), the absorption coefficient (µ a ) and the reduced scattering coefficient (µ' s ) of both layers with deviations that are typically less than 10% for all parameters. Our method is aimed at absolute measurements of hemoglobin concentration and saturation in cerebral and extracerebral tissue of adult human subjects, where the top layer (layer 1) represents extracerebral tissue (scalp, skull, dura mater, subarachnoid space, etc.) and the bottom layer (layer 2) represents cerebral tissue. Human subject measurements have shown a significantly greater total hemoglobin concentration in cerebral tissue (82±14 µM) with respect to extracerebral tissue (30±7 µM). By contrast, there was no significant difference between the hemoglobin saturation measured in cerebral tissue (56%±10%) and extracerebral tissue (62%±6%). To our knowledge, this is the first time that an inversion procedure in the frequency domain with six unknown parameters with no other prior knowledge is used for the retrieval of the optical coefficients and top layer thickness with high accuracy on two-layered media. Our absolute measurements of cerebral hemoglobin concentration and saturation are based on the discrimination of extracerebral and cerebral tissue layers, and they can enhance the impact of NIRS for cerebral hemodynamics and oxygenation assessment both in the research arena and clinical practice.

Joint derivation method for determining optical properties based on steady-state spatially resolved diffuse reflectance measurement at small source-detector separations and large reduced albedo range: theory and simulation

Accurate determination of the optical properties (the absorption coefficient μ(a) and the reduced scattering coefficient μ(s) (')) of tissues is very important in a variety of diagnostic and therapeutic procedures. Optical diffusion theory is frequently used as the forward model for describing the photon transfer in media with large reduced albedos (a(')) and in large source-detector separations (SDS). Several other methods (PN approximation, hybrid diffusion-P3 approximation) have also been published that describe photon transfer in media with low a(') or small SDSs. We studied the theoretical models for the steady-state spatially resolved diffuse reflectance measurement to accurately determine μ(a) and μ(s) (') at large a(') range but small SDSs. Instead of using a single model, a joint derivation method is proposed. The developed method uses one of the best aforementioned theoretical methods separately in five ranges of a(') determined from several forward models. In the region of small SDSs (the range between 0.4 and 8 mm) and large a(') range (between 0.5 and 0.99), the best theoretical derivation model was determined. The results indicate that the joint derivation method can improve the derivation accuracy and that a(') range can be determined by the steady-state spatially resolved diffuse reflectance measurement.

Analytic evaluation of diffuse fluence error in multilayer scattering media with discontinuous refractive index

A simple analytic method of estimating the error involved in using an approximate boundary condition for diffuse radiation in two adjoining scattering media with differing refractive indices is presented. The method is based on asymptotic planar fluences and enables the relative error to be readily evaluated without recourse to Monte Carlo simulation. Three examples of its application are considered: (1) evaluating the error in calculating the diffuse fluences at a boundary between two media with differing refractive index and dissimilar scattering properties, (2) the dependence of the relative error in a multilayer medium with discontinuous refractive index on the ratio of the reduced scattering coefficient to the absorption coefficient μ(s)'/μ(a), and (3) the parametric dependence of the error in the radiant flux J(s) at the surface of a three-layer medium. The error is significant for strongly forward-biased scattering media with non-negligible absorption and is cumulative in multilayered media with refractive index increments between layers.

Electromagnetic imaging of separable obstacle problem

The practical problem of imaging scatterers that are separable from the known obstacles is addressed. Using such a priori information, the obstacle is regarded as a known scatterer rather than part of the background and can be excluded from the retrieving process by reformulating the cost function. As a result, the proposed method transforms the problem into an inverse scattering problem with homogeneous background, and avoids the computationally intensive calculation of Green's function for inhomogeneous background (bases of the physical model of the problem). Meanwhile, the factors that influence the imaging quality for such kind of problem are also analyzed. Various difficult numerical examples are presented to show the good performance of our method. In addition, a data set of scattering experiments from the Institut Fresnel is tested to verify the validity of our method.

Influence of boundary conditions on photon diffusion through an interface between two turbid media with different refractive indices

The photon migration in two semi-infinite highly scattering media with different refractive indices is studied in the diffusion approximation for two sets of boundary conditions at the interface. In commonly used boundary conditions, the ratio of the intensity (fluence rate) to the squared refractive index is assumed continuous across an interface and the normal component of flux is required to be continuous. However, a more rigorous approach shows that the boundary condition for the intensity may be different. As was shown by Aronson [J. Opt. Soc. Am. A12, 2532 (1995)], the ratio of the intensity to the squared refractive index undergoes a jump across an interface that is proportional to the diffuse flux. A diffusion model with an instantaneous point source that can be solved analytically for both sets of boundary conditions is considered. The analytical solutions are derived and compared with the results of Monte Carlo simulations that take into account the reflections and refractions at the interface according to Fresnel's formulas. It is shown that the analytical solutions with the Aronson boundary condition for intensity match the Monte Carlo results better than the solutions with a continuous ratio of the intensity to the squared refractive index.

Light diffusion in N-layered turbid media: steady-state domain

We deal with light diffusion in N-layered turbid media. The steady-state diffusion equation is solved for N-layered turbid media having a finite or an infinitely thick N'th layer. Different refractive indices are considered in the layers. The Fourier transform formalism is applied to derive analytical solutions of the fluence rate in Fourier space. The inverse Fourier transform is calculated using four different methods to test their performance and accuracy. Further, to avoid numerical errors, approximate formulas in Fourier space are derived. Fast solutions for calculation of the spatially resolved reflectance and transmittance from the N-layered turbid media ( approximately 10 ms) with small relative differences (<10(-7)) are found. Additionally, the solutions of the diffusion equation are compared to Monte Carlo simulations for turbid media having up to 20 layers.

Optical measurement of adipose tissue thickness and comparison with ultrasound, magnetic resonance imging, and callipers

Near-infrared spectroscopy is used to quantify the subcutaneous adipose tissue thickness (ATT) over five muscle groups (vastus medialis, vastus lateralis, gastrocnemius, ventral forearm and biceps brachii muscle) of healthy volunteers (n=20). The optical lipid signal (OLS) was obtained from the second derivative of broad band attenuation spectra and the lipid absorption peak (lambda=930 nm). Ultrasound and MR imaging as well as mechanical calliper readings were taken as reference methods. The data show that the OLS is a good predictor for ATT (<16 mm) with absolute and relative errors of <0.8 mm and <24%, respectively. The optical method compares favourably with calliper reading. The finding of a non-linear relationship of optical signal vs. ultrasound is explained by a theoretical two-layer model based on the diffusion approximation for the transport of photons. The crosstalk between the OLS and tissue hemoglobin concentration changes during an incremental cycling exercise was found to be small, indicating the robustness of OLS. Furthermore, the effect of ATT on spatially-resolved spectroscopy measurements is shown to decrease the calculated muscle hemoglobin concentration and to increase oxygen saturation.

Image reconstruction method for a two-layer tissue structure accounts for chest-wall effects in breast imaging

We develop a new tomographic imaging reconstruction algorithm for a two-layer tissue structure. Simulations and phantom experiments show more accurate reconstruction of target optical properties compared with those results obtained from a semi-infinite tissue model for layered structures. This improvement is mainly attributed to the more accurate estimation of background optical properties and more accurate estimation of weight matrix for imaging reconstruction by considering the light propagation effect in the second layer. Clinical results of breast lesions are also presented to demonstrate the utility of this new imaging algorithm.

Time-domain Green functions for diffuse light in two adjoining turbid half-spaces

Propagation of light emitted by an instantaneous source located above a plane interface between two semi-infinite turbid media is considered using the diffusion approximation. Green functions are derived for an instantaneous line source and an instantaneous point source by the method of Bellman et al. [Philos. Mag. 40, 297 (1949)], which is based on integral transforms. Both two-dimensional and three-dimensional Green functions for diffuse light have been obtained in the form of single integrals that allow for fast calculation of the specific intensity in the whole space. The influence of the optical parameters of the two media (diffusion coefficients, absorptions, and refractive indices) on the shapes of the contour lines of the specific intensity is analyzed.

Analytical solution for light propagation in a two-layer tissue structure with a tilted interface for breast imaging

Reflectance measurement of breast tissue is influenced by the underlying chest wall, which is often tilted as seen by the detection probe. We develop an analytical solution of light propagation in a two-layer tissue structure with tilted interface and refractive index difference between the layers. We validate the analytical solution with Monte Carlo simulations and phantom experiments, and a good agreement is seen. The influence of varying the tilting angle of the interface on the reflectance is discussed for two types of layered structures. Further, we apply the developed analytical solution to obtain the optical properties of breast tissue and chest wall from clinical data. Inverse calculation using the developed solution applied to the data obtained from Monte Carlo simulations shows that the optical properties of both layers are obtained with higher accuracy as compared to using a simple two-layer model ignoring the interface tilt. This is expected to improve the accuracy in estimating the optical properties of breast tissue, thus enhancing the accuracy of optical tomography of breast tumors.

Rapid simulation of steady-state spatially resolved reflectance and transmittance profiles of multilayered turbid materials

We present a technique for efficiently computing the reflection and transmission of light by arbitrary systems of turbid layers. To approximate the steady-state reflectance and transmittance without the need to solve difficult boundary conditions, we convolve the reflectance and transmittance profiles of individual layers. We extend single-slab boundary conditions to handle index-of-refraction mismatches between turbid slabs and account for interlayer scattering by applying methods similar to Kubelka-Munk theory in frequency space. We demonstrate good agreement between the reflectance and the transmittance predicted by our model and numerical Monte Carlo methods and show that the far-source reflectance and transmittance of multilayered turbid materials are dominated by interlayer scattering.

Green functions for diffuse light in a medium comprising two turbid half-spaces

A review of Green functions for diffuse light in two semi-infinite scattering and absorbing half-spaces separated by a plane interface is presented. The frequency-domain Green functions for an intensity-modulated point source are derived within the diffusion approximation by the Hankel transform with respect to the variable in the plane of the interface. Green functions for a line source and a plane source parallel to the interface are obtained from the three-dimensional Green functions by the method of descent. Green functions for a steady state are obtained as a limit of zero modulation frequency. Connection of the frequency-domain Green functions with the time-domain Green functions is shown by use of the Fourier transform in time. The influence of the relative optical parameters, namely, the ratios of diffusion coefficients, absorption coefficients, and refractive indices of the two media on the shape of the contour lines of the specific intensity, is shown for the continuous and intensity-modulated point sources.

Time-resolved multidistance near-infrared spectroscopy of the adult head: intracerebral and extracerebral absorption changes from moments of distribution of times of flight of photons

We report on multidistance time-resolved diffuse reflectance spectroscopy of the head of a healthy adult after intravenous administration of a bolus of indocyanine green. Intracerebral and extracerebral changes in absorption are deduced from moments (integral, mean time of flight, and variance) of the distributions of times of flight of photons (DTOFs), recorded simultaneously at four different source-detector separations. We calculate the sensitivity factors converting depth-dependent changes in absorption into changes of moments of DTOFs by Monte Carlo simulations by using a layered model of the head. We validate our method by analyzing moments of DTOFs simulated for the assumed changes in absorption in different layers of the head model.

Green functions for diffuse photon-density waves generated by a line source in two nonabsorbing turbid media in contact

Diffuse photon-density waves generated by an instantaneous line source that is parallel to the interface between two semi-infinite turbid media are studied by use of the diffusion approximation. For two nonabsorbing media the Green functions for diffuse light are obtained based on the Green functions for temperature fields that were derived with the Cagniard-de Hoop method. The boundary conditions for diffuse light take into account the discontinuity in the specific intensity at the interface between two media with different refractive indices. The results of the calculations of the specific intensities and the gradient lines for different sets of parameters are presented.

Asymptotic behavior and inverse problem in layered scattering media

The main challenge of noninvasive optical biopsy is to obtain an accurate value of the optical coefficients of an encapsulated organ (muscle, brain, etc.). The idea developed by us is that some interesting information could be deduced from the long-time behavior of the reflectance function. This asymptotic behavior is analyzed for layered media in the framework of the diffusion approximation. A new method is derived to obtain accurate values for the optical parameters of the deepest layers. This method is designed to work in a specific long-time regime that is still within the scope of standard time-of-flight experiments but far from being included in the mathematically defined asymptotic region. The limits of this method, linked to the cases where the asymptotic behavior is no longer governed by the deepest layer, are then discussed.

Solution of the time-dependent diffusion equation for layered diffusive media by the eigenfunction method

An exact solution of the time-dependent diffusion equation for the case of a two- and a three-layered finite diffusive medium is proposed. The method is based on the decomposition of the fluence rate in a series of eigenfunctions and upon the solution of the consequent transcendental equation for the eigenvalues obtained from the boundary conditions. Comparisons among the solution of the diffusion equation and the results of Monte Carlo simulations show the correctness of the proposed model.

Analytical approximate solutions of the time-domain diffusion equation in layered slabs

Time-domain analytical solutions of the diffusion equation for photon migration through highly scattering two- and three-layered slabs have been obtained. The effect of the refractive-index mismatch with the external medium is taken into account, and approximate boundary conditions at the interface between the diffusive layers have been considered. A Monte Carlo code for photon migration through a layered slab has also been developed. Comparisons with the results of Monte Carlo simulations showed that the analytical solutions correctly describe the mean path length followed by photons inside each diffusive layer and the shape of the temporal profile of received photons, while discrepancies are observed for the continuous-wave reflectance or transmittance.

Non-invasive determination of muscle blood flow in the extremities from laser Doppler spectra

We investigate theoretically the non-invasive determination of blood flow in muscles of the extremities using laser Doppler measurements. Laser Doppler spectra are calculated using Monte Carlo simulations and solutions of the correlation diffusion equation. The extremities are modelled as a two-layered turbid medium. The first layer represents the skin and subcutaneous fat layer and the second layer the muscle. It is shown that the absolute root-mean-square velocity of the blood in the muscle layer can be accurately derived in many practical cases if the laser Doppler spectra are measured at a distance which is sufficiently far from the source, and if the optical properties of the muscle are simultaneously determined.