## Abstract

Diffuse optical correlation methods were adapted for three-dimensional (3D) tomography of cerebral blood flow (CBF) in small animal models. The image reconstruction was optimized using a noise model for diffuse correlation tomography which enabled better data selection and regularization. The tomographic approach was demonstrated with simulated data and during *in-vivo* cortical spreading depression (CSD) in rat brain. Three-dimensional images of CBF were obtained through intact skull in tissues deep (∼ 4 mm) below the skull surface.

©2006 Optical Society of America

## 1. Introduction

Normal brain function depends on delivery of oxygen and glucose, and on clearance of the by-products of metabolism. Thus, an understanding of the normal and pathologic conditions of oxygen supply and consumption, and measurement of blood flow is important for clinical applications [1]. To this end a variety of tools have been developed to image cerebral blood flow, but all of these techniques have limitations. For example, PET-based [2] and MRI-based [3–5] cerebral blood flow measurements are expensive and sometimes lack the spatio-temporal resolution important for animal studies. Laser speckle imaging [6–9] and scanning laser-Doppler flowmetry [10] have high spatial resolution, but are limited to two-dimensional (2D) mapping of cerebral blood flow and require the skull to be removed or thinned during the study.

The focus of this paper is the diffuse optical method. Diffuse optical imaging and spectroscopy is a growing subfield of biomedical optics whose aim is to investigate physiology millimeters to centimeters below the tissue surface [11–15]. Thus far diffuse optical methods have been used in research and clinical settings for measurement of blood volume, blood oxygen saturation [16–29], and, to a lesser extent, blood flow [28–34]. The basics of diffuse correlation tomography (DCT) have already been developed and tested in phantom studies [35–37], and two-dimensional image slices have been obtained below the tissue surface in a three-dimensional (3D) rat brain ischemia model [28]. However, to our knowledge, 3D *in-vivo* images of dynamic changes in cerebral blood flow using the diffuse correlation method have not been demonstrated.

The principle of diffuse correlation tomography is similar to that of diffuse optical tomography (DOT), which is well established for mapping three-dimensional tissue optical properties [11,14]. However, in practice, blood flow imaging using DCT is harder to obtain due to its sensitivity to measurement noise and data selection. In this contribution we describe an analysis to optimize data selection and image reconstruction for DCT, and we use this optimization scheme to derive 3D tomographic images of flow dynamics from an *in-vivo* model of cortical spreading depression (CSD). CSD is a wave of excitation and depolarization of neuronal cells that spreads radially with a speed of 2-5 mm/min over the cerebral cortex [38]. It leads to temporary loss of specific cell function and is involved in some clinical disorders, including migraine, cerebrovascular disease, head injury and transient global amnesia. Its mechanism and physiology have recently been reviewed by Gorji [39] and Somjen [40]. Cortical spreading depression is also accompanied by robust blood flow changes [9, 10], making it a good model for testing the feasibility of diffuse optical imaging of blood flow.

The paper is structured as follows. In Section 2, we describe the basic theory of diffuse correlation spectroscopy (DCS) and the image reconstruction algorithm we use for three-dimensional tomography of blood flow. In Section 3, we introduce a noise model for the DCS measurements and test its validity with experiment (readers not interested in the noise model can skip this section without loss of continuity). After we show explicitly how an optimal data set and regularization parameters can be obtained (Section 4), computer simulated data is generated and used to determine optimal parameters for the image reconstruction (Section 5). Our findings are then employed to reconstruct three dimensional (3D) *in-vivo* images of relative cerebral blood flow (rCBF) in a rat brain CSD model (Section 6). A discussion of the results follows these sections.

## 2. DCS theory and Image reconstruction method

In the near infrared wavelength range (∼ 650–950 nm range), the propagation of the unnormalized electric field temporal auto-correlation function, *G*
_{1}(**r**, * τ*) = 〈

*E*(

**r**,

*t*)

*E*

^{*}(

**r**,

*t*+

*))〉*

*τ**, inside tissue can be accurately described using a diffusion equation [37]:*

_{t}Here, *E*(**r**,*t*) is the electric field at position **r** and time *t*, ^{*} denotes complex conjugate, * τ* is the correlation delay time and 〈〉

*denotes time average.*

_{t}*μ*and

_{a}*μ*

_{s}^{́}are the absorption and reduced scattering coefficients of the tissue respectively,

*v*is the light speed in the media,

*D*≅

*v*/3

*μ*

_{s}^{́}is the light diffusion constant,

*k*

_{0}is the optical wave-vector,

*α*is the percentage of light scattering events from moving scatterers (eg., blood cells), 〈Δ

*r*

^{2}(

*)〉 is the mean-square displacement of the moving scatters in time*

*τ**, and*

*τ**Sδ*

^{3}(

**r**-

**r**

_{s}) is the point source term located at

**r**

_{s}. This differential equation approach [37] is formally equivalent to the original integral formulation (termed diffusing-wave spectroscopy [41,42]), but is particulary attractive for investigation of heterogeneous media [35–37]. Its solution in the semi-infinite homogeneous geometry is [43],

where *K*
^{2}(*τ*) = 3*μ _{a}*

*μ*

_{s}^{́}+

*μ*

_{s}^{́2}

*k*

^{2}

_{0}

*α*〈Δ

*r*

^{2}(

*τ*)〉,

*r*

_{1}= ∣

**r**-

**r**

_{s}+

*z*

_{tr}**n**̂∣,

*r*

_{2}= ∣

**r**-

**r**

_{s}- (

*z*+ 2

_{tr}*z*)

_{b}**n**̂∣,

**n**̂ is the unit normal to the boundary pointing away from the turbid medium,

*z*= 1/

_{tr}*μ*

_{s}^{́}is the distance into the media where the collimated source is considered isotropic,

*z*is the extrapolation distance from the sample boundary as determined by the mismatch in the indices of refraction; here we use

_{b}*z*= 2/3

_{b}*μ*

_{s}^{́}to be consistent with literature [44]. For the important case of random ballistic flow in the tissue vasculature, 〈Δ

*r*

_{2}(

*τ*)〉 =

*V*

^{2}

*τ*

^{2}. Here

*V*

_{2}is the second moment of the cell velocity distribution. For the case of diffusive motion, 〈Δ

*r*

^{2}(

*τ*)} = 6

*D*

_{b}*τ*, where

*D*is an effective diffusion coefficient of the moving scatterers. Generally, we have found that both of these models fit our

_{b}*in-vivo*data [30], but the latter model often provides better quality fits. Furthermore, we have found that relative changes in

*αD*correspond quite well to relative changes in blood flow measured by other techniques [31,32, 45–47].

_{b}In experiments, the temporal field auto-correlation function is not measured directly. Instead, the intensity fluctuations within a single speckle area are detected using a single mode fiber and a photon counting detector. A custom made correlator board then uses the output from the detector to calculate the temporal intensity auto-correlation functions of the scattered light, *g*
_{2}(**r**, *τ*) = 〈*I*(**r**,0)*I*(**r**, *τ*)〉/〈*I*(**r**, 0)〉^{2}. The normalized temporal field auto-correlation function *g*
_{1}(**r**, *τ*) = *G*
_{1}(**r**, *τ*)/*G*
_{1}(**r**, 0), in turn, is linked to *g*
_{2}(**r**, *τ*) through the Siegert relation [48]:

*β* is a parameter that depends on the detection optics and is approximately inversely proportional to the number of speckles within the detection area.

Diffuse correlation tomography uses DCS measurements from many source-detector pairs to construct a blood perfusion image. When the dynamic properties of the media are heterogeneous, the measured auto-correlation functions contain contributions from all volume elements inside the medium. Within the Rytov approximation [35], we write *g*
_{1}(**r**
_{s},**r _{d}**,

*τ*) =

*g*

_{1,0}(

**r**

_{s},

**r**,

_{d}*τ*)

*exp*(

*ϕ*(

_{s}**r**

_{s},

**r**,

_{d}*τ*)), where

*g*

_{1,0}(

**r**

_{s},

**r**,

_{d}*τ*) is the contribution of the homogeneous background, and

*ϕ*(

_{s}**r**

_{s},

**r**,

_{d}*τ*) accounts for perturbations due to the heterogeneities of dynamic and static optical properties. Then, following the procedure of Kak and Slaney [49], and assuming, for simplicity, that changes in absorption, Δ

*μ*= 0 cm

_{a}^{-1}, and scattering, Δ

*μ*

_{s}^{́}= 0 cm

^{-1}, are negligible, we can derive a matrix equation, which relates the measured perturbations in

*g*

_{1}(

**r**,

*τ*) to the heterogeneous blood flow variations, i.e.

Here *N* is the number of sample voxels, *i* is from 1 to *M*(*M* is the number of measurements), *W _{ij}* links the flow perturbation in the

*j*voxel (Δ(

^{th}*αD*(

_{b}**r**

_{j}))) to the

*i*source-detector measurement pair (

^{th}*ϕ*(

_{s}**r**

_{si},

**r**,

_{di}*τ*)).

*W*can be calculated analytically from the correlation propagation model [35,37],i.e.

_{ij}here, *H*(**r _{j}**,

**r**,

_{di}*τ*) is the Green’s function for the homogeneous correlation diffusion equation. By solving this inverse problem, the spatial distribution of the heterogeneous flow properties, Δ(

*αD*(

_{b}**r**)), can be obtained. Generally, the inverse problem is ill-posed. Therefore, in order to stabilize the image reconstruction of Δ(

*αD*(

_{b}**r**)), regularization is necessary. We use Tikhonov regularization [50] as follows:

where *T* indicates a transpose. The regularization parameter *λ* is made spatially variant to reduce source-detector artifacts at the surface plane,i.e.

where *z* is the depth of each voxel, *z _{max}* is the maximum depth in the reconstruction geometry,

*λ*

*= 10*

_{e}*λ*

*is chosen to produce even image noise as a function of depth [51]. The inverse (*

_{c}*W*∙

*W*+

_{T}*λ*∙

*I*)

^{-1}is obtained by Singular Value Decomposition (SVD). In order to achieve the best image quality,we carefully have studied the optimization of delay time

*τ*and regularization parameter

*λ*

*for the diffuse correlation problem. This is discussed in detail in Section 4.*

_{c}## 3. Noise model for DCS

In order to derive meaningful optimization schemes to guide applications of DCT, a proper estimate of experimental noise must be made. However, in contrast to the problem of diffusive waves [52–54], which measures light intensity and phase shift variations caused by propagation of photon fluence rate through tissues, the noise model for correlation experiments is non-trivial. A noise model suitable for photon correlation measurements was previously developed for the single scattering limit [55, 57]. Here, we have adapted the noise model developed by Koppel [57] for fluorescence correlation spectroscopy (FCS) in the single scattering limit and tested its feasibility in diffuse correlation experiments, wherein photons experience multiple scattering.

In a typical DCS experiment, the normalized field auto-correlation function decays approximately exponentially, i.e. *g*
_{1}(*τ*) = exp(-Γ*τ*)[56]. The experiment and experimental configuration is characterized by the decay rate Γ, the correlator bin time interval *T*, the bin index *m* corresponding to the delay time *τ*, the average number of photons 〈*n*〉 within bin time *T*(i.e. 〈*n*〉 = *IT*, where *I* is the detected photon count rate), the total averaging time *t* and the parameter *β* as described in Section 2. Following the steps in reference [57] (see Appendix), the standard deviation (σ(*τ*), noise) of the measured correlation function (*g*
_{2}(*τ*) - 1) at each delay time (*τ*) is estimated to be

$$+2{\u3008n\u3009}^{-1}\beta \left({1+e}^{-2\Gamma T}\right)+{\u3008n\u3009}^{-2}\left({1+e}^{-2\Gamma \tau}\right){]}^{\frac{1}{2}}.$$

We designed a simple experiment to test the accuracy of this noise model (Fig. 1). A single source-detector pair, separated by 1 cm, is placed into an Intralipid phantom. A long coherence length (∼ 50 m) laser provides light to the phantom. The input power is then adjusted manually using an optical attenuator in order to test the noise-model under different signal-to-noise ratio (SNR) conditions. The light was detected by a photon counting avalanche photo-diode (APD) whose output was fed into a correlator board to calculate the normalized intensity autocorrelation function *g*
_{2}(*τ*). The instrument is described in further detail in Section 6.

One hundred DCS curves were collected for each input power, and the standard deviation of the fluctuations at each *τ* was calculated and plotted (dots in Fig. 2(a)). The solid lines in Fig. 2(a) are the calculated noise using equation 8 with all the input parameters obtained from experiments: *β* was obtained from the intercept of the correlation curve; Γ was obtained by fitting the experimental data with the exponentially decaying function; the averaging time *t* = 1 s was kept constant for all measurements; photon count rates were recorded by the correlator board; the binning time interval *T* and bin index *m* were fixed on the correlator board. As shown in Fig. 2(a), the measured noise decreases as the delay time *τ* increases. The “steps” in the figure are due to the multi-tau arrangement of the correlator [58]. On our correlator board, the bin time interval is *T* = 160 ns for the first 32 channels and is doubled every 16 channels thereafter. Figure 2(a) shows that noise drops when the detected light intensity increases, as expected from equation 8. On the timescales of interest (*τ* ∈ {10^{-6} s,10^{-3} s}), the noise model provides a good estimate of the measurement noise in DCS, although it sightly overestimates the noise at large *τ* when the photon count rate is high. The failure of the noise model at large delay time has also been observed in FCS studies [59,60].

Figure 2(b) shows the signal-to-noise ratio (SNR) of the measured correlation curves, *SNR* = (*g*
_{2}(*τ*) - 1)/σ(*τ*), at different light intensities. In Fig. 2(a) we see that the DCS measurement noise decreases as the delay time *τ* increases. However, the “signal” drops even faster as *τ* increases. As a result, the signal-to-noise ratio of the DCS measurement still decreases as the delay time increases. The signal-to-noise ratio can be improved by increasing the light intensity collected by the detector, as well as increasing the averaging time *t* (data not shown), but the SNR curves will have same general form.

After studying the noise in (*g*
_{2}(*τ*) - 1), we have developed a technique to estimate the noise in the perturbation *ϕ _{s}*(

*r*,

_{si}*r*,

_{di}*τ*) for our image optimization. The noise of

*ϕ*can be derived from the relations in equations 3, 4 and 8 as:

_{s}This results in a perturbation noise that is higher at large delay time *τ*. We will use equation 9 in the following section to derive the upper limit of the normalized image noise.

## 4. Optimization Criteria

In this section we describe in detail how to optimize data selection for the image reconstruction and how to choose the optimal regularization parameter to reduce image artifacts.

#### 4.1. Choosing the optimal data set

Generally, we collect the entire correlation curve for each DCS measurement, which has ∼ 200 data points, each corresponding to a different delay time, *τ*. DCT uses DCS measurements at different source-detector pairs; usually one data point from each correlation curve is picked for the image reconstruction. Here, we introduce a parameter *n*, which defines the point where the field auto-correlation function decreases to exp(-*n*) of its initial value, i.e. *g*
_{1}(*τ*) = exp(-*n*)*g*
_{1}(0) (shown in Fig. 3). From Eq. (2), it is easy to show that *τ* and *n* have the following relationship,

For each *n*, the delay time *τ* is calculated using the above formula for each source-detector pair, and is used to calculate the elements in the weight matrix *W*. The condition number of *W*(i.e. *N _{c}* =

*ϕ*/

_{max}*ϕ*, where

_{min}*ϕ*is the largest singular value of

_{max}*W*and

*ϕ*is the smallest)is often used to characterize the weight matrix. The larger the condition number, the bigger the error amplification after the system is inverted. It can be shown that [61],

_{min}where ∥ ∥ denotes the two-norm of a vector,i.e. $\Vert {\varphi}_{s}\Vert =\sqrt{{\displaystyle {\sum}_{i}{\varphi}_{si}^{2}}}$. The upper limit of the normalized image noise (∥σ_{Δ(αDb)}∥/∥Δ(*αD _{b}*)∥) can be estimated by computing the product of the normalized measurement noise (∥σ

*∥/∥*

_{ϕs}*∥) and the condition number of the weight matrix (*

_{ϕs}*N*). Therefore, by calculating the upper limit of the normalized image noise for different

_{c}*n*, the minimal point can be determined and the set of these minimal points defines the optimal data set for image reconstruction.

#### 4.2. Choosing optimal regularization parameter

The optimization of the regularization parameter *λ*
* _{c}* is achieved by conducting a standard L-Curve analysis [62]. After a data set is chosen, a perfusion image Δ(

*αD*(

_{b}**r**)) can be reconstructed with regularization parameter

*λ*

*. The solution norm*

_{c}*η*(λ

*) = ∥Δ(*

_{c}*αD*(

_{b}**r**))∥, which provides a measure of the fluctuation in the reconstructed images, and the normalized residual norm,

*ξ*(

*λ*

*) = ∥*

_{c}*W*∙ Δ(

*αD*(

_{b}**r**)) -

*ϕ*∥/∥

_{s}*ϕ*∥, which shows the quality of the data fit, are then calculated. Generally, we want to minimize both of these values to reduce image fluctuations and obtain a good fit. If we calculate

_{s}*η*(

*λ*

*) and*

_{c}*ξ*(λ

*) for different λ*

_{c}*and plot them on the x- and y-axis, an “L” shape curve will be obtained as shown in Fig. 5. The optimized*

_{c}*λ*

*is obtained at the elbow of the L-curve as the best compromise between improved data fitting and reduced image noise (minimizing the norm of the reconstructions). In our study, the curvature at each point of the L-Curve was calculated, the maximum curvature point was found and was considered as optimal.*

_{c}## 5. Simulation Results

In this section, we use simulated data to compare the image quality using different data sets and different reconstruction schemes. The simulation geometry is on the same scale as our small animal models. Thus our conclusions can be directly used to improve image quality for our *in-vivo* small animal studies.

As shown in Fig. 4, a spherical object with radius of 0.2 cm, whose center is located at (0.4 cm, 0.2 cm, 0.4 cm), sits in a homogeneous background. The dynamic property of the object is 10% lower than the background (i.e. *αD _{b}* = 0.9 × 10

^{-8}cm

^{2}/s,

*αD*= 1 × 10

_{b}^{-8}cm

^{2}/s), while the static optical properties of the sphere and background are the same (

*μ*= 0.1 cm

_{a}^{-1},

*μ*

_{s}^{́}= 8 cm

^{-1}). Twenty-five sources and 16 detectors are placed in the z = 0 cm plane and cover a region ranging from -0.6 cm to 0.6 cm in both x and y dimensions. The analytical solution for a spherical perturbation [35, 54] is used to simulate noise-free measurement data for each source-detector pair. DCS measurement noise is then calculated based on Eq. (8) and added to the noise free data with a normal distribution.

The noise in *ϕ _{s}* is estimated using Eq. (9) and $\frac{\parallel {\sigma}_{{\varphi}_{s}}\parallel}{\parallel {\varphi}_{s}\parallel}\bullet {N}_{c}$ at different

*n*, and is plotted in Fig.5(a). The optimal data set is found at

*n*= 1, which results from a balance between the image reconstruction model and the measurement noise. In Fig. 5(b), L-Curves at different

*n*are plotted to help to choose the optimal regularization parameter

*λ*

*. The curvature at each point along the curve is calculated and the maximum curvature point is considered as the turning point of the L-curve, wihch is the optimum of*

_{c}*λ*

*. Note also that the L-Curve of the*

_{c}*n*= 1 data set was the closest to the origin, which also indicates the advantage of using a data set from

*n*= 1 for image reconstruction.

In Fig. 6, we compare the reconstructed images using data from *n* = 1 with our simulations. The reconstructed 3D images cover the region (x: -0.8 cm-0.8 cm, y: -0.8 cm-0.8 cm,z: 0 cm-0.8 cm)with 1 mm^{3} voxels. For simplicity,images located every 2 mm along the z direction are shown in the figure (from left to right). The depth for each layer is marked as the title for each column. Figure 6(a) illustrates the simulation (Sim) geometry and points to the object location from the images. Figure 6(b) shows reconstructed images using data directly from the noisy raw data (Direct Raw-data Reconstruction, DRR). The reconstructed object can be seen centered at a displaced layer (z = 0.6 cm), although many image artifacts are clearly seen in the top layers. We then smooth the simulation data by fitting it with the semi-infinite solution for the diffusion equation through the minimization of *χ*
^{2} = ∥*g*
_{2m}(*τ*) - *g*
_{2c}(*τ*)∥,where *g*
_{2m}(*τ*) and *g*
_{2c}(*τ*) are the measured and calculated intensity autocorrelation curves respectively, and we use the data from the fitted curve to reconstruct the images (Fig. 6(c)) (Smoothed Fitting Reconstruction, SFR). Image artifacts from the top layers are greatly reduced by this smoothing procedure. Moreover, the reconstructed image quality can be further improved if we use the noise information and fit the data by minimizing ${\chi}^{2}=\parallel \frac{{g}_{2m}\left(\tau \right)-{g}_{2c}\left(\tau \right)}{\sigma \left(\tau \right)}\parallel $ (Noise Fitting Reconstruction, NFR). After weighting the data points at different delay times *τ* by the correct noise, the latter fitting procedure preserves more information in the raw data, as well as effectively reduces the artifacts in the images reconstructed with data from the fitted correlation curves (Fig. 6(d)). All the images in Fig. 6(b), (c), (d) were reconstructed using Eq. (6). The differences are the data used for the image reconstruction, as described above in detail.

We compare the reconstructed images more quantitatively in Fig. 7. As has been discussed in the literature [63], the point spread function (PSF) of the diffuse optical imaging techniques is large, and the reconstructed values are usually underestimated. However, if we calculate the volume weighted rCBF for the reconstructed object, as shown in Fig. 7(a), the Noise Fitting Reconstruction (NFR) gives an object very close to the simulation. We also calculate the distance from the center of the simulation to the center of our reconstructed object (Fig. 7(b)). The center of the object reconstructed from NFR is displaced only about 1 mm from our simulation geometry and provides the best location accuracy among all three reconstruction schemes. We note that the form/parameter of the regularization can affect the location of the reconstructed object. Here, we have kept both constant throughout, However, a detailed discussion of these effects is beyond the scope of this paper. In Fig.7(c), we compare the contrast-to-noise ratio (CNR) of the reconstructed images by calculating the following [64],

Here, the region of interest (ROI) is defined as the continuous region where the rCBF change is more than $\frac{1}{2}$ of the maximum change within the reconstructed object. $\overline{r{\mathrm{CBF}}_{\mathrm{ROI}}}\phantom{\rule{.2em}{0ex}}\mathrm{and}\phantom{\rule{.2em}{0ex}}\overline{r{\mathrm{CBF}}_{\mathrm{bg}}}$ are the mean values, σ* _{ROI}* and σ

*are the fluctuations,*

_{bg}*ω*=

_{ROI}*V*/(

_{ROI}*V*+

_{ROI}*V*) and

_{bg}*ω*=

_{bg}*V*/(

_{bg}*V*+

_{ROI}*V*) are the volume weights of the ROI and the background separately. Images with high contrast-to-noise ratio are better for identifying the region of interest. For the example shown, NFR, once again, gives the highest CNR of the three reconstruction schemes. Images obtained using noise information in the fitting/smoothing process provide the best result.

_{bg}We have tried the optimization and reconstruction procedures for different simulation scales, different optode configurations and different perturbations location and values. The optimization results were found to be consistent with the findings reported here. We have also tried using multiple data points simultaneously from the correlation curve measured at each source-detector pair for the image reconstruction, but it has not as yet lead to any significant improvement in the image quality.

In summary, from computer simulations we find that the data set from *n* = 1 is optimal for the image reconstruction. We also find that use of data from the fitted correlation curves can improve image quality, and use of noise information in the fitting process gives a better fit and further improves the image quality. These findings are applicable to cases wherein a relative image is reconstructed by comparison to a secondary measurement (e.g. a baseline or a reference sample), and when reconstructing absolute images using numerical solutions from a homogeneous background as the reference.

## 6. *In-vivo* 3D flow imaging of cortical spreading depression on rat brain

A portable, relatively fast (several seconds per frame), large field of view instrument was constructed for imaging blood flow changes during cortical spreading depression (CSD) in the rat [28, 30]. A non-contact probe with a grid-like pattern of source/detector fibers was developed for this purpose (Fig. 8). The probe was mounted on the film-plane of a regular 35 mm camera body, which acted as a light sealed, robust box to hold the lens system and probe. The depth-of-focus of the camera lens reduces the motion artifacts along the optical axis of the lens. This non-contact probe enabled manipulations to the animal without movement of the probe, thereby avoiding some common experimental artifacts. Furthermore, crossed-polarizers (OFR Inc., NJ, USA) were used to reduce surface reflections from the tissue. Light from a continuous wave, long coherence length laser source (Model TC40, SDL Inc., San Jose, CA, USA) operating at 800 nm was coupled to the non-contact probe through optical switches (Di-Con, CA, USA) in order to time-share the source positions. Eight fast, photon-counting APDs (SPCMAQR-14, Perkin-Elmer, Canada) with low dark current were used in parallel as DCS detection units. The output of the APDs were fed into a custom built, fast, 9-channel correlator board (Correlator.com, Bridgewater, NJ, USA) to calculate the intensity auto-correlation function *g*
_{2}(*τ*). The whole system was automated and controlled by a desktop computer. Three source and eight detector positions were used for DCT measurements and a full frame was acquired every ∼ 6.5 seconds.

Adult male Sprague-Dawley rats weighing 300-325 g were fasted overnight with free access to water. They were anesthetized with a 1 – 1.5% halothane in a 70% nitrous oxide, 30% oxygen mixture. Catheters were placed into the femoral artery for monitoring of arterial blood pressure. Body temperature was maintained at 37±0.5*°*C by a controlled heating pad. The animals were tracheotomized, mechanically ventilated, and the head was fixed on a custom stereotaxic frame. Blood gases were obtained frequently and the respirator was adjusted in order to keep the blood gases within the normal physiological range. The scalp was retracted to avoid additional complications due to the fur. A 2 mm burr-hole was made over the frontal cortex of the right hemisphere leaving the dura intact. CSD was evoked by placing a 1 mm^{3} filter paper soaked in 2 mol/L potassium chloride (KCl) onto the dura for the duration of the desired induction of CSD waves (∼ 30 min). The paper was changed rapidly every 15 minutes. The setup is illustrated in Fig. 9. After measuring 5 – 10 CSD waves, the KCl was removed and the brain was washed with saline. A total of six animals were studied and gave similar results, but we present here reconstructed images from one representative animal.

In order to reconstruct rCBF images during CSD, baseline measurements from each source-detector pair were used as references to calculate the perturbations using Eq. (4). Background optical properties were kept constant as *μ _{a}* = 0.1 cm

^{-1}and

*μ*

_{s}^{́}= 15 cm

^{-1}, while average

*αD*= 4.5 × 10

_{b}^{-8}cm

^{-2}/s was calculated from baseline measurements at different source-detector pairs as the background dynamic property (see Section 7 for the discussion about the influence of

*μ*,

_{a}*μ*

_{s}^{́}changes during CSD to our image reconstructions). After the weight matrix

*W*was built, reconstructed flow images were optimized following the descriptions in Sections 4 and 5. Since blood flow changes during CSD were large, reconstructed images using the linear Rytov approximation usually underestimated flow values, but the position mappings should be accurate [65]. Thus, we scaled the rCBF images obtained directly from the reconstruction to match the bulk rCBF changes obtained from spectroscopy measurements over the matching brain regions.

Figure 10 shows the reconstructed rCBF images in different layers of the rat brain during CSD. The burr-hole is located at the top, slightly to the left of the midline (∼ x = 0 mm). Each panel shown is averaged over 0.5 mm in the depth direction (z) starting with the top of the skull (z = 0 mm), and at 1 mm, 2 mm and 3 mm from top to bottom, respectively. Images (from left to right) are shown roughly every 20 seconds from immediately before KCl was applied (t ≈ 26 s) until the end of the first CSD peak. The panel titles indicate the corresponding time point in this figure. Along the surface of the cortex (∼ 1 mm deep), a strong increase in blood flow appears from the top and proceeds to the bottom. It is quite significant that minimal activity is visible in the top panel (which corresponds mostly to the skull) and in the bottom panel (which penetrates below the cortex).

The time series images for the layer at 1 mm depth illustrates the spreading of two CSD waves in the cortex (Fig. 11. A movie showing the rCBF changes at different brain layers during CSD is also provided as supporting media.). After the first peak, there is a sustained decrease in blood flow (∼ 3 min) which covers most of the image area. The sustained decrease has been observed previously [39] and is compatible with inhibition of neuronal activity. Three regions of interest were selected and the rCBF changes within them are plotted in Fig. 12(a). The propagation of the CSD waves can be clearly identified from the delay between each curve. Figure 12(b) shows the dependence of maximal rCBF changes on depth using the data from the second region of interest (ROI-2) as in Fig. 11. The maximal change occurs at 1 mm (just below the skull) which corresponds to the surface of the cortex. The peak spreads ∼ 0.5 mm above and below the cortical surface as expected from the broadening due to the diffuse nature of photons. There is no significant change at the surface (z = 0 mm) and in the deep region (z = 3 mm). Note that the *in-vivo* images appear to have less “cross-contamination” across layers at different depths compared to the simulation results. We believe this effect is real and is due to the localization of CSD blood flow responses to the thin two-dimensional layer of rat brain cortex. Clearly, three-dimensional tomographic *in-vivo* relative blood flow information is revealed.

## 7. Discussion

The measurement noise of DCS depends on only a few parameters which can be obtained experimentally, and therefore can be estimated based on the model discussed in Section 3. Furthermore, the upper limit of the reconstructed image noise in DCT can be estimated using the DCS noise information at each source-detector pair. This, in turn, can be minimized by correct selection of the data set for the image reconstruction. We have parameterized the autocorrelation curve (parameter n, see Section 4, we note here *n* does not have to be integer) and *n* = 1, corresponding to *g*
_{1}(*τ*) = *exp*(-1), was found to be the optimal data point for the image reconstruction. Previously, we reported that the condition number of the weight matrix decreased as we increase the parameter *n* [66]. By contrast, the noise model teaches that a data set derived from small *n*(e.g. *n* = 0.25), where the noise in *ϕ _{s}* is low, does not improve image quality because the condition number of the weight matrix is large, i.e. the measurement noise is amplified to intolerable levels when the weight matrix is inverted; however, using a data set derived from large

*n*, where the condition number is small, does not help either, because the noise in the Rytov scattered data (

*ϕ*) is increased in this regime. The interplay of these two effects is balanced by using a data set from

_{s}*n*= 1, i.e. both the condition number and the measurement noise in

*ϕ*are optimized. The development of this model was a key theoretical aspect of the present work.

_{s}The noise model enabled us to calculate the theoretical signal-to-noise-ratio (SNR) of the DCS measurement, i.e. (*g*
_{2}(*τ*) - 1)/σ(*τ*), which in turn provides practical guidelines for experimental design. In contrast to diffuse optical near-infrared spectroscopy (NIRS) wherein fiber bundles can be used to increase detection area, diffuse correlation spectroscopy typically employs single mode fibers with diameters ∼ 6 *μ*m to detect intensity fluctuations within a single speckle area. This limits the amount of light that can be collected at large separations (e.g. > 3 cm) by the flow detectors. (Although it is possible to employ multiple detectors at nearly the same position to improve signal-to-noise ratio.) As a result, the SNR of the DCS measurement is greatly reduced when the light intensity is low, as demonstrated in Fig. 2(b). Recently, the use of few-mode fibers for DCS detection have been reported [34]; in this case more light can be collected from a few speckles. However, the value of *β*(see Section 2, equation 3)is inversely proportional to the number of speckles within the detection area, and so would decrease by the same amount, and the product *β*〈*n*〉 is about the same as if a single-mode fiber was used. Therefore the noise model suggests that the SNR of the measurement would not improve when multi-mode/few-mode fibers are used. This observation was confirmed experimentally (data not shown). Ultimately, in order to increase the signal-to-noise ratio, it is necessary to adjust the averaging time, the input light power, and/or the number of fibers/detectors working in parallel.

The development of the noise model also proved to be useful in fitting the DCS data. The best fitting curves are derived when an error can be assigned to each data point from the collected DCS curve. By weighting all the data points with the appropriate estimated measurement noise in the definition of *χ*
^{2}[67],e.g. ${\chi}^{2}=\parallel \frac{{g}_{2m}\left(\tau \right)-{g}_{2c}\left(\tau \right)}{\sigma \left(\tau \right)}\parallel ,$ each data point from the curve is better used in the fitting. We have shown that if the estimated noise information is considered, the reconstructed image quality is thus improved (see Section 5). From a phantom study, we also observed a smaller standard deviation in fitted *αD _{b}*, when the estimated noise was used in the fitting (data not shown), especially when the measured DCS curves were noisy. This will be important when using more complex models to fit spectroscopy data, for example, in brain measurements where multiple layers should be considered [34,43,46].

Finally we examine some of the assumptions used for rCBF diffuse correlation tomography image reconstructions. One of our assumptions was that the static optical properties (*μ _{a}*,

*μ*

_{s}^{́}) do not change during activation. This assumption may be incorrect for

*in-vivo*animal studies. Kohl

*et al*reported

*in-vivo*dynamic oxygenation and scattering changes during cortical spreading depression [68], for example, and found that the magnitude of the optical property changes were relatively small (ΔHbO

_{2}∼ +15

*μ*M, ΔHb ∼ -7

*μ*M and Δ

*μ*

_{s}^{́}∼ 1 cm

^{-1}). Compared to the baseline brain optical properties we measured and used in our image reconstruction, these relative changes in both

*μ*and

_{a}*μ*

_{s}^{́}are less than 7%. To this end, we changed the global optical properties by the same amount during CSD, and no significant changes in the reconstructed rCBF images were observed (i.e. changes in the image voxels were less than 10%). In practice, it is desirable to carry out frequency domain diffuse optical tomography measurements concurrently with the diffuse correlation tomography measurement for

*in-vivo*studies. It will then be possible to reconstruct absorption and scattering images from DOT first, and use them for calculating DCT Green’s functions, thereby reducing the optical property influence on blood flow image reconstructions. Furthermore, by combining the DOT and DCT images, it is possible to image the cerebral metabolism rate of oxygen (CMRO

_{2}) in three dimensions [28, 69].

Over the past thirty years, there has been great interest in measuring cerebral blood flow, oxygen consumption and metabolic responses during cortical spreading depression [9, 10, 70–73]. The optical imaging technique we describe in this paper provides reliable three-dimensional *in-vivo* images of rCBF during CSD with a relatively fast frame rate (∼ 0.15 Hz) and moderate transverse and depth resolution (∼ 0.5 mm). The relative blood flow changes we observed during CSD, an initial strong increase followed by a sustained decrease, is consistent with observations from other techniques such as laser Doppler flowmetry (1D) [70], scanning laser Doppler imaging [10], and laser speckle imaging (2D) [9]. The strong increase of rCBF is believed to be coupled to the increase of oxygen consumption [71,74], wihch in turn can also be directly measured using the diffuse optical imaging methods as discussed above. On the other hand, in addition to its capability for providing three dimensional blood flow images, the non-invasive nature of our technique is desirable for studying brain activity *in-vivo*. We have recently extended this technique for local measurements of rCBF in adult human brain [29] and it is conceivable to adapt methods developed in this paper for regional imaging of relative blood flow in a variety of human tissues.

## 8. Conclusion

In summary, we have demonstrated tomographic three-dimensional relative blood flow images using diffuse correlation measurements. A noise model for DCS measurements was introduced and its accuracy in the multiple scattering regime was studied by experiment and simulation. Optimized data sets and regularization parameters for the image reconstruction were derived, and the optimal data set was achieved at time-points wherein the field auto-correlation function *g*
_{1}(*τ*) decreases to 1/*e* of its initial value. Our findings were then employed in the study of the cortical spreading depression in a rat brain, and three-dimensional *in-vivo* blood flow images during CSD were obtained using the diffuse optical correlation tomography technique.

## Appendix

In this section, we briefly outline the derivation of the noise model for DCS measurements as in Eq. (8) following the steps in [57].

In reality, the measured intensity auto-correlation function *G*̂_{2}(*τ*) is calculated as

where *τ* is the delay time, *T* is the correlator bin time interval, *n*(*iT*) is the number of photon counts in time *iT*, *N* is the total number of the measurements (exposure time *t* = *NT*). Herein, we will use 〉 〈 to represent the true value of a parameter, and ̂ to denote the experimental estimate of the true value. We begin by defining a parameter

where *n*̂= ∑^{N}_{i=1}
*n*(*iT*), and

$$\phantom{\rule{.8em}{0ex}}\approx {\u3008n\u3009}^{2}\left(1-2\frac{\u3008n\u3009-\hat{n}}{\u3008n\u3009}\right)$$

$$\phantom{\rule{.8em}{0ex}}=2\u3008n\u3009\hat{n}-{\u3008n\u3009}^{2}.$$

The noise of the normalized intensity autocorrelation function (*g*
_{2}(*τ*) - 1) can be obtained by normalizing the variance of *S*̂(*τ*), as

Using Eqs. (13), (14), and (15), the variance of *S*̂(*τ*) can be written as

$$\phantom{\rule{3.4em}{0ex}}=\mathrm{var}(\frac{1}{N}\sum _{i=1}^{N}n\left(\mathrm{iT}\right)[n\left(\mathrm{iT}+\tau \right)-2\u3008n\u3009]).$$

If we define

Eq. (17) can be expanded as [75]

Using the fact in photon statistics [76],

where *n*(*iT*) and *I*(*iT*) are the photon count and the light intensity in time *iT* seperately, and writing the higher order correlation functions as sums of products of the second order correlation functions [77], an analytical expression of the variance of *S*̂(*τ*) can be derived for exponentially decayed (*g*
_{2}(*τ*) - 1) = *β*
*e*
^{-2Γr} after tedious math

$$+2{\u3008n\u3009}^{3}\beta \left(1+{e}^{-2\Gamma \tau}\right)+{\u3008n\u3009}^{2}\left(1+{\beta e}^{-\Gamma \tau}\right)].$$

Consequently, the noise (σ(*τ*)) of the normalized intensity autocorrelation function (*g*
_{2}(*τ*)-1) can be expressed as

$${+2{\u3008n\u3009}^{-1}\beta \left(1+{e}^{-2\Gamma \tau}\right)+{\u3008n\u3009}^{-2}\left(1+{\beta e}^{-\Gamma \tau}\right)]}^{\frac{1}{2}},$$

as shown in Eq. (8). And in the limit of Γ*T* ≤ 1,which is satisfied in most of our experiments,

## Acknowledgments

We thank Alper Corlu, Goro Nishimura, Joseph P. Culver, David A. Boas and Douglas J. Durian for useful discussions and gratefully acknowledge funding from the National Institute of Health grants 2–RO1–HL–57835–04, and NS-033785.

## References and links

**1. **A. Zauner and J. P. Muizelaar, *Head Injury, Chapter 11* (Chapman and Hall, 1997).

**2. **R. S. J. Frackowiak, G. L. Lenzi, T. Jones, and J. D. Heather, “Quantitative Measurement of Regional Cerebral Blood-Flow and Oxygen-Metabolism in Man Using O-15 and Positron Emission Tomography - Theory, Procedure, and Normal Values,” J. Comput. Assist. Tomogr. **4**, 727–736 (1980). [CrossRef] [PubMed]

**3. **D. S. Williams, J. A. Detre, J. S. Leigh, and A. P. Koretsky, “Magnetic resonance imaging of perfusion using spin inversion of arterial water.” Proc. Natl. Acad. Sci. U. S. A. **89**, 212–216 (1992). [CrossRef] [PubMed]

**4. **J. A. Detre and D. C. Alsop, “Perfusion magnetic resonance imaging with continuous arterial spin labeling: methods and clinical applications in the central nervous system,” Eur. J. Radiol. **30**, 115–124 (1999). [CrossRef] [PubMed]

**5. **G. Zaharchuk, J. Bogdanov, A. A., J. J. Marota, M. Shimizu-Sasamata, R. M. Weisskoff, K. K. Kwong, B. G. Jenkins, R. Weissleder, and B. R. Rosen, “Continuous assessment of perfusion by tagging including volume and water extraction (CAPTIVE): a steady-state contrast agent technique for measuring blood flow, relative blood volume fraction, and the water extraction fraction,” Magn. Reson. Med. **40**, 666–678. (1998). [CrossRef] [PubMed]

**6. **A. K. Dunn, H. Bolay, M. A. Moskowitz, and D. A. Boas, “Dynamic imaging of cerebral blood flow using laser speckle,” J. Cereb. Blood Flow Metab. **21**, 195–201 (2001). [CrossRef] [PubMed]

**7. **A. K. Dunn, A. Devor, H. Bolay, M. L. Andermann, M. A. Moskowitz, A. M. Dale, and D. A. Boas, “Simultaneous imaging of total cerebral hemoglobin concentration, oxygenation, and blood flow during functional activation,” Opt. Lett. **28**, 28–30 (2003). [CrossRef] [PubMed]

**8. **T. Durduran, M. G. Burnett, G. Yu, C. Zhou, D. Furuya, A. G. Yodh, J. A. Detre, and J. H. Greenberg, “Spa-tiotemporal Quantification of Cerebral Blood Flow During Functional Activation in Rat Somatosensory Cortex Using Laser-Speckle Flowmetry,” J. Cereb. Blood Flow Metab. **24**, 518–525 (2004). [CrossRef] [PubMed]

**9. **C. Ayata, H. K. Shin, S. Salomone, Y. Ozdemir-Gursoy, D. A. Boas, A. K. Dunn, and M. A. Moskowitz, “Pronounced hypoperfusion during spreading depression in mouse cortex,” J. Cereb. Blood Flow Metab. **24**, 1172–1182 (2004). [CrossRef] [PubMed]

**10. **A. N. Nielsen, M. Fabricius, and M. Lauritzen, “Scanning laser-Doppler flowmetry of rat cerebral circulation during cortical spreading depression,” J. Vasc. Res. **37**, 513–522 (2000). [CrossRef]

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

**12. **A. G. Yodh and D. A. Boas, *Biomedical Photonics* (CRC Press, 2003). Chapter Functional Imaging with Diffusing Light.

**13. **D. A. Boas, M. A. Franceschini, A. K. Dunn, and G. Strangman, “Non-Invasive imaging of cerebral activation with diffuse optical tomography,” in *Optical Imaging of Brain Function*,
R. Frostig, ed. (CRC Press, 2002). [CrossRef]

**14. **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]

**15. **A. H. Hielscher, “Optical tomographic imaging of small animals,” Curr. Opin. Biotechnol. **16**, 79–88 (2005). [CrossRef] [PubMed]

**16. **A. Villringer and B. Chance, “Non-invasive optical spectroscopy and imaging of human brain function,” Trends Neurosci. **20**, 435–442 (1997). [CrossRef] [PubMed]

**17. **G. Gratton, M. Fabiani, P. M. Corballis, D. C. Hood, M. R. Goodman-Wood, J. Hirsch, K. Kim, D. Friedman, and E. Gratton, “Fast and localized event-related optical signals (EROS) in the human occipital cortex: comparisons with the visual evoked potential and fMRI.” Neuroimage **6**, 168–180 (1997). [CrossRef] [PubMed]

**18. **B. W. Pogue and K. D. Paulsen, “High-resolution near-infrared tomographic imaging simulations of the rat cranium by use of apriori magnetic resonance imaging structural information,” Opt. Lett. **23**, 1716–1718 (1998). [CrossRef]

**19. **D. A. Benaron, S. R. Hintz, A. Villringer, D. Boas, A. Kleinschmidt, J. Frahm, C. Hirth, H. Obrig, J. C. van Houten, E. L. Kermit, W. F. Cheong, and D. K. Stevenson, “Noninvasive functional imaging of human brain using light,” J. Cereb. Blood Flow Metab. **20**, 469–477 (2000). [CrossRef] [PubMed]

**20. **D. A. Boas, D. H. Brooks, E. L. Miller, C. A. DiMarzio, M. Kilmer, R. J. Gaudette, and Q. Zhang, “Imaging the body with diffuse optical tomography,” IEEE Signal Process. Mag. **18**, 57–75(2001). [CrossRef]

**21. **D. M. Hueber, M. A. Franceschini, H. Y. Ma, Q. Zhang, J. R. Ballesteros, S. Fantini, D. Wallace, V. Ntziachristos, and B. Chance, “Non-invasive and quantitative near-infrared haemoglobin spectrometry in the piglet brain during hypoxic stress, using a frequency-domain multidistance instrument,” Phys. Med. Biol. **46**, 41–62. (2001). [CrossRef] [PubMed]

**22. **A. Bluestone, G. Abdoulaev, C. Schmitz, R. Barbour, and A. Hielscher, “Three-dimensional optical tomography of hemodynamics in the human head,” Opt. Express **9**, 272–286 (2001). [CrossRef] [PubMed]

**23. **J. C. Hebden, A. Gibson, R. M. Yusof, N. Everdell, E. M. C. Hillman, D. T. Delpy, S. R. Arridge, T. Austin, J. H. Meek, and J. S. Wyatt, “Three-dimensional optical tomography of the premature infant brain,” Phys. Med. Biol. **47**, 4155–4166 (2002). [CrossRef] [PubMed]

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

**25. **T. Wilcox, H. Bortfeld, R. Woods, E. Wruck, and D. A. Boas, “Using near-infrared spectroscopy to assess neural activation during object processing in infants,” J. Biomed. Opt. **10**, 011,010 (2005). [CrossRef] [PubMed]

**26. **E. Gratton, V. Toronov, U. Wolf, M. Wolf, and A. Webb, “Measurement of brain activity by near-infrared light,” J. Biomed. Opt. **10**, 011,008 (2005). [CrossRef] [PubMed]

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

**28. **J. P. Culver, T. Durduran, T. 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]

**29. **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]

**30. **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]

**31. **G. Q. Yu, T. Durduran, C. Zhou, H. W. Wang, M. E. Putt, H. M. Saunders, C. M. Sehgal, E. Glatstein, A. G. Yodh, and T. M. Busch, “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]

**32. **G. Q. Yu, T. Durduran, G. Lech, C. Zhou, B. Chance, R. E. Mohler, and A. G. Yodh, “Time-dependent blood flow and oxygenation in human skeletal muscles measured with noninvasive near-infrared diffuse optical spec-troscopies,” J. Biomed. Opt. **10**, 024,027-1-12 (2005). [CrossRef] [PubMed]

**33. **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]

**34. **J. Li, G. Dietsche, D. Iftime, S. E. Skipetrov, G. Maret, T. Elbert, B. Rockstroh, and T. Gisler, “Noninvasive detection of functional brain activity with near-infrared diffusing-wave spectroscopy,” J. Biomed. Opt. **10**, 1–12 (2005). [CrossRef]

**35. **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-Opt. Image Sci. Vis. **14**, 192–215 (1997). [CrossRef]

**36. **M. Heckmeier, S. E. Skipetrov, G. Maret, and R. Maynard, “Imaging of dynamic heterogeneities in multiple-scattering media,” J. Opt. Soc. Am. A-Opt. Image Sci. Vis. **14**, 185–191 (1997). [CrossRef]

**37. **D. A. Boas, L. E. Campbell, and A. G. Yodh, “Scattering and Imaging with Diffusing Temporal Field Correla-tions,” Phys. Rev. Lett. **75**, 1855–1858 (1995). [CrossRef] [PubMed]

**38. **A. A. P. Leao, “Spreading depression of activity in cerebralcortex,” J. Neurophysiol. **7**, 359–390 (1944).

**39. **A. Gorji, “Spreading depression: a review of the clinical relevance,” Brain Res. Rev. **38**, 33–60 (2001). [CrossRef] [PubMed]

**40. **G. G. Somjen, “Mechanisms of spreading depression and hypoxic spreading depression-like depolarization,” Physiol. Rev. **81**, 1065–1096 (2001). [PubMed]

**41. **G. Maret and P. Wolf, “Multiple light scattering from disordered media. The effect of brownian motion of scat terers,” Z. Phys. B. **65**, 409–413 (1987). [CrossRef]

**42. **D. Pine, D. Weitz, P. Chaikin, and Herbolzheimer, “Diffusing-wave spectroscopy,” Phys. Rev. Lett. **60**, 1134–1137 (1988). [CrossRef] [PubMed]

**43. **D. Boas, “Diffuse PhotonProbes of StructuralandDynamicalProperties of Turbid Media: Theory and Biomed-ical Applications,” *Ph.D.,*University of Pennsylvania (1996).

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

**45. **C. Menon, G. M. Polin, I. Prabakaran, A. Hsi, C. Cheung, J. P. Culver, J. F. Pingpank, C. S. Sehgal, A. G. Yodh, D. G. Buerk, and D. L. Fraker, “An integrated approach to measuring tumor oxygen status using human melanoma xenografts as a model,” Cancer Res. **63**, 7232–7240 (2003). [PubMed]

**46. **T. Durduran, “Noninvasive measurements of tissue hemodynamics with hybrid diffuse opticalmethods,” *Ph.D.,*University of Pennsylvania (2004).

**47. **G. Yu, T. F. Floyd, T. Durduran, C. Zhou, J. J. Wang, J. M. Murphy, and A. G. Yodh, “Concurrent Optical-MRI Measurement of Limb Blood Flow/Perfusion,” Opt. Lett. in prep (2005). [PubMed]

**48. **S. Rice, “Mathematical analysis of random noise,” in *Noise and Stochastic Processes*,
N. Wax, ed., p. 133 (Dover, NewYork, 1954).

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

**50. **S. R. Arridge, “OpticalTomography in medicalimaging,” Inverse Probl. **15**, R41–R93 (1999). [CrossRef]

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

**52. **M. A. Oleary, D. A. Boas, B. Chance, and A. G. Yodh, “Refraction of diffuse photon density waves,” Phys. Rev. Lett. **69**, 2658–2661 (1992). [CrossRef]

**53. **D. A. Boas, M. A. Oleary, B. Chance, and A. G. Yodh, “Scattering and Wavelength Transduction of Diffuse Photon Density Waves,” Phys. Rev. E **47**, R2999–R3002 (1993). [CrossRef]

**54. **D. A. Boas, M. A. Oleary, 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. U. S. A. **91**, 4887–4891 (1994). [CrossRef] [PubMed]

**55. **K. Schatzel, “Noise in photon-correlation and photon structure functions,” Optica ACTA **30**, 155–166 (1983). [CrossRef]

**56. **
The solution to the correlationdiffusion equation,i.e., Eq.(2), is amore accurate description forg *g*_{1}(*τ*). However, when the delay time *τ* is small $\left(\tau \ll \frac{{3}_{{\mu}_{a}}}{{\mu}_{s}^{\u0301}{k}_{0}^{2}\alpha {D}_{b}}\right),$*g*_{1}(*t*) can be simplified as an exponential decay function. On the other hand, we have compared the noise calculated from Eq. (8) assuming exponential decay and the noise calculated numerically using exact the semi-infinite solution as input. No significant difference was observed.

**57. **D. E. Koppel, “Statistical accuracy influorescence correlation spectroscopy,” Phys. Rev. A**10**, 1938–1945 (1974). [CrossRef]

**58. **K. Schatzel, M. Drewel, and S. Stimac, “Photon-Correlation Measurements at Large Lag Times - Improving Statistical Accuracy,” J. Mod. Opt. **35**, 711–718 (1988). [CrossRef]

**59. **U. Meseth, T. Wohland, R. Rigler, and H. Vogel, “Resolution of fluorescence correlation measurements,” Bio-phys. J. **76**, 1619–1631 (1999).

**60. **T. Wohland, R. Rigler, and H. Vogel, “The standard de viationinfluorescence correlation spectroscopy,” Biophys. J. **80**, 2987–2999 (2001). [CrossRef] [PubMed]

**61. **S. H. Friedberg, A. J. Insel, and L. E. Spence, *Linear algebra*, 3rd ed. (Prentice Hall, 1997).

**62. **P. Hansen, “Analysis of discreteill-posed problems by means of the L-curve,” SIAM Rev. **34**, 561–580 (1992). [CrossRef]

**63. **J. P. Culver, R. Choe, M. J. Holboke, L. Zubkov, T. Durduran, A. Slemp, V. Ntziachristos, B. Chance, and A. G. Yodh, “Three-dimensionaldiffuse opticaltomography in the parallelplane transmission geometry: Evaluation of a hybrid frequency domain/continuous wave clinical system for breastimaging,” Med. Phys. **30**, 235–247 (2003). [CrossRef] [PubMed]

**64. **X. M. Song, B. W. Pogue, S. D. Jiang, M. M. Doyley, H. Dehghani, T. D. Tosteson, and K. D. Paulsen, “Automated region detection based on the contrastto-noise ratio in near-infrared tomography,” Appl Optics **43**, 1053–1062 (2004). [CrossRef]

**65. **M. A. O’Leary, “Imaging with diffuse photon density waves,” *Ph.D.,*Unversity ofPennsylvania (1996).

**66. **C. Zhou, T. Durduran, G. Yu, and A. G. Yodh, “Optimizing Image reconstruction of tissue blood flow by diffuse correlation tomography,” in *Photonics West, SPIE*, vol. 4955-43, pp. 287-95 (San Jose, CA, 2003).

**67. **P. E. Greenwood and M. S. Nikulin, *A guide to chi-squared testing*, Wiley series in probability and statistics. Applied probab lity and statistics (New York: John Wiley & Sons, 1996).

**68. **M. Kohl, U. Lindauer, U. Dirnagl, and A. Villringer, “Separation of changes in light scattering and chromophore concentrations during cortical spreading depression in rats,” Opt. Lett. **23**, 555–557 (1998). [CrossRef]

**69. **R. D. Hoge, J. Atkinson, B. Gill, G. R. Crelier, S. Marrett, and G. B. Pike, “Investigation of BOLD signal dependence on cerebral blood flow and oxygen consumption: the deoxyhemoglobin dilution model,” Magn. Reson. Med. **42**, 849–63 (1999). [CrossRef] [PubMed]

**70. **M. Lauritzen and M. Fabricius, “Real time laser-Doppler perfusion imaging of cortical spreading depression in rat neocortex,” Neuroreport **6**, 1271–1273 (1995). [CrossRef] [PubMed]

**71. **A. Mayevsky and H. R. Weiss, “Cerebral Blood-Flow and Oxygen-Consumption in Cortical Spreading Depression,” J. Cereb. Blood Flow Metab. **11**, 829–836 (1991). [CrossRef] [PubMed]

**72. **A. Mayevsky and B. Chance, “Repetitive Patterns of Metabolic Changes During Cortical Spreading Depression of Awake Rat,” Brain Res. **65**, 529–533 (1974). [CrossRef] [PubMed]

**73. **J. Sonn and A. Mayevsky, “Effects of brain oxygenation on metabolic, hemodynamic, ionic and electrical responses to spreading depression in the rat,” Brain Res. **882**, 212–216 (2000). [CrossRef] [PubMed]

**74. **L. D. Lukyanov and J. Bures, “Changes in PO2 Due to Spreading Depression in Cortex and Nucleus Caudatus of Rat,” Physiologia Bohemoslovaca **16**, 449–455 (1967).

**75. **W. B. Davenport and W. L. Root, *Random Signals and Noise* (McGraw-Hill, 1958).

**76. **C. D. Cantrell, “N-Fold Photonelectric counting statistics of gaussian light,” Phys. Rev. A **1**, 672–685 (1970). [CrossRef]

**77. **P. A. Lemieux and D. J. Durian, “Investigating non-Gaussian scattering processes by using nth-order intensity correlation functions,” J. Opt. Soc. Am. A-Opt. Image Sci. Vis. **16**, 1651–1664 (1999). [CrossRef]