1. Introduction

Non-invasive imaging of human brain function based on absorption of near-infrared light has, in the past years, seen a steady enhancement of temporal and spatial resolution (for a recent review see [1]). Using measurements of the time-of-flight of transmitted photons for multiple source-receiver pairs and sophisticated data inversion including anatomical information, full-brain maps of regional oxygen saturation in the infant brain have been measured entirely non-invasively with a spatial resolution of about 1 cm [2]. While the strong scattering of light in tissue blurs absorption contrasts measured by near-infrared spectroscopy, multiple scattered light carries rich information when the coherence length of the illumination source is large enough to produce a speckle pattern on the tissue surface. This speckle pattern contains information on the positions of all the scatterers, such as red blood cells, membranes and subcellular organelles, in the volume swept by the diffuse photon cloud. Microscopic motions of the scatterers in the tissue over distances much smaller than the wavelength of light, λ, then lead to fluctuations of the speckle pattern. This is the basis of diffusing-wave spectroscopy (DWS; also called diffuse correlation spectroscopy), the extension of quasi-elastic light scattering to multiple scattering of light [3, 4]. The central quantity measured by DWS is the temporal autocorrelation function g ⁽¹⁾(r,τ) = 〈E ^*(r,t)E(r,t+τ)〉/〈|E(r,t)|²〉 of the scattered electric field E(r,t) measured at the position r at time t, which is related to the mean-squared displacement 〈Δr ²(τ)〉 of the scatterers within the time τ by

Here k ₀ = 2π/λ is the wavenumber of light in the medium, l ^* is the transport mean free path length and P(r, s) is the normalized distribution of photon path lengths s at the position r for a source located at the origin. The latter quantity can be directly measured by a time-of-flight experiment.

Very recently, DWS has been used to detect functional activation in the human brain through intact scalp and skull, revealing a significantly accelerated decay of g ⁽¹⁾(r, τ) during contralateral stimulation of the somatomotor cortex by a finger opposition protocol [5, 6]. Analyzing measured field autocorrelation functions with a homogeneous semi-infinite 1-layer model, using an empirical correction factor accounting for the underestimation of the cortical diffusion coefficient by the 1-layer model, resulted in a functional acceleration of the DWS signal by about 40 % which was attributed to the functional enhancement of the cortical blood flow rate [5]. Modelling the field autocorrelation function by a 3-layer model reflecting the different dynamical and optical parameters in the scalp, skull and cortex, Li et al. obtained a functional increase of the cortical diffusion coefficient by about 38% [6]. These functional accelerations of the speckle fluctuations observed in DWS are consistent with the vasodilation-induced increase of the blood flow velocity during functional activation [7], and are much larger than the functional changes of blood volume or oxygen saturation measured by NIRS [8], which makes DWS a potentially attractive method to non-invasively monitor cerebral perfusion. Nevertheless, the observation that DWS autocorrelation functions measured over the human cortex are not described by directed or random flow, but rather by Brownian motion [6], has raised questions about the validity of modelling DWS data with simple single- or multilayer models for the head with fully diffusive photon transport in scalp, skull and cortex. Indeed, Monte Carlo simulations of light propagation in tissue phantoms show that the diffusion approximation for the path-length distribution P(r,s) breaks down in the presence of a extended non-scattering inclusion representing, e.g., the cerebrospinal fluid (CSF) layer surrounding the cortex [9, 10, 11]. Given that the field autocorrelation function depends, by Eq. (1), not only on the dynamics of the scatterers, but also on the path-length distribution P(r, s), the presence of the cerebrospinal fluid layer might thus affect the shape of the field autocorrelation function g ⁽¹⁾(r, τ) and the determination of parameters characterizing the cortical dynamics, such as the cortical diffusion coefficient, from experimental data.

In this publication, we study the influence of a non-scattering layer on the field autocorrelation function g ⁽¹⁾(r, τ) from head-like multilayer systems, using experimental data from a tissue phantom, Monte Carlo simulation and an extension of the analytical correlation-diffusion equation to non-diffusive boundaries in the absence of refractive index mismatches. The comparison of simulations with analytical theory on phantoms with optical parameters representative for the human head shows very good agreement provided the modified boundary conditions are taken into account. While well described by the analytical theory for large source-receiver distances, experimental data from the multilayer phantom show slight disagreement at short source-receiver distances which might arise from (i) refractive index mismatches or (ii) low-order scattering not taken into account by the diffusion theory. The analysis of simulated DWS data from a (2+1)-layer model (2 diffusive and 1 non-scattering layers) agrees very well with the analytical (2+1)-layer theory over a wide range of diffusion coefficients characterizing the cortical dynamics, allowing to retrieve cortical diffusion coefficients to within 20%. The simple 2-layer diffusion theory neglecting the non-scattering CSF layer agrees reasonably well with the simulated (2+1)-layer data, but underestimates the cortical diffusion coefficient by about 40%.

The paper is organized as follows: in Section 2.1, we derive an expression for the field autocorrelation function g ⁽¹⁾(r, τ) for the backscattering geometry with both point-like source and receiver in the presence of a non-scattering layer located between two scattering layers and briefly describe the Monte Carlo (MC) simulation procedure used to obtain g ⁽¹⁾(r, τ) for a multilayer medium. Experimental details are given in Section 2.2. The influence of the thickness and absorption coefficient of the non-scattering layer on g ⁽¹⁾(r, τ) is studied with simulation data and theory in Section 3.1, and experiments on a multilayer phantom are compared with theory in Section 3.2.

2. Materials and methods

2.1. Theory and simulations

A. Diffusion model

We will restrict ourselves to the case of the backscattering geometry in a laterally infinite medium where the point-like source and receiver are separated by the distance ρ. The diffusion approximation allows to model the propagation of light in a multilayer medium as shown in Fig. 1. A clear, non-scattering layer with thickness d is located between two scattering layers (layer 1, thickness Δ₁, and layer 2 with thickness Δ₂).

 

Fig. 1. Anatomical magnetic resonance image of human head structure (left) modelled by a laterally infinite multilayer geometry (right). CSF: cerebrospinal fluid layer. Source (S) and receiver (R) are assumed to be point-like.

Download Full Size | PDF

The diffusion approximation is not valid in a non-scattering medium such as the CSF. Nevertheless, by accounting for its presence by the boundary conditions between this non-scattering layer and the scattering layers, the diffusion equation can be recovered, as was shown for the case of the transmission geometry with plane-wave illumination [12]. Similar to the diffusion equation for the radiance, the (unnormalized) temporal autocorrelation function

satisfies the diffusion equation

where α ²(τ) = 3μ _a μ′_s + $k_0^2$ μ${\prime }_{\mathrm{s}}^2$〈Δr ²(τ)〉, with μ _a the absorption coefficient and μ _s′ = 1/l ^* the reduced scattering coefficient. The source located at r′ = {ρ′ = 0,z′} inside the first layer at depth z′ = 1/μ′_{s 1} has a time-independent strength s ₀. Eq. (3) is conveniently solved using the Fourier transform of G(r, τ) with respect to the transverse coordinate ρ

whose solution for the mth layer is [6]

where ${\beta }_m^2$(q, τ) = ${\alpha }_m^2$(τ) + q ². The constants A_m and B_m are determined by the boundary conditions. Near the surface (z = 0) the boundary conditions are [6]

where z_m is the extrapolation length of layer m [6]. Here we have divided the first layer into 2 sub-layers: layer 0 for 0 < z < z′ and layer 1 for z′ < z < Δ₁. Layer 2 (z > Δ₁ + d) is separated from layer 1 by the non-scattering layer of thickness d. Light which is scattered from layer 1 into the non-scattering layer may travel ballistically and will thus be converted to a diffusing photon in layer 2 with an angle-dependent transfer probability 0 ≤ f(q = |q|) ≤ 1. The same holds for the transfer of photons from layer 2 to layer 1. The boundary conditions including the transfer probability f(q) can thus be written as

The Fourier transform of the measured field autocorrelation function at the surface is given by

where for Δ₂ → ∞, i.e. a semi-infinite layer 2,

and

In order to determine the transfer probability f(q), we use the formalism given in [13], assuming that the sample is infinite in the direction parallel to the surface. We find (see Appendix)

where J ₀(x) is the zeroth order Bessel function of the first kind and μ _a the absorption coefficient of the non-scattering layer. In the case of a non-scattering layer without absorption, Eq. (14) simplifies to

where K ₁(x) is the modified Bessel function of the second kind. Note that for vanishing thickness of the non-scattering layer (d = 0), f(q) = 1 and the diffusive boundary conditions of [14] are recovered. The field autocorrelation function G ₀(r,τ) at the position r = {ρ, z = 0} at the surface of the sample is then obtained by the inverse Fourier transform of Ĝ ₀(q,z = 0,τ) [6],

The field autocorrelation function g ⁽¹⁾(r, τ) is then obtained by normalization of G ₀(r, τ) with its value at τ = 0.

In order to test the validity of this model, we carried out Monte Carlo simulations that allow to look at the case with vanishing refractive index mismatch between the non-scattering and the scattering layers, in which case the extrapolation lengths are given by z_m = 2/(3μ′_sm). The next section briefly reviews the basic scheme of the Monte Carlo simulation procedure.

B. Monte Carlo simulations

Monte Carlo simulations of the field autocorrelation function g ⁽¹⁾(r, τ) were performed by simulating the light propagation [15] in the multilayer geometry of Fig. 1 (with P = 3 layers and a pixel size of 0.05cm). Each layer is characterized by an absorption coefficient μ _a, a reduced scattering coefficient μ′_s, and its thickness. The refractive index of the layers was chosen to be uniform, n _med = 1.33. At each scattering event, the photon direction is randomly chosen such that the Henyey-Greenstein phase function with anisotropy factor g = 0 is sampled for all layers. Photons entering or leaving the non-scattering layer pass through the layer without being refracted or scattered. Simulations were carried out with N = 4 × 10⁶ photons.

The field autocorrelation function measured at the surface at the position r = {ρ,z = 0} can be written as

where w_α (r) is the weight of photon path α (with M_α steps) starting at the origin and ending at the position r, determined by the absorption between successive scattering events [16]. The quantity $g_{\alpha }^{\left(1\right)}$(τ) is the field autocorrelation function for photon path α which is given by

where q_αl is the magnitude of the scattering vector of the lth scattering event in path α and 〈Δ$r_j^2$(τ)〉 is the mean-square displacement within time τ in layer j. The weighting factor h_αlj is unity if the scattering event l of path α is within layer j, and zero otherwise. The motion of the scatterers is modelled by Brownian diffusion, 〈Δ$r_j^2$(τ)〉 = 6D_jτ, characterized by a diffusion coefficient D_j for each layer.

2.2. Experimental

Experiments were carried out at a temperature T = 295 K on a phantom consisting of a cylindrical vat made of stainless steel (diameter 15 cm) whose 2 fluid compartments were separated by a glass window with thickness d = 0.05cm and refractive index ≈ 1.4 (see Fig. 2). Fluid layers consisted of turbid aqueous suspensions (refractive index n _med = 1.33) of polystyrene latex spheres. Particle diffusion coefficients, as measured by quasi-elastic light scattering on dilute suspensions, were D ₁ = 9.2 × 10^-9cm²/s for layer 1 (representing scalp and skull) and D ₂ = 1.26 × 10^-8cm²/s for layer 2 (representing the cortex). Modelling the dynamics within the skull by Brownian motion is motivated by the observation that DWS data from humans probing scalp and skull show no sign of static scattering (detectable by a reduction of the intercept of the intensity autocorrelation function [17]), reflecting the fact that in both layers the perfusion is strong enough to provide a fully fluctuating speckle pattern. Thicknesses of layers 1 and 2 were Δ₁ = 0.88 cm and Δ₂ = 7 cm, respectively. Layer 2 is sufficiently thick that for the source-receiver distances used here, it can be assumed to be semi-infinite. Light from a diode laser with wavelength λ ₀ = 802 nm (Toptica TA-100) was delivered into a point-like source by a multimode fiber. Multiple scattered light was collected at distances 0.5cm ≤ ρ ≤ 4 cm from the source by a 6-mode fiber [18] and detected by an avalanche photodiode (Perkin-Elmer SPCM-AQR-15-FC). The normalized autocorrelation function g ⁽²⁾(r, τ) of the measured photon count rate was computed by a correlator (ALV-5000E). In order to eliminate artefacts, the top liquid layer was covered with a rigid black plastic sheet serving both to reduce internal reflections at the free liquid surface, and to suppress surface fluctuations. Through small holes in the cover, the source and receiver fibers were immersed about 1–2 mm into the liquid in order to ensure stable optical coupling.

Reduced scattering coefficients of layers 1 and 2 were determined from measurements of the DWS signal in bulk samples of the respective latex suspensions in glass containers. Both source and receiver fibers were immersed about 2.1 cm into the suspensions, with similar distances to the walls of the container. Measured g ⁽¹⁾(r, τ) were analyzed using the solution of the diffusion equation for the infinite geometry [19], the independently determined values of the diffusion coefficients D _1,2 and assuming the absorption coefficient μ _a1,2 = 0.0223 cm^-1 of water [20]. Measured values μ′_{s 1} = 8cm^-1 and μ′_{s 2} = 17.24cm^-1 agreed to within 4% with calculations based on Mie theory [21]. Anisotropy factors g ₁ = 0.94 and g ₂ = 0.92 for layers 1 and 2, respectively, were calculated using Mie theory [21].

 

Fig. 2. Schematic cross-section (not to scale) of the multilayer phantom. Scalp and skull are represented by layer 1 and the cortex by layer 2. The cerebrospinal fluid layer is placed between layers 1 and 2. The plastic sheet covering the top layer is used to suppress artefacts in the DWS signal. Fibers are immersed about 1–2mm into the top liquid to ensure stable optical coupling. Valve-equipped inlets allow for bubble-free filling of the lower compartment.

Download Full Size | PDF

3. Results

3.1. Comparison between simulation and theory

A comparison of solutions of the correlation-diffusion equation (3) and simulated field autocorrelation functions for the multilayer geometry of Fig. 1 is shown in Fig. 3 for source-receiver distances ρ = 10 mm and ρ = 20 mm. We note the very good agreement between theory and simulations both with and without the (non-absorbing) non-scattering layer between layers 1 and 2, provided the modified boundary conditions Eqs. (9)–(10) and the expression Eq. (15) for the transfer probability f(q) are used. Similar to the DWS signal in transmission geometry, the decay of g ⁽¹⁾(r, τ) is shifted towards shorter times as the source-receiver distance is increased, reflecting the increasing contribution of long photon paths to the DWS signal. The presence of a non-scattering layer separating the turbid layers 1 and 2 results in a slowing-down of the field autocorrelation function. This is due to the fact that when photons are scattered into the non-scattering layer, there is a finite probability that they travel ballistically far outside the banana-shaped acceptance volume spanned by source and receiver [19], thus reducing the contribution of long photon paths to the DWS signal. This is indeed reflected by the relatively weak effect of the non-scattering layer on g ⁽¹⁾(r, τ) at the shortest source-receiver distance ρ = 10mm where the typical depth probed by DWS, ρ/2 = 5 mm, is smaller than the thickness of the top layer, Δ₁ = 8.8mm.

Figure 4 shows that when for fixed source-receiver distance the thickness of the non-scattering (non-absorbing) layer is increased, the decay of g ⁽¹⁾(r, τ) is increasingly slowed down, reflecting the larger probability for the photon to escape the acceptance volume by being scattered into the non-scattering layer, which results in a reduced contribution of long paths to g ⁽¹⁾(r, τ).

Figure 5 shows the comparison between the diffusion model and the MC simulations as one changes the absorption of the non-scattering layer at constant thickness d = 0.1 cm, calculating the transfer probability f(q) by Eq. (14). As would be expected for increased absorption loss, the decay of g ⁽¹⁾(r, τ) is slowed down as the absorption coefficient of the non-scattering layer is increased. This effect is particularly pronounced for source-receiver distances large enough that the acceptance volume intersects the non-scattering layer. Absorption within the non-scattering layer thus counteracts the ballistic propagation of photons within the non-scattering layer and tends to restore diffusive behavior.

3.2. Comparison between theory and experiments

As shown in Fig. 6, the field autocorrelation functions measured from the multilayer phantom agree very well with the predictions of the correlation-diffusion equation for a (2+1)-layer medium with boundary conditions modified to account for the presence of a non-absorbing, non-scattering layer. For increasing source-receiver distance ρ, the decay of the field autocorrelation function shifts monotonically towards shorter times, as would be expected for the increasing weight of long photon paths to g ⁽¹⁾(r, τ). We find that the presence of the opaque cover and the immersion of the fibers into the top layer is essential for observing the increasingly fast decay of g ⁽¹⁾(r, τ) with increasing source-receiver distance. While the agreement between theory and experiment is excellent for ρ > 1.5 cm, it is less so at shorter distances, in particular for long times (τ > 10^-4 s). We think that the observed discrepancies arise from contributions of low-order scattering which are not correctly accounted for by the diffusion theory.

4. Discussion

As shown in Sec. 3.1, MC simulations and calculations based on the correlation-diffusion equation agree closely even in the presence of a non-scattering layer provided that the boundary conditions for the correlation-diffusion equation are corrected for the non-diffusive correlation fluxes at the boundaries with the non-scattering layer. For vanishing refractive index mismatch between the layers, we find that deviations between simulations and theory become noticeable (yet still small) at long lag times τ (corresponding to short photon path lengths s) and increasing thickness of the non-scattering layer. On the other hand, when the absorption coefficient of the non-scattering layer is finite, the agreement of theory with simulation is significantly improved, reflecting the fact that diffusive propagation is restored by suppression of long ballistic photon paths within the non-scattering layer.

 

Fig. 3. Field autocorrelation functions g ⁽¹⁾(r, τ) for a 2-layer geometry calculated from the solution of the correlation-diffusion Eq. (3) (lines) and from simulation data (symbols), at different source-receiver distances ρ = 10mm (a) and ρ = 20mm (b). Squares: g ⁽¹⁾(r, τ) in the presence of a non-absorbing, non-scattering layer with thickness d = 5 mm separating the top and bottom turbid layers; circles: g ⁽¹⁾(r, τ) without the non-scattering layer. Parameters for the calculations are μ′_{s 1} = 8.0 cm^-1, μ _a1 = 0.0223 cm^-1, thickness Δ₁ = 0.88 cm for the top layer 1, and μ′_{s 2} = 17.24 cm^-1, μ _a2 = 0.0223 cm^-1, and Δ₂ = 15cm for the bottom layer 2. Particle diffusion coefficients for layers 1 and 2 are D ₁ = 0.92 × 10^-8 cm²/s and D ₂ = 1.26 × 10^-8cm²/s. The laser wavelength is λ ₀ = 802nm, the refractive index is n _med = 1.33. For the simulations, the anisotropy factor is g = 0.

Download Full Size | PDF

Experimental data from a (2+1)-layer phantom agree well with the predictions of the correlation-diffusion theory with boundary conditions Eqs. (6)–(10). In particular the agreement between theory and the data at the large source-receiver distances ρ = 1.5cm is remarkable given that the theory curves contain no adjustable parameters. This indicates that for source-receiver distances large enough to probe the cortex, the refractive index mismatch between layers 1 and 2 and the intervening non-scattering layer (which is present in the experiment, but neglected in the theory) only marginally affects the accuracy of the diffusion theory, provided the presence of the non-scattering layer is accounted for by the transfer probability f(q). The discrepancy between theory and experiment for short source-receiver distances (ρ ≤ 1cm) and long times is due to the breakdown of the diffusion approximation for short photon paths. Enhanced agreement between experiment and theory for short photon paths could be achieved either by simulating the light transport with realistic values of the anisotropy factors, or by solving the radiative transfer equation [22], which for the present multilayer geometry is however beyond the scope of the present work.

 

Fig. 4. Influence of the thickness d of the non-scattering layer thickness on the field autocorrelation functions g ⁽¹⁾(r,τ). Symbols denote MC simulation results: d = 0.5cm (□), d = 0.3cm (■), and d = 0.1cm (◆). Lines: solutions of the correlation-diffusion equation. The source-receiver distance is ρ = 2cm. Optical and dynamical parameters are as in Fig. 3.

Download Full Size | PDF

In analyzing DWS data obtained from the human head through intact scalp and skull, the thickness of the non-scattering cerebrospinal fluid layer is, in general, difficult to estimate. The question thus arises how much the neglect of the CSF layer in the analysis of DWS data will affect the determination of cortical diffusion coefficients. We have thus simulated field autocorrelation functions over a large range of cortical diffusion coefficients 10^-9cm²/s ≤ D ₂ ≤ 10^-6cm²/s for the optical and geometrical parameters of the phantom experiment (representative of the situation in the human head). From the simulated data we retrieved the cortical diffusion coefficient D ₂ by fitting the diffusion theory to the MC data, using a Levenberg-Marquardt optimization routine. As shown in Fig. 7, the analysis of the MC data with a 2-layer model, neglecting the presence of the non-scattering layer, underestimates the cortical diffusion coefficient by about 40%. When the (2+1)-layer diffusion theory is used, the agreement with the simulation is enhanced: the error of the retrieved cortical diffusion coefficient D ₂ increases from very small values at D ₂ ≤ 5 × 10^-8 cm²/s to about 20 % at D ₂ = 10^-6 cm²/s. The underestimation of the cortical diffusion coefficient is mainly due to a discrepancy between simulation and theory at long times due to the breakdown of the diffusion approximation. Using only the short-time data with g ⁽¹⁾(r, τ) > 0.1 for the fitting, the error in the retrieved D ₂ can be made very small (data not shown). However, this cutoff is arbitrary and should be used with caution.

The sensitivity of the field autocorrelation function to the presence of a non-scattering layer for the present backscattering geometry is in marked contrast to the small effect of a non-scattering layer in transmission geometry with plane-wave illumination [12]. This is due to the fact that the contribution of long photon paths to the field autocorrelation function is reduced by the presence of the CSF, both by the transfer probability f(q) < 1 and by absorption, while the propagation of detected photons diffusing on short paths through layer 1 only is only marginally modified. In the transmission geometry, on the other hand, g ⁽¹⁾(r, τ) is dominated by photon paths close to the line of sight between source and receiver, resulting in an only minor distortion of the path length distribution function by a non-scattering layer.

 

Fig. 5. Field autocorrelation functions g ⁽¹⁾(r, τ) for different values of the absorption coefficient μ _a of the non-scattering layer with thickness d = 0.1 cm. Symbols denote MC simulation results: μ _a = 0.01 cm^-1 (◇), μ _a = 0.5cm^-1 (●), and μ _a = 1 cm^-1 (□). Lines: solutions of the correlation-diffusion equation. The source-receiver distance is ρ = 2 cm. Other optical and dynamical parameters are as in Fig. 3.

Download Full Size | PDF

5. Conclusions

We studied the influence of a non-scattering layer (mimicking the cerebrospinal fluid [CSF] surrounding the brain) in a head-like tissue phantom on the field autocorrelation function g ⁽¹⁾(r, τ) of multiple scattered light by means of a diffusion theory, Monte Carlo simulations and experiments for the backscattering geometry used in DWS experiments on the human head. For a configuration with a point source and a point receiver, boundary conditions for the correlation-diffusion equation were introduced accounting for the presence of the non-scattering layer. Field autocorrelation functions g ⁽¹⁾(r, τ) predicted by this modified diffusion theory are in very good agreement with Monte Carlo simulations in the absence of refractive index mismatches between the layers. The comparison with experimental data indicates that contributions from low-order scattering and, possibly, refractive index mismatches between the layers play a role for source-receiver distances comparable to about twice the distance between surface and the cerebrospinal fluid layer. The proper treatment of low-order scattering and refractive index mismatches might thus be important for the quantitative analysis of DWS data from the human head measured at short source-receiver distances. Data measured at larger source-receiver distance which are able to probe the cortex are, on the other hand, rather insensitive to refractive index mismatches. The analysis by the diffusion theory qualitatively reproduces the simulated DWS data for the optical and dynamical properties used here even when the presence of the non-scattering layer is neglected. Nevertheless, the coupling between superficial and deep tissue layers in the DWS signal and the ensuing complexity of the DWS signal for a multilayer system makes it difficult to chart the region in parameter space where the diffusion theory satisfyingly describes the simulations for other combinations of μ′_s, μ _a, D ₁ and D ₂. Methods for selectively probing deep tissue layers by DWS might ultimately help obviate these complications.

 

Fig. 6. Field autocorrelation function g ⁽¹⁾(r, τ) measured from a (2+1)-layer phantom with diffusive dynamics in the turbid layers 1 and 2 (symbols) and predictions of the correlation-diffusion equation with modified boundary conditions for a non-absorbing, non-scattering layer (lines), for source-receiver distances ρ = 1.0cm (◇), ρ = 1.5cm (□), and ρ = 2cm (△). Optical parameters are μ′_{s 1} = 8cm^-1 and μ′_{s 2} = 17.24cm^-1, μ _a1,2 =0.0223cm^-1, and diffusion coefficients are D ₁ = 9.2 × 10^-10cm²/s for layer 1 and D ₂ = 1.26 × 10^-9cm²/s for layer 2. Layer thicknesses are Δ₁ = 0.88cm and Δ₂ = 15 cm for layers 1 and 2, respectively, and d = 0.05 cm for the non-scattering layer.

Download Full Size | PDF

Appendix: Boundary conditions

We start from Eq. (23) in Ref [13], using the autocorrelation function G(r, τ) instead of the average intensity U(r). In our case, since we have no light sources at the boundary (the right-hand side of Eq. (24) of Ref. [13] vanishes) the total normal correlation flux is given by

(19)$${퓙}_m(\mathbf{r},\tau )=-{퓓}_m^{\prime }{\nabla }_{\mathbf{n}}G(\mathbf{r},\tau )$$

where 𝓓′_m = 3(μ _am + μ′_sm)]^-1 is the reduced photon diffusion coefficient of the mth layer, and m is the unit vector normal to the interface. Eq. (36) of Ref [13] is written with G(r, τ) as

(20)$$G(\mathbf{r},\tau )=C_m{퓙}_m(\mathbf{r},\tau )+\frac{1}{\pi }{\int }_S\left[G(\mathbf{r}\prime ,\tau )+\frac{R_J}{R_U}{퓙}_m(\mathbf{r}\prime ,\tau )\right]𝓖\left(\mathbf{r}-\mathbf{r}\prime \right)dS\prime$$

where S is the surface the irradiation originates from, which in our case is the xy plane, C_m = (2 - R_J )/R_U = 2 since R_U = ${\int }_0^1$[1 - |R _0→1(θ)|]² cosθ d(cosθ) = 1/2 and R_J = 3 ${\int }_0^1$[1 -|R _0→1(θ)|]² cos² θ d(cosθ) = 1 in the case without refractive index mismatch (R _0→1(θ), the reflection coefficient with incidence angle θ from medium 0 to medium 1, is equal to 1) from Eq. (30) of Ref [13]. Then the boundary conditions at z = Δ₁ are

(21)$$G_1\left(\mathbf{\rho },z={\Delta }_1,\tau \right)=-2{퓓}_1^{\prime }\frac{\partial G_1\left(\mathbf{\rho },z={\Delta }_1,\tau \right)}{\partial z}+$$
$$\frac{1}{\pi }{\int }_S{d}^2{\mathbf{\rho }}^{\prime }\left[G_2\left({\mathbf{\rho }}^{\prime },z={\Delta }_1+d,\tau \right)+2{퓓}_2^{\prime }\frac{\partial G_2\left({\mathbf{\rho }}^{\prime },z={\Delta }_1+d,\tau \right)}{\partial z}\right]𝓖\left(\mathbf{\rho }-{\mathbf{\rho }}^{\prime }\right)$$

and at z = Δ₁ + d:θ

(22)$$G_2\left(\mathbf{\rho },z={\Delta }_1+d,\tau \right)=2{퓓}_2^{\prime }\frac{\partial G_1\left(\mathbf{\rho },z={\Delta }_1+d,\tau \right)}{\partial z}+$$
$$\frac{1}{\pi }{\int }_S{d}^2{\mathbf{\rho }}^{\prime }\left[G_1({\mathbf{\rho }}^{\prime },z={\Delta }_1,\tau )+(-2{퓓}_1^{\prime }\frac{\partial G_1\left({\mathbf{\rho }}^{\prime },z={\Delta }_1,\tau \right)}{\partial z})\right]퓖\left(\mathbf{\rho }-{\mathbf{\rho }}^{\prime }\right)$$

The quantity 𝓰(ρ - ρ′) is given by Eq. (18) of Ref. [13],

(23)$$𝓖\left(\mathbf{\rho }-{\mathbf{\rho }}^{\prime }\right)=\frac{\exp (-{\mu }_a\mid \mathbf{\rho }-{\mathbf{\rho }}^{\prime }\mid )}{{\mid \mathbf{\rho }-{\mathbf{\rho }}^{\prime }\mid }^2}\cos \theta \cos {\theta }^{\prime }$$

where μ _a is the absorption coefficient of the non-scattering layer. Defining ρ = |ρ - ρ′| then $\cos \theta =\cos {\theta }^{\prime }=\frac{d}{\sqrt{d^2+{\rho }^2}}$. Thus

(24)$$퓖\left(\rho \right)=d^2\frac{(-{\mu }_a\sqrt{\rho +d^2})}{{\left({\rho }^2+d^2\right)}^2}.$$

Taking the Fourier transform of Equations (21) and (22) we obtain, after some calculations,

(25)$${\hat{G}}_1\left(\mathbf{q},z={\Delta }_1,\tau \right)+2{퓓}_1^{\prime }\frac{\partial {\hat{G}}_1\left(\mathbf{q},z={\Delta }_1,\tau \right)}{\partial z}=f\left(q\right)\left[{\hat{G}}_2\left(\mathbf{q},z={\Delta }_1+d,\tau \right)+2{퓓}_2^{\prime }\frac{\partial {\hat{G}}_2\left(\mathbf{q},z={\Delta }_1+d,\tau \right)}{\partial z}\right]$$

(26)$${\hat{G}}_2\left(\mathbf{q},z={\Delta }_1+d,\tau \right)-2{퓓}_2^{\prime }\frac{\partial {\hat{G}}_2\left(\mathbf{q},z={\Delta }_1+d,\tau \right)}{\partial z}=f\left(q\right)\left[{\hat{G}}_1\left(\mathbf{q},z={\Delta }_1,\tau \right)-2{퓓}_1^{\prime }\frac{\partial {\hat{G}}_1\left(\mathbf{q},z={\Delta }_1,\tau \right)}{\partial z}\right]$$

which are equivalent to Eqs. (9) and (10). The transfer probability f(q) is given by

(27)$$f\left(q\right)=\frac{1}{\pi }\int d^2\mathbf{\rho }\frac{d^2}{{\left({\rho }^2+d^2\right)}^2}\exp (-{\mu }_a\sqrt{\rho +d^2})\exp (i\mathbf{q}\bullet \mathbf{\rho })$$

which reduces to Eq. (14), and, for μ _a = 0, to Eq. (15).

 

Fig. 7. Error of the cortical diffusion coefficient D ₂ from a (2+1)-layer geometry retrieved by analysis of g ⁽¹⁾(r, τ) with a 2-layer diffusion model with and without a CSF layer (◆ and □, respectively). Data for D ₂ < 10^-8 cm²/s (not shown) differ strongly from the true values, possibly due to round-off errors introduced by the numerical Fourier transform in Eq. (16). The thickness of the CSF layer is d = 0.3cm, and the source-receiver spacing is ρ = 2cm. Other parameters are as in Fig. 3.

Download Full Size | PDF

Acknowledgments

We thank S. Eiden and her group for providing the latex samples, and B. Awiszus for the MRI recording. This work is funded by the Deutsche Forschungsgemeinschaft (DFG), the Lan-desstiftung Baden-Württemberg, the Center for Applied Photonics (CAP) Konstanz and the Zentrum für Wissenschaftlicher Nachwuchs (ZWN) at the University of Konstanz.

References and links

1. A. P. Gibson, J. C. Hebden, and S. R. Arridge, “Recent advances in diffuse optical imaging,” Phys. Med. Biol. 50, R1–R43 (2005). [CrossRef]   [PubMed]  

2. J. C. Hebden, A. P. Gibson, T. Austin, R. M. Yusof, N. Everdell, D. T. Delpy, S. R. Arridge, J. H. Meek, and J. S. Wyatt, “Imaging changes in blood volume and oxygenation in the newborn infant brain using three-dimensional optical tomography,” Phys. Med. Biol. 49, 1117–1130 (2004). [CrossRef]   [PubMed]  

3. G. Maret and P. E. Wolf, “Multiple light scattering from disordered media: The effect of Brownian motion of scatterers,” Z. Phys. B 65, 409–413 (1987). [CrossRef]  

4. D. J. Pine, D. A. Weitz, P. M. Chaikin, and E. Herbolzheimer, “Diffusing-Wave Spectroscopy,” Phys. Rev. Lett. 60, 1134–1137 (1988). [CrossRef]   [PubMed]  

5. T. Durduran, G. Yu, M. G. Burnett, J. A. Detre, J. H. Greenberg, J. Wang, C. Zhou, and A. G. Yodh, “Diffuse optical measurement of blood flow, blood oxygenation, and metabolism in a human brain during sensorimotor cortex activation,” Opt. Lett. 29, 1766–1768 (2004). [CrossRef]   [PubMed]  

6. J. Li, G. Dietsche, D. Iftime, S. E. Skipetrov, G. Maret, T. Elbert, B. Rockstroh, and T. Gisler, “Non-Invasive Detection of Functional Brain Activity with Near-Infrared Diffusing-Wave Spectroscopy,” J. Biomed. Opt. 10, 044002-1–12 (2005). [CrossRef]  

7. H. Ito, K. Takahashi, J. Hatazawa, S.-G. Kim, and I. Kanno, “Changes in Human Regional Cerebral Blood Flow and Cerebral Blood Volume During Visual Stimulation Measured by Positron Emission Tomography,” J. Cereb. Blood Flow Metab. 21, 608–612 (2001). [CrossRef]   [PubMed]  

8. M. Wolf, U. Wolf, V. Toronov, A. Michalos, L. A. Paunescu, J. H. Choi, and E. Gratton, “Different Time Evolution of Oxyhemoglobin and Deoxyhemoglobin Concentration Changes in the Visual and Motor Cortices during Functional Stimulation: A Near-Infrared Spectroscopy Study,” NeuroImage 16, 704–712 (2002). [CrossRef]   [PubMed]  

9. M. Firbank, S. R. Arridge, M. Schweiger, and D. T. Delpy, “An investigation of light transport through scattering bodies with non-scattering regions,” Phys. Med. Biol. 41, 767–783 (1996). [CrossRef]   [PubMed]  

10. E. Okada, M. Firbank, M. Schweiger, S. R. Arridge, M. Cope, and D. T. Delpy, “Theoretical and experimental investigation of near-infrared light propagation in a model of the adult head,” Appl. Opt. 36, 21–31 (1997). [CrossRef]   [PubMed]  

11. H. Dehghani, D. T. Delpy, and S. R. Arridge, “Photon migration in non-scattering tissue and the effects on image reconstruction,” Phys. Med. Biol. 44, 2897–2906 (1999). [CrossRef]  

12. F. Scheffold, S. E. Skipetrov, S. Romer, and P. Schurtenberger, “Diffusing-wave spectroscopy of nonergodic media,” Phys. Rev. E 63, 061404-1–11 (2001). [CrossRef]  

13. J. Ripoll, M. Nieto-Vesperinas, S. R. Arridge, and H. Dehghani, “Boundary conditions for light propagation in diffusive media with nonscattering regions,” J. Opt. Soc. Am. A 17, 1671–1681 (2000). [CrossRef]  

14. S. E. Skipetrov and R. Maynard, “Dynamic multiple scattering of light in multilayer turbid media,” Phys. Lett. A 217, 181–185 (1996). [CrossRef]  

15. L. V. Wang, S. L. Jacques, and L. Zheng, “MCML - Monte Carlo modeling of light transport in multi-layered tissues,” Comput. Meth. Progr. Biomed. 47, 131–146 (1995). [CrossRef]  

16. D. A. Boas and A. G. Yodh, “Spatially varying dynamical properties of turbid media probed with diffusing temporal light correlation,” J. Opt. Soc. Am. A 14, 192–215 (1997). [CrossRef]  

17. P. N. Pusey and W. van Megen, “Dynamic Light Scattering by Nonergodic Media,” Physica A 157, 705–741 (1989). [CrossRef]  

18. T. Gisler, H. Rüger, S. U. Egelhaaf, J. Tschumi, P. Schurtenberger, and J. Rička, “Mode-selective dynamic light scattering: theory versus experimental realization,” Appl. Opt. 34, 3546–3553 (1995). [CrossRef]   [PubMed]  

19. S. Feng, F. Zeng, and B. Chance, “Photon migration in the presence of a single defect: a perturbation analysis,” Appl. Opt. 34, 3826–3837 (1995). [CrossRef]   [PubMed]  

20. L. Kou, D. Labrie, and P. Chylek, “Refractive indices of water and ice in the 0.65- to 2.5-μm spectral range,” Appl. Opt. 32, 3531–3540 (1993). [CrossRef]   [PubMed]  

21. A. Ishimaru, Wave propagation and scattering in random media (IEEE Press, New York, 1997).

22. R. Carminati, R. Elaloufi, and J.-J. Greffet, “Beyond the Diffusing-Wave Spectroscopy Model for the Temporal Fluctuations of Scattered Light,” Phys. Rev. Lett. 92, 213903-1–4 (2004). [CrossRef]   [PubMed]