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

Louis Gagnon; Michėle Desjardins; Julien Jehanne-Lacasse; Louis Bherer; Frédéric Lesage

doi:10.1364/OE.16.015514

1. Introduction

The diffusion of near-infrared photons in biological tissue [1, 2, 3, 4, 5] provides a powerful tool to study brain activity [6, 7, 8, 9, 10, 11]. Diffuse correlation spectroscopy (DCS) [12, 13] is an optical imaging technique that enables noninvasive measurements of cerebral blood flow (CBF) in tissue with a relatively high temporal resolution (0.1 s). Its major drawback is a low spatial resolution (1 cm) due to the diffusive nature of photon propagation in tissue. In particular, hemodynamic responses [14, 15, 16, 17, 18], photodynamic therapy [19, 20], blood flow monitoring in tumors [21, 22] and muscle flow responses [23, 24] have been studied in the past with this technique.

The theory behind DCS is based on the correlation transport [12, 25, 26], which can be approximated by a diffusion equation under certain conditions [5, 13, 27, 28, 29]. By assuming that light propagates in an homogeneous medium, the solution of the diffusion equation is simple and as been widely used in the litterature [5, 30]. For in vivo studies applied to the human brain, this assumption does not hold due to the layered structure of the head having different optical properties for each layer [31, 32, 33]. To overcome this problem, depth-sensitive reconstruction and regularization techniques have been developed [18] and tested for the localization of brain activation in the rat. Although this technique is adequate, it requires a large number of measurements (25 sources and 16 detectors in [18]) in order to reconstruct the activation in the tridimensionnal volume.

Another approach is to take into account the layered structure of the head in the model. Solutions of the diffusion equation for layered geometries [34, 35, 36] have been developed and used to fit experimental data in time resolved spectroscopy [37, 38]. Significant differences between properties of the superficial layers of the head and the brain [31, 39] have been obtained by using a two-layered solution to analyze in vivo measurements. In DCS, similar solutions have been developed and validated with Monte Carlo (MC) simulations [5, 17, 30, 40]. A two-layered solution has been used to analyze experimental measurements done on a phantom [30] representing the human head. Also, a three-layered solution has been tested with MC simulations and used to analyse in vivo measurements during brain activity [17]. Both of these studies reported good agreement between the real and the recovered flow parameters. However, the sensitivity of the recovery procedure to the assumed parameters of the model has not been reported in a quantitative way, neither has a quantitative comparison of the efficiencies of the homogeneous versus the multi-layered model to recover flow.

Since errors on the assumed properties of the medium could introduce systematic errors on the recovered flow, we propose in this paper to work with flow changes, which is the quantity of interest in DCS, rather than the absolute flows. We first study the sensitivity of the two-layered recovery procedure to the optical properties of the medium and the thickness of the first layer which are assumed parameters in the model. We then present a quantitative comparison between the ability of the homogeneous and the two-layered model to recover flow changes in DCS, which has not been reported previously.

The derived solution is first validated with MC simulations for two-layered geometries. Then realistic MC simulations done on a segmented MRI anatomical data set and analyzed with the two-layered model are presented. Finally, experimental results obtained using a dynamic double-layer phantom containing a solution of Liposin are analyzed with the two-layered solution.

2. Theory

2.1. Analytical solution

In DCS, the quantity of interest is the electric field temporal autocorrelation function

G_{1} (r, τ) = 〈 E (r, t) E^{*} (r, t + τ) 〉 .

The brackets used here denote an average over time t for an ergodic system and τ is the correlation time. Boas et al. [13] have shown that using the standard diffusion approximation, the correlation transport equation is reduced to the correlation diffusion equation

[D \nabla^{2} - μ_{a} - \frac{1}{3} μ_{s}' k_{0}^{2} 〈 Δ r^{2} (τ) 〉] G_{1} (x, y, z, τ) = - S (x, y, z)

where D=1/3(µ_a+µ′_s) is the diffusion coefficient, c is the velocity of light in the medium, µ′_s and µ_a are respectively the reduced scattering and the absorption coefficient, k ₀ is the light wave vector in the medium, S is the source and 〈Δr ² (τ)〉=6D_Bτ is the mean-squared displacement of the scattering particle in time t, where D_B is the Brownian diffusion coefficient. We assumed that the movement of the particles in the turbid medium is Brownian since it has been reported to provide better fits than gaussian flows on in vivo measurements [30].

We derive the steady-state solution of Eq. (2) for a two-layer medium based on the solution of the time domain diffusion equation developed by Kienle et al. [35, 36, 37]. We assume an infinitely thin beam incident onto the turbid two-layer medium. This beam is scattered isotropically in the upper layer at a depth of z=z ₀ where z₀=1/(µ′_s1+µ _a1). Denote X_i and Gⁱ ₁ (τ) respectively the optical property X and the unnormalized electric field autocorrelation function for layer i. We further assume that the Brownian movement is independent in the first and in the second layer which means that particles can not pass from one layer to the other in the medium. Eq. (2) becomes

\begin{matrix} [D_{1} \nabla^{2} - μ_{a 1} - 2 μ_{s 1}^{'} k_{0}^{2} D_{B 1} τ] G_{1}^{1} (x, y, z, τ) = - δ (x, y, z - z_{0}) & 0 \leq z \leq l, \end{matrix}

\begin{matrix} [D_{2} \nabla^{2} - μ_{a 2} - 2 μ_{s 2}^{'} k_{0}^{2} D_{B 2} τ] G_{1}^{2} (x, y, z, τ) = 0 & l \leq z \end{matrix}

where ℓ is the thickness of the upper layer. D _B1 and D _B2 is the Brownian diffusion coefficient for the first and second layer respectively.

Following the procedure of Kienle et al. [35, 36, 37] and by assuming cylindrical symmetry, the solution of these equations in the first layer is given in the Fourier domain by

{\tilde{G}}_{1}^{1} (s, z, τ) = \frac{\sinh [α_{1} (z_{b} + z_{0})]}{D_{1} α_{1}} \times \frac{D_{1} α_{1} \cosh [α_{1} (ℓ - z)] + D_{2} α_{2} \sinh [α_{1} (ℓ - z)]}{D_{1} α_{1} \cosh [α_{1} (ℓ + z_{b})] + D_{2} α_{2} \sinh [α_{1} (ℓ + z_{b})]}

where α ² _j=(D_js ²+µ_aj+2cµ′_sj k ² ₀ D_Bj)/D_j,s is the radial spatial frequency and

z_{b} = \frac{1 + R_{eff}}{1 - R_{eff}} 2 D_{1} .

R_eff represents the fraction of photons that is internally diffusely reflected at the boundary. It has been calculated in the literature and is about 0.493 for a refractive index of 1.4 [41].

The spatial Fourier inversions, over s, have to be done numerically since no analytical solutions are available. To avoid numerical errors, hyperbolic functions in Eq. (5) are expanded to their exponential form. The two-dimensional Fourier inversion of the previous expression is given by

G_{1}^{1} (ρ, z, τ) = \frac{1}{2 π} \int_{0}^{\infty} {\tilde{G}}_{1}^{1} {(s, z, τ)}_{1} s J_{0} (s ρ) d s

where J ₀ is the Bessel function of zeroth order. The Hankel transform in Eq. (7) is done numerically by using a Gauss-Laguerre quadrature of 5000 points. The nodes and the weights are calculated with MATLAB using the algorithm presented in [42].

2.2. Monte Carlo simulations

The method behind MC simulations for DCS is exactly the same as the one described in [5, 13]. Only the most salient points are described here. It is based on the fact that the normalized correlation function of a field that scatters once from a dilute suspension of noninteracting uncorrelated particles is given by

g_{1}^{s} (τ) = \exp [- \frac{1}{6} q^{2} 〈 Δ r^{2} (τ) 〉]

with q the norm of the momentum transfer imparted by the scattering event and defined by q=|k _out-k _in|=2k ₀ sin (θ/2) with θ the angle between the incoming and the scattered direction. For the case of multiple scattering events, Eq. (8) is summed over each event. By noting that q ²=2k ₀ (1-cosθ), it is easy to compute q ² in the simulation using the normalized wave vector and the definition of the scalar product for two such vectors k_xik_xo+k_yik_yo+k_zik_zo=cosθ. By substituting the definitions of q ² and 〈Δr ² (τ)〉 for Brownian motion in Eq. (8) and summing over each events, we finally come to the following expression for a single photon γ experiencing multiple scattering events over its pathlength in a nonabsorbing medium:

g_{1}^{γ} (τ) = Σ_{j = 1}^{n} \exp [- 2 k_{0}^{2} (1 - \cos θ_{j}) D_{B j} τ] .

Here n denotes the number of individual scattering events j and D_Bj represents the Brownian diffusion coefficient in the region where the event j occurred. In the case of an absorbing medium such as biological tissue, another factor accounting for absorption needs to be added to Eq. (9):

g_{1}^{γ} (τ) = Σ_{i = 1}^{m} Σ_{j = 1}^{n} \exp [- 2 k_{0}^{2} (1 - \cos θ_{j}) D_{B j} τ - μ_{a i} L_{i}]

where m denotes the number of different tissue types i experienced by the photon over its trajectory and L_i is the pathlength in the tissue type i.

The MC simulation propagates and follows the trajectory of one photon at a time in a medium containing different tissue types with each of them having its own optical properties (µ_a, µ′_s, n, g) and its own dynamical property (D_B). The parameter g here refers to the anisotropy parameter used to compute the Henyey-Greenstein phase function [43] for a scattering event. When a photon reaches a detector, Eq. (10) is computed and at the end of the simulations, all the g^γ ₁ (τ) functions for each photon are summed and normalized to obtain g ₁ (τ), the normalized autocorrelation function. Such MC simulations have been done and validated with theory in literature [44, 45, 46]. Many algorithms exist to propagate photons in diffuse media with MC techniques [47, 48, 49, 50, 51] and the one used in this work is a modified version of [50] where the algorithm is described in detail.

2.3. Experimental data

Experimentally, the autocorrelation function of the output of a light detector is the autocorrelation function of the intensity of the light

G_{2} (τ) = 〈 I (t) I (t + τ) 〉 .

For gaussian processes, G ₂ (τ) and G ₁ (τ) defined previously are related by the Siegert relation[52, 53, 54]

G_{2} (τ) = {〈 I 〉}^{2} + β ∣ G_{1} (τ) ∣^{2}

where 〈I〉 is the ensemble-averaged intensity. The parameter β varies with the experimental setup since it depends on the coherence length, the stability and the intensity of the laser [5, 55, 56]. By normalizing to the intensity of light we get:

g_{2} (τ) = 1 + β ∣ g_{1} (τ) ∣^{2}

where g1 (τ) is the normalized electric field autocorrelation function given previously by Eq. (1)

g_{1} (τ) = \frac{G_{1} (τ)}{〈 I 〉} .

3. Methodology

3.1. Simulations

As described previously, the code tMCimg [50] has been modified to take into account the dynamic properties of the turbid medium. We simulated 10⁸ photons for each simulation. This number has been reported to give good photon statistics in media with optical properties equivalent to biological tissue [50, 57, 58]. Simulations took about two hours on a 3 GHz dual-core processor. The resolution of the three-dimensional volume was 1 mm×1 mm×1 mm for each simulation performed for this work. MATLAB with the function lsqcurvefit was used to fit the simulated data with the theoretical model.

Simulations were run in a two-layered medium which consisted of a rectangular slab over a semi-infinite medium. In order to test what factors cause errors on the recovery procedure, we introduced errors in optical properties assumed by the fits. We studied the effect of a wrong absorption and reduced scattering coefficients in one layer at a time and also that of a wrong thickness for the first layer.

3.2. Human head

DCS is a technique which is mainly used to measure relative flow changes. We were interested in measuring the efficiency of the two-layered model to extract relative cerebral blood flow (CBF) changes. To build a medium representing the human head, a Siemens Trio 3T Magnetic Resonance Imaging (MRI) scanner with a MPRage sequence was used to get an anatomical image. The exam was performed on a 26-year-old female volunteer and written consent from the subject was obtained. This exam was approved by the ethic comity of the Institut Universitaire de Gériatrie de Montréal. This anatomical data set was segmented into skin, skull, cerebrospinal fluid (CSF) and finally, white and gray matter of brain tissue. The CSF, white and gray matter were extracted from the whole image with SPM5 (UCL). The skin and the skull were then extracted semi-automatically with a MATLAB algorithm. Optical properties were then assigned to each tissue type as done in refs. [57, 58]. The optical properties used in the simulation are listed at Table 1. The optical source and detectors were located over the motor region as shown in Fig. (1).

3.3. Experimental setup

The experimental setup is similar to the one built by Durduran et al. [30]. It consists of a 810 nm laser (CrystaLaser, Reno, NV, USA), a 300 microns core mutli-mode optical fiber to inject light and four monomode fibers used to collect it. Four photon-counting avalanche photodiodes were connected to an 8 channel autocorrelator (Correlator.com, Bridgewater, NJ, USA).

3.4. Phantom

The dynamic phantom was built with acrylic glass and consisted of two compartments filled with a solution containing 1 % of Liposin, as illustrated in Fig. (2). The division between the

Table 1. Optical properties used in the MC simulations of the segmented human head. [57]

View Table | View all tables in this article

two compartments was built with a thin plastic membrane to avoid back-reflections and was removable so we were able to adjust the width of the first compartment. The source and detector optical fibers were fixed to the phantom with a piece of rubber that was glued to the acrylic glass with silicone. To validate the acquisition system, a small magnet stalk was placed in the second compartment and the whole phantom was positioned over a magnetic agitator plate to modify the flow as desired in the second compartment. In order to introduce realistic flow changes, we removed the phantom from the agitator and introduced glycerol in the Liposin solution of the second compartment to raise its viscosity and reduce its Brownian diffusion coefficient [14].

Fig. 1. Localisation of the optical probe on the MRI anatomical data. The human head was segmented in skin, skull, CSF, white matter and gray matter. Source and detectors are illustrated by red and green arrows respectively.

Download Full Size | PDF

Fig. 2. Schematic view of the dynamical phantom. The division which separates the two compartments is mobile allowing different thicknesses for the first layer. The two compartments are filled with a 1 % Liposin solution. The optical fibers are fixed to the phantom by a rubber band glued to its side.

Download Full Size | PDF

4. Results

4.1. Monte Carlo simulations

A comparison between Monte Carlo simulations and the analytical model was first done to validate the analytical solution. The comparison is shown on Fig. (3). MC simulations were done on a two-layered medium. The thickness of the first layer was 10 mm while source-detector (S-D) distances ranged from 10 to 30 mm. The absorption coefficient µ_a and the effective scattering coefficient µ′_s were respectively 0.0074 and 1.28 mm^-1 for both layers while the Brownian diffusion coefficient D_B was 10^-8 mm²/s in the first layer and 10^-6 mm²/s for the other semi-infinite part of the medium.

Having validated the analytical model, we then moved to the problem of parameter estimation. The recovery procedure has been tested on two separateMC data sets and computed at two distances. Table 2 shows the Brownian diffusion coefficient simulated and the recovered values using the two-layered analytical model in the fitting procedure. In these fits, the thickness of the first layer and the optical properties were assumed to be known. The S-D distances were chosen to be 15 and 20 mm given that these represent realistic distances for in vivo measurements [16, 30]. For comparison, results recovered using a homogeneous fit are also provided.

Since absolute flow recovery is hard to achieve experimentally, we tested the efficency of the two-layered model to recover relative flow changes. We simulated an increase in the dynamic properties of the second layer by performing two different simulations. In the first one (baseline), we set the Brownian diffusion coefficient to 10^-8 in the first layer and 10^-7 in the second one. In the second simulation (increased flow), we kept the D_B to 10^-8 in the first layer but we set it to 1.5×10^-7 in the second layer which corresponds to a flow increase of 50 %. We estimated the flow change by comparing the D_B recovered from the fits of the baseline and the increased flow simulations. We applied the procedure both with the two-layered and the homogeneous model. The recovered flow changes are presented in Fig. (4).

Fig. 3. Comparison of MC simulations with the analytical model of Eqs. (5) and (7) on a two-layered medium. The simulations are represented by points while the analytical model by solid lines. The thickness of the first layer was 10 mm and the source-detector (S-D) distances varied from 10 to 30 mm. The optical properties (µ_a=0.0074 mm^-1 and µ′_s=1.28 mm^-1) were constant for the two layers while the Brownian diffusion coefficient D_B was 10^-8 mm²/s and 10^-6 mm²/s for the first and second layer respectively.

Download Full Size | PDF

4.2. Sensitivity of the fitting procedure

Since the absorption and the reduced scattering coefficients are usually unknowns and have to be measured experimentally, we tested the errors introduced by assuming wrong absorption or scattering in the medium. The upper part of Fig. (5) shows the error on the recovered flow changes assuming different absorption coefficients µ_a ranging from 0.005 to 0.015 while the real value was 0.01 mm^-1. The procedure has been tested assuming wrong absorption in each layer independently. For these simulations, the effective scattering coefficient µ′_s was 1.28 mm^-1 for both layers and the thickness of the first layer was 10 mm. The S-D distance was 15 mm. The same procedure has been applied for the reduced scattering coefficient µ′_s and the resulting errors are presented in the lower part of Fig. (5). In these last simulations, µ _a was set to 0.01 mm^-1 in both layers.

The sensitivity to the error on the thickness of the first layer was also tested. Flow changes have been recovered on a MC data set assuming different thicknesses for the first layer, ranging from 9 to 11 mm while the real thickness was 10 mm. In this simulation, µ_a=0.01 mm^-1 and µ′_s=1.28 mm^-1 for both layers while the S-D distance was 15 mm. Results are presented in Fig. (6).

4.3. Human head simulations

In order to test the efficiency of the two-layered model to extract changes in cerebral blood flow (CBF) in vivo, simulations on a segmented human head have been performed. For these simulations, the skin, skull, cerebrospinal fluid (CSF), the white and the gray matter have been segmented from the anatomical MRI data set. Simulations have been performed over the motor

Table 2. Brownian diffusion coefficients recovered using the two-layered analytical model in the fitting procedure. In these fits, the thickness of the first layer and the optical properties was assumed to be known. In these simulations, µ′_s=1.28 mm^-1 for both layers and µ_a was respectively 0.008 and 0.0176 mm^-1 for the first and second layers. The thickness of the first layer was 10 mm. Results obtained using a homogeneous model are also shown for comparison.

View Table | View all tables in this article

region of the brain and S-D distances ranged from 10 to 30 mm. Optical properties used in the simulations (taken from [57]) are shown in Table 1.

We performed two different simulations. In the first one (baseline), we set the Brownian diffusion coefficient to 10^-8 in the superficial layers (skin, skull and CSF) and to 10^-7 in the brain. In the second one (increased flow), we let it to 10^-8 in the superficial layers but we increase it to 1.5×10^-7 in the brain which corresponds to a 50 % increase. We fitted both the baseline and the increased flow simulated data with the two-layered and the homogeneous model and compared the increase in the recovered parameters between the two simulations. The changes in the fitted D_B for each layer of the two-layered model and for the homogeneous D_B are reported in Fig. (7).

4.4. Phantom validation

Experimental data were taken on the two-layered dynamical phantom at rest and during induced agitation in the second layer. For these measurements, the liquid located in the first layer was at rest while one in the second layer was agitated with a magnetic agitator. The obtained curves are presented in Fig. (8). The thickness of the first layer of the phantom was 10 mm and the S-D distances were 10 and 20 mm.

To provide realistic flow changes, we used glycerol to increase the viscosity of the Liposin solution [14]. We first removed the mobile division of the phantom in order to work with an homogeneous medium. We then introduced glycerol and monitored the D _B variations with DCS measurements and analyzed them with the homogeneous model. We measured that a glycerol concentration of 6 % decreases the D _B value by 37 %. After putting back the mobile division, we took two measurements: one with the Liposin solution in both compartments and the other one with the second compartment filled with the Liposin solution containing 6 % of glycerol. Those measurements were analyzed both with the homogeneous and the two-layered model. The recovered flow changes are presented in Fig. (9).

5. Discussion

In this study, our aim was to quantify the efficiency of DCS applied to brain imaging when the layered geometry of the human head is taken into account. We wanted to evaluate whether using a two-layered model could improve the efficiency to recover changes in CBF, by reducing partial volume error.

Fig. 4. Changes in D_B recovered from MC simulations both with the two-layered and the homogeneous model. These changes were estimated by recovering the D_B of two simulations. In the first one (baseline), we set the Brownian diffusion coefficient to 10^-8 in the 10 mm thick first layer and to 10^-7 in the second one. In the second simulation (increased flow) we kept the D_B to 10^-8 in the first layer but increased it to 1.5×10^-7 in the second layer which corresponds to an increase of 50 %. The source-detector (S-D) distances varied from 10 to 30 mm and the optical properties (µ_a=0.01 mm^-1 and µ′_s=1.28 mm^-1) were constant for the two layers in both simulations (baseline and increased flow).

Download Full Size | PDF

5.1. Simulations and fits

As seen in Fig. (3), theMC simulations agree well with the analytical solution of the model. The presence of the second layer introduces a change in the early decay of the correlation function. This change is small for short S-D distances (10 mm) but increases as it becomes larger. This was expected due to the weak contribution of the second layer for short S-D distances since most of the photons that reach the detector have traveled mainly in the first layer. A similar increase in contribution of the second layer has been observed by reducing the thickness of the first layer while keeping the S-D distance constant.

We also observe in Fig. (3) that the agreement between the MC simulations and the analytical model is better for long S-D distances. This was expected since the diffusion approximation, which the analytical model is based on, is not as accurate for short S-D distances. Our MC code has also been tested with an homogenous medium and compared with an homogeneous analytical model [30]. The agreement was as good as for the two-layered model. All these results validate our MC code with the two-layered analytical model.

5.2. Recovery procedure

In simulations, the errors on the recovered D_B with the two-layered model presented in Table 2 are less than 9 % which is comparable to errors on the recovered optical parameters in TD spectroscopy [37, 38, 31]. A fit with an homogeneous model is also shown for comparison in the last column. In this case, the Brownian diffusion coefficients recovered with the homogeneous fit are closer to the first layer values even for S-D distances of 20 mm. This suggests that an homogeneous fitting procedure is more sensitive to the dynamic properties of the superficial layers. For a two-layered medium with a larger flow in the second layer, one underestimates the Brownian diffusion coefficient of the second layer when fitting the data with an homogeneous model.

Fig. 5. Error on the recovered flow change due to the estimation of the absorption µ_a and reduced scattering µ′_s coefficient for each layer. The thickness of the first layer was 10 mm and the S-D distance was 15 mm. Top) Error on the assumed absorption. Bottom) Error on the assumed scattering.

Download Full Size | PDF

Fig. 6. Error on the recovered flow changes due to a wrong estimation of the thickness of the first layer. For this simulation, the real absorption coefficient was 0.01 mm^-1 and the effective scattering coefficient µ′_s was 1.28 mm^-1 for both layers. The thickness of the first layer was 10 mm and the S-D distance was 15 mm.

Download Full Size | PDF

Fig. 7. Recovered CBF changes in each layer with the two-layered model following an increase of 50 % of the Brownian diffusion coefficient in the brain tissue (white and gray matter). Again, the source and the detectors were located over the motor region (Fig. (1)).

Download Full Size | PDF

Fig. 8. Experimental DCS measurements done on a two-layered dynamical phantom. The set up consisted in two detectors located at 10 mm from the source and two others at 20 mm. The blue curves were obtained while the liquid was at rest in the phantom and the red ones when agitation was induced in the second layer by the magnetic agitator. The thickness of the first layer of the phantom was 10 mm while the S-D distances were 10 and 20 mm.

Download Full Size | PDF

Fig. 9. Decreases in recovered D_B measured on the phantom after the introduction of the Liposin solution containing 6 % of glycerol in the second compartment. The data was analyzed both with the homogeneous and the two layered model. The thickness of the first layer of the phantom was 10 mm while the S-D distances were 10 and 20 mm. The real D_B decrease of 37 % has been measured experimentally on an homogeneous phantom using the homogeneous model to analyze the data.

Download Full Size | PDF

The errors on the recovered D_B increase rapidly when the assumed optical properties are wrongly estimated which makes absolute flow measurements hard to recover. To bypass this difficulty, we focussed our study on the recovery of flow changes, which is the quantity of interest in DCS. The recovered flow changes presented in Fig. (4) show that the two-layered model performs better than the homogeneous one when estimating flow changes in the second layer. The recovered values with the homogeneous model are 5 to 10 times lower than the real simulated increase (50 %) while the changes recovered with the two-layered model are only few percent below 50 %. This result raises questions about using an homogeneous model to measure changes in CBF in the human brain given the layered structure of the tissue.

5.3. Sensitivity

The sensitivity of the relative flow change recovery procedure to the optical properties was computed. The errors introduced on the recovered flow change in the first layer are illustrated by circles in Fig. (5). Since absorption and scattering do not affect the shape of the correlation function in regions where the first layer contributes the most (long correlation times), a smaller error is introduced in the fitting procedure for this layer. Errors in the second layer are shown in the same figure. An error in the absorption coefficient of either the first or the second layer introduces a larger error in the recovered flow change in the second layer compared to the first layer. Since the contribution of the second layer in the correlation function is low at those short S-D distances, a small change in the absorption of the medium introduces large changes in the contribution of this layer. Separately, an error in the absorption of the second layer will induce a larger error in the recovered flow change and again, due to the low contribution of this layer in the correlation function. All recovered errors are below 5 % except for an error of 50 % in the absorption coefficient of the second layer. On the other hand, errors on the reduced scattering coefficient induce only small errors in the recovered flow change. The scattering coefficient multiplies the Brownian diffusion coefficient in the analytical solution and an error on its assumption mainly rescales the recovered D_B. Since we measured the relative change between two curves, baseline and flow condition, a shift of both curves does not induce a large error.

Finally, the flow change recovered in the second layer is more sensitive to the first layer thickness than the one in the first layer as seen in Fig. (6). Again, this is due to the lower contribution of the second layer in the correlation function. The errors are below 15 % for a ±1 mm thickness error.

5.4. Human head simulations

Since DCS has been mainly used to measure relative CBF changes [30], we performed a simulation to measure the efficiency of the two-layered model to extract relative changes in dynamic properties of the human brain. The results presented in Fig. (7) show that the changes in flow recovered with the two-layered model are closer to the real simulated change (50 %). This was expected since a two-layered model reduces the partial volume error. The homogeneous model underestimates the flow change by assuming that the flow was constant everywhere in the medium. A correction factor has been computed with simulations on a two-layered medium in the literature [16, 30] to take into account this effect and the obtained value was 5. However, our results with homogeneous fits in Fig. (7) show that this correction factor should be higher (10-20). The two-layered model reduces this artefact by considering regions with different dynamic properties in the medium. We also see in Fig. (7) that flow changes recovered for the first layer are small (about a few %) which is close to the simulated value (0 %). These results confirm that the two-layered model improves the quantification of DCS measurements and suggest that it performs better than an homogeneous model to fit in vivo measurements when looking at relative changes in CBF.

5.5. Experimental data

Experimental measurements shown in Fig. (8) agree with simulated ones of Fig. (3). For short S-D distances, an increase in flow in the second layer makes the early part of the autocorrelation function decrease faster. The long correlation time decays are not affected since the contribution of photons from the second layer in the autocorrelation fonction is weak for long correlation time. By increasing the S-D distance, a larger part of the autocorrelation function is affected by the flow increase in the second layer, as predicted by the 30 mm simulation of Fig. (3). The effect is observed to be larger for the experimental data since the absorption in the Liposin solution was low which in turn increases the contribution of the second layer.

The flow changes recovered with the two-layered model for the two configurations (with and without glycerol) presented in Fig. (9) are in good agreement with the real value. However, the flow changes recovered with the homogeneous model are 7 to 15 times lower than the real change. The recovered values for the first layer only show small changes. These experimental results are in agreement with the ones obtained by MC simulations and suggest that a two-layered model should be used to measure CBF changes in the brain in order to reduce partial volume errors.

Crosstalk between the two layers in the two-layered recovery seems larger for large S-D distances. This is due to the lower contribution of the first layer in the 20 mm S-D measurement since most of the photons travel in the second layer. However, this crosstalk is small since the flow changes recovered for the first layer are below 3 %.

In this study, we assumed that the optical properties (µ_a, µ′_s) of the medium were not altered by the introduction of glycerol as in [14]. In order to verify this assumption, a time resolved system could be used to extract the absorption and the effective scattering coefficient of the Liposin solution before and after the introduction of glycerol.

5.6. Comparison with DCT

Diffuse correlation tomography (DCT) has been used to reconstruct images of blood flow at different depths from many S-D measurements [18]. The procedure relies on the Rytov approximation [59] and a depth-sensitive regularization technique [60]. The two-layered model used in this work also allows depth-sensitive recovery. The main difference between the two methods is the number of measurements needed which is lower with the two-layered approach. Moreover, DCT provides a flow reconstruction in a thin slice while two-layered models assume that the second layer is semi-infinite.

6. Conclusion

We have provided a quantitative comparison between an homogeneous and a two-layered model to extract efficiently flow changes in layered geometries. We first tested the two-layered solution with MC simulations and showed good agreement between the analytical model and the simulated data. A fit-based recovery procedure has been tested and errors of only few % have been obtained for flow change measurements. Our work indicates that the homogeneous model is more sensitive to superficial layers and that the two-layered model performs better at estimating flow changes in deep layers. We built a two-layered dynamical phantom and we took experimental DCS measurements while changing the flow in the second layer. Clear distinctions between the change of flow in each of the two layers have been obtained by fitting the experimental data with the two-layered analytical model.

We investigated the possibility of applying this model to in vivo measurements by performing MC simulation on a MRI human head anatomical data set. We showed that using a two-layered model to analyze relative changes in CBF performs better than the homogeneous case. The two-layered model improves the quantification of DCS measurements by reducing partial volume errors. This study suggests that a two-layered model should be used to quantify relative changes in flow measurements.

Acknowledgments

We thank David Boas for helpful comments and discussions. This work has been supported by the Natural Sciences and Engineering Research Council of Canada (NSERC). L. Gagnon and M. Desjardins are supported by the NSERC Postgraduate Scholarship and F. Lesage by the NSERC Discovery Grant.

References and links

1. M. S. Patterson, B. Chance, and C. Wilson, “Time resolved reflectance and transmittance for the noninvasive measurement of tissue optical properties,” Appl. Opt. 28, 2331–2336 (1989). [CrossRef] [PubMed]

2. B. J. Tromberg, L. O. Svaasand, T. Tsay, and R. C. Haskell, “Properties of photon density waves in multiplesscattering media,” Appl. Opt. 32, 607–616 (1993). [CrossRef] [PubMed]

3. D. Boas, M. A. O’Leary, B. Chance, and A. G. Yodh, “Scattering and wavelength transduction of diffuse photon density waves,” Phys. Rev. E 47 (1993). [CrossRef]

4. D. Boas, M. A. O’Leary, B. Chance, and A. G. Yodh, “Scattering of diffuse photon density waves by spherical inhomogeneities within turbid media: Analytic solution and applications,” Proc. Natl. Acad. Sci. USA 91, 4887–4891 (1994). [CrossRef] [PubMed]

5. D. A. Boas, “Diffuse photon probes of structural and dynamical properties of turbid media : theory and biomedical applications,” Ph.D. thesis, University of Pennsylvania (1996).

6. A. Yodh and B. Chance, “Spectroscopy and Imaging with Diffusing Light,” Phys. Today 48, 34–40 (1995). [CrossRef]

7. A. Villringer, C. Hock, L. Schleinkofer, and U. Dirnagl, “Near Infrared Spectroscopy (NIRS): A New Tool to Study Hemodynamic Changes During Activation of Brain Function in Human Adults,” Neurosci. Lett. 154, 101–104 (1993). [CrossRef] [PubMed]

8. Y. Hoshi and M. Tamura, “Multichannel near-infrared optical imaging of human brain activity,” J. Appl. Physiol. 75, 1842–1846 (1993). [PubMed]

9. H. Obrig and A. Villringer, “Beyond the visible - imaging the human brain with light,” J. Cereb. Blood Flow Metab. 23, 1–18 (2002). [PubMed]

10. M. Franceschini and D. Boas, “Noninvasive measurement of neuronal activity with near-infrared optical imaging,” NeuroImage 21, 372–386 (2004). [CrossRef] [PubMed]

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

12. D. Boas, L. E. Campbell, and A. G. Yodh, “Scattering and Imaging with Diffusing Temporal Field Correlations,” Phys. Rev. Lett.75 (1995). [CrossRef] [PubMed]

13. D. Boas and A. G. Yodh, “Spatially varying dynamical properties of turbid media probed with diffusing temporal light correlation,” J. Opt. Soc. Am. A14 (1997).

14. C. Cheung, J. P. Culver, K. Takahashi, J. H. Greenberg, and A. G. Yodh, “In vivo cerebrovascular measurement combining diffuse near-infrared absorption and correlation spectroscopies,” Phys. Med. Biol. 46, 2053–2065 (2001). [CrossRef] [PubMed]

15. J. P. Culver, T. Durduran, D. Furuya, C. Cheung, J. H. Greenberg, and A. G. Yodh, “Diffuse optical tomography of cerebral blood flow, oxygenation, and metabolism in rat during focal ischemia,” J. Cereb. Blood Flow Metab. 23, 911–924 (2003). [CrossRef] [PubMed]

16. 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]

17. J. Li, G. Dietsche, D. Iftime, S. E. Skipetrov, G. Maret, T. Elbert, B. Rockstroh, and T. Gisler, “Noninvasive detection of functionnal brain activity with near-infrared diffusing-wave spectroscopy,” J. Biomed. Opt.10, 044,002 (2005). [CrossRef] [PubMed]

18. C. Zhou, G. Q. Yu, D. Furuya, J. H. Greenberg, A. G. Yodh, and T. Durduran, “Diffuse optical correlation tomography of cerebral blood flow during cortical spreading depression in rat brain,” Opt. Express 14, 1125–1144 (2006). [CrossRef] [PubMed]

19. G. Yu, T. Durduran, H. W. Wang, C. Zhou, E. M. Putt, H. M. Saunders, C. M. Sehgal, E. Glatstein, A. G. Yodh, and T. M. Bush, “Noninvasive Monitoring of Murine Tumor Blood Flow During and After Photodynamic Therapy Provides Early Assessment of Therapeutic Efficacy,” Clin. Cancer Res. 11, 3543–3552 (2005). [CrossRef] [PubMed]

20. G. Yu, T. Durduran, C. Zhou, T. C. Zhu, J. C. Finlay, T. M. Busch, S. B. Malkowicz, S. M. Hahn, and A. G. Yodh, “Real-time In Situ Monitoring of Human Prostate Photodynamic Therapy with Diffuse Light,” Photochem. Photobiol. 82, 1279–1284 (2006). [CrossRef] [PubMed]

21. T. Durduran, R. Choe, G. Yu, C. Zhou, J. C. Tchou, B. J. Czerniecki, and A. G. Yodh, “Diffuse optical measurement of blood flow in breast tumors,” Opt. Lett. 30, 2915–2917 (2005). [CrossRef] [PubMed]

22. U. Sunar, H. Quon, T. Durduran, J. Zhang, J. Du, C. Zhou, G. Yu, R. Choe, A. Kilger, R. Lustig, L. Loevner, S. Nioka, B. Chance, and A. G. Yodh, “Non-invasive diffuse optical measurement of blood flow and blood oxygenation for monitoring radiation therapy in patients with head and neck tumors,” J. Biomed. Opt.11, 064,021 (2006). [CrossRef]

23. G. Yu, T. Durduran, G. Lech, C. Zhou, B. Chance, E. R. Mohler, and A. G. Yodh, “Time-dependant Blood Flow and Oxygenation in Human Skeletal Muscle Measured with Noninvasive Near-infrared Diffuse Optical Spectroscopies,” J. Biomed. Opt.10, 024,027 (2005). [CrossRef] [PubMed]

24. G. Yu, T. F. Floyd, T. Durduran, C. Zhou, J. Wang, J. A. Detre, and A. G. Yodh, “Validation of diffuse correlation spectroscopy for muscle blood flow with concurrent arterial spin labeled perfusion MRI,” Opt. Express15 (2007). [CrossRef] [PubMed]

25. B. J. Ackerson, R. L. Dougherty, N. M. Reguigui, and U. Nobbman, “Correlation transfer: application of radiative transfer solution methods to photon correlation problems,” J. Thermophys. Heat Transfer 6, 577–588 (1992). [CrossRef]

26. R. L. Dougherty, B. J. Ackerson, N. M. Reguigui, F. Dorri-Nowkoorani, and U. Nobbman, “Correlation transfer: development and application,” J. Quant. Spectrosc. Radiat. Transfer 52, 713–727 (1994). [CrossRef]

27. S. Chandrasekhar, Radiative Transfer (Dover, New York, 1960).

28. A. Ishimaru, Wave propagation and Scattering in Random Media, vol. 1 (New York : Academic, 1978).

29. S. R. Arridge, “Optical tomography in medical imaging,” Inverse Problem 15, R41–R93 (1999). [CrossRef]

30. T. Durduran, “Non-Invasive Measurements of Tissue Hemodynamics with Hybrid Diffuse Optical Methods,” Ph.D. thesis, University of Pennsylvania (2004).

31. L. Gagnon, C. Gauthier, J. Selb, D. A. Boas, R. D. Hoge, and F. Lesage, “Double layer estimation of intra- and extra-cerebral hemoglobin concentration with a time-resolved system,” J. Biomed. Opt.13 (2008). [CrossRef] [PubMed]

32. L. Gagnon, J. Selb, D. A. Boas, R. D. Hoge, and F. Lesage, “Measurements of Hemoglobin Concentrations in the Human Forehead Using Time-Resolved Reflectance,” in Biomedical Optics, vol. 1 of OSA Technical Digest (CD) (Optical Society of America, 2008). Paper BSuE73.

33. F. Lesage, L. Gagnon, and M. Dehaes, “Diffuse Optical-MRI fusion and applications,” Proc. SPIE 6850, C1–C11 (2008).

34. F. Martelli, A. Sassaroli, S. Del-Bianco, Y. Yamada, and G. Zaccanti, “Solution of the time-dependent diffusion equation for layered diffusive media by the eigenfucntion method.” Phys. Rev. E 67, 056,623 (2003). [CrossRef]

35. A. Kienle, M. S. Patterson, N. Dognitz, R. Bays, G. Wagnieres, and H. van den Bergh, “Noinvasive determination of the optical properties of two-layered turbid media,” Appl. Opt. 37, 779–791 (1998). [CrossRef]

36. A. Kienle, T. Glanzmann, G. Wagnieres, and H. van den Bergh, “Investigation of two-layered turbid media with time-resolved reflectance,” Appl. Opt. 37, 6852–6862 (1998). [CrossRef]

37. A. Kienle and T. Glanzmann, “In vivo determination of the optical properties of muscle with time-resolved reflectance using a layered model,” Phys. Med. Biol. 44, 2689–2702 (1999). [CrossRef] [PubMed]

38. F. Martelli, S. D. Bianco, G. Zaccanti, A. Pifferi, A. Torricelli, A. Bassi, P. Taroni, and R. Cubeddu, “Phantom validation and in vivo application of an inversion procedure for retrieving the optical properties of diffusive layered media from time-resolved reflectance measurements.” Opt. Lett.29 (2004). [CrossRef] [PubMed]

39. M. Choi, V. Wolf, U. Toronov, C. Wolf, D. Polzonetti, L. P. Hueber, R. Safonova, A. Gupta, W. Michalos, E. Mantulin, and Gratton, “Noninvasive determination of the optical properties of adult brain: near-infrared spectroscopy approach,” J. Biomed. Opt. 9, 221–229 (2004). [CrossRef] [PubMed]

40. C. Zhou, “In vivo optical imaging and spectroscopy of cerebral hemodynamics,” Ph.D. thesis, University of Pennsylvania (2007).

41. R. C. Haskell, L. O. Svaasand, T. T. Tsay, T. C. Feng, M. McAdams, and B. J. Tromberg, “Boundary conditions for the diffusion equation in radiative transfer,” J. Opt. Soc. Am. A 11, 2727–2741 (1994). [CrossRef]

42. G. H. Golub and J. H. Welsch, “Calculation of Gauss quadrature rules,” Math. Comput. 23, 221–230 (1969). [CrossRef]

43. L. G. Henyey and J. Greenstein, “Diffuse radiation in the galaxy,” Astrophys. J. 93, 70–83 (1941). [CrossRef]

44. A. A. Middleton and D. S. Fisher, “Discrete scatterers and autocorrelations of multiply scattered light,” Phys. Rev. B 43, 5934–5938 (1991). [CrossRef]

45. D. J. Durian, “Accuracy of diffusing-wave spectroscopy theories,” Phys. Rev. E 51, 3350–3358 (1995). [CrossRef]

46. M. H. Koelink, F. F. M. de Mul, J. Greve, R. Graaff, A. C. M. Dassel, and J. G. Aarnoudse, “Laser Doppler blood flowmetry using two wavelengths: Monte Carlo simulations and measurements,” Appl. Opt. 33, 3549–3558 (1994). [CrossRef] [PubMed]

47. R. Graaff, M. H. Koelink, F. F. M. de Mul, W. G. Zijlstra, A. C. M. Dassel, and J. G. Aarnoudse, “Condensed Monte Carlo simulations for the description of light transport,” Appl. Opt. 32, 426–434 (1993). [CrossRef] [PubMed]

48. S. L. Jacques and L. Wang, “Monte Carlo modeling of light transport in tissues,” in Optical-Thermal Response of Laser-Irradiated Tissue, pp. 73–100 (Plenum, New York,1995).

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

50. D. Boas, J. Culver, J. Stott, and A. Dunn, “Three dimensional Monte Carlo code for photon migration through complex heterogenous media including the adult human head,” Opt. Express10 (2002). [PubMed]

51. L. Wang and H. I. Wu, Biomedical Optics: Principles and Imaging (John Wiley & Sons, Hoboken, 2007).

52. S. O. Rice, “Mathematical analysis of random noise,” in Noise and Stochastic Processes, p. 133 (Dover, New York, 1954).

53. 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]

54. D. J. Pine, D. A. Weitz, P. M. Chaikin, and E. Herboltzheimer, “Diffusing-wave spectroscopy,” Phys. Rev. Lett. 60, 1134–1137 (1988). [CrossRef] [PubMed]

55. P. N. Pusey, J. M. Vaughan, and D. V. Willet, “Effect of spatial incoherence of the laser in photon-counting spectroscopy,” J. Opt. Soc. Am. 73, 1012–1017 (1983). [CrossRef]

56. T. Bellini, M. A. Glaser, N. A. Clark, and V. Degiorgio, “Effects of finite laser coherence in quasielastic multiple scattering,” Phys. Rev. A 44, 5215–5223 (1991). [CrossRef] [PubMed]

57. G. Strangman, M. A. Franceschini, and D. A. Boas, “Factors affecting the accuracy of near-infrared spectroscopy concentration calculations for focal changes in oxygenation parameters,” NeuroImage 18, 865–879 (2003). [CrossRef] [PubMed]

58. T. Huppert, R. Hoge, A. M. Dale, M. Franceschini, and D. Boas, “Quantitative spatial comparison of diffuse optical imaging with blood oxygen level-dependent and arterial spin labeling-based functional magnetic resonance imaging,” J. Biomed. Opt.11 (2006). [CrossRef] [PubMed]

59. A. Kak and M. Slaney, Principles of Computerized Tomographic Imaging (IEEE Press, New York, 1988).

60. B. W. Pogue, T. O. McBride, J. Prewitt, U. L. Osterberg, and K. D. Paulsen, “Spatially variant regularization improves diffuse optical tomography,” Appl. Opt. 38, 2950–2961 (1999). [CrossRef]

S-D (mm)	layer	real D_B (mm²/s)	recovered D_B (mm²/s)	error (%)	homogeneous D_B (mm²/s)
15	1 ^st	1e-008	9.34e-009	6.56	1.43e-008
	2 ^nd	1e-006	1.05e-006	5.43
20	1 ^st	1e-008	9.46e-009	5.39	1.94e-008
	2 ^nd	1e-006	1.08e-006	8.40

S-D (mm)	layer	real D_B (mm²/s)	recovered D_B (mm²/s)	error (%)	homogeneous D_B (mm²/s)
15	1 ^st	1e-008	9.34e-009	6.56	1.43e-008
	2 ^nd	1e-006	1.05e-006	5.43
20	1 ^st	1e-008	9.46e-009	5.39	1.94e-008
	2 ^nd	1e-006	1.08e-006	8.40

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

Abstract

1. Introduction

2. Theory

2.1. Analytical solution

2.2. Monte Carlo simulations

2.3. Experimental data

3. Methodology

3.1. Simulations

3.2. Human head

3.3. Experimental setup

3.4. Phantom

4. Results

4.1. Monte Carlo simulations

4.2. Sensitivity of the fitting procedure

4.3. Human head simulations

4.4. Phantom validation

5. Discussion

5.1. Simulations and fits

5.2. Recovery procedure

5.3. Sensitivity

5.4. Human head simulations

5.5. Experimental data

5.6. Comparison with DCT

6. Conclusion

Acknowledgments

References and links

Cited By

Figures (9)

Tables (2)

Equations (15)

Optics Express

layer	µ^′_s (mm-1)	µ_a (mm-1)
skin	0.74	0.0164
skull	0.81	0.0115
CSF	0.01	0.0017
white matter	1.16	0.017
gray matter	1.16	0.017