Open Access
Add to BibSonomyAbstract
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 + τ))〉t, inside tissue can be accurately described using a diffusion equation [37]:
Here, E(r,t) is the electric field at position r and time t, * denotes complex conjugate, τ is the correlation delay time and 〈〉t denotes time average. μa and μ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 + ztr n̂∣, r 2 = ∣r - r s - (ztr + 2zb)n̂∣, n̂ is the unit normal to the boundary pointing away from the turbid medium, ztr = 1/μs ́ is the distance into the media where the collimated source is considered isotropic, zb is the extrapolation distance from the sample boundary as determined by the mismatch in the indices of refraction; here we use zb = 2/3μ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(τ)} = 6Db τ, where Db is an effective diffusion coefficient of the moving scatterers. Generally, we have found that both of these models fit our in-vivo data [30], but the latter model often provides better quality fits. Furthermore, we have found that relative changes in αDb correspond quite well to relative changes in blood flow measured by other techniques [31,32, 45–47].
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,rd,τ) = g 1,0(r s,rd,τ) exp(ϕs(r s,rd,τ)), where g 1,0(r s,rd,τ) is the contribution of the homogeneous background, and ϕs(r s,rd,τ) 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, Δμa = 0 cm-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), Wij links the flow perturbation in the jthvoxel (Δ(αDb(r j))) to the ith source-detector measurement pair (ϕs(r si, rdi, τ)). Wij can be calculated analytically from the correlation propagation model [35,37],i.e.
here, H(rj, rdi, τ) is the Green’s function for the homogeneous correlation diffusion equation. By solving this inverse problem, the spatial distribution of the heterogeneous flow properties, Δ(αDb(r)), can be obtained. Generally, the inverse problem is ill-posed. Therefore, in order to stabilize the image reconstruction of Δ(αDb(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, zmax is the maximum depth in the reconstruction geometry, λ e = 10λ c is chosen to produce even image noise as a function of depth [51]. The inverse (W ∙ WT + λ ∙ 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 λ c for the diffuse correlation problem. This is discussed in detail in Section 4.
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.
Fig. 1. Experimental setup used to test the accuracy of the noise model. One source-detector pair, with a 1 cm separation, is placed into an Intralipid phantom. A long coherence length laser (∼ 50 m) is provides light to the phantom. In order to test the noise-model under different signal-to-noise ratio conditions, the input power is adjusted manually using an optical attenuator connected to the input fiber. The light is detected by a photon counting APD the output of which is fed into a correlator board to calculate the normalized intensity auto-correlation function g 2(τ). g 2(τ) is then collected and saved in a desktop computer. A hundred g 2(τ) curves are measured under the same conditions. The measurement noise, plotted in Fig. 2(a), is calculated as the standard deviation of the fluctuation at each delay time τ.
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
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].
Fig. 2. (a) Comparison of measured noise (dots) and calculated noises using the model (solid lines). All input parameters for the noise model are obtained from experiments. Measurement noise decreases as the delay time τ increases. The “steps” are due to the multi-tau arrangement of our correlator. (b) Signal-to-Noise Ratio (SNR) comparison of the measured correlation curves and the model predictions. Although the measurement noise decreases as the delay time τ increases, the SNR of the DCS measurements also decreases because the “signal” drops even faster as τ increases. (kcps = kilo-counts per second)
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(rsi,rdi,τ) for our image optimization. The noise of ϕs can be derived from the relations in equations 3, 4 and 8 as:
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,
Fig. 3. Example of the field auto-correlation function g 1(τ). Points corresponding to different values of n are notated n is defined as the point where g 1(τ) decreases to exp(-n) of its initial value (i.e. we write g 1(τ) = exp(-n)g 1(0)).
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. Nc = ϕmax/ϕmin, where ϕmax is the largest singular value of W and ϕmin 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],
where ∥ ∥ denotes the two-norm of a vector,i.e. . The upper limit of the normalized image noise (∥σΔ(αDb)∥/∥Δ(αDb)∥) can be estimated by computing the product of the normalized measurement noise (∥σϕs∥/∥ϕs∥) and the condition number of the weight matrix (Nc). Therefore, by calculating the upper limit of the normalized image noise for different 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 Δ(αDb(r)) can be reconstructed with regularization parameter λ c. The solution norm η(λc) = ∥Δ(αDb(r))∥, which provides a measure of the fluctuation in the reconstructed images, and the normalized residual norm, ξ(λ c) = ∥W ∙ Δ(αDb(r)) - ϕs∥/∥ϕ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 η(λ c) 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.
Fig. 4. Simulation geometry: (a) A spherical object with a radius of 0.2 cm is placed in a homogeneous background at three distances from the source/detector plane (0.4 cm, 0.2 cm, 0.4 cm). The dynamic property of the object is 10% lower than the background (αDb, = 0.9 × 10-8 cm2/s, αDb0 = 1 × 10-8 cm2/s). The static optical properties of the sphere and background are the same (μa = 0.1 cm-1, μs ́ = 8 cm-1). (b) 25 sources and 16 detectors are placed at the z = 0 cm plane and cover a region ranging from -0.6 cm to 0.6 cm in both the x and y dimensions. An analytical solution is used for the simulation. Measurement noise is calculated and added to the simulated data with a normal distribution.
Fig. 5. Choice of the optimal data set and optimal regularization parameter: (a) The normalzed image noise plotted as a function of n. The optimal data set is obtained at n = 1, where the upper limit of the normalized image noise is a minimum. (b) L-Curves with noisy data at different n are plotted to identify the optimal regularization parameter λ. The n = 1 curve is the closest to the origin, which also indicates the advantage of using a data set with n = 1 for image reconstruction.
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. αDb = 0.9 × 10-8cm2/s, αDb = 1 × 10-8 cm2/s), while the static optical properties of the sphere and background are the same (μa = 0.1 cm-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.
Fig. 6. Reconstructed images using data with n = 1. 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 mm3 voxels. Images at every 2 mm along the z direction are shown (from left to right). The depth for each layer is marked for each column. (a) Simulation (Sim) geometry. (b) Reconstructed images using data directly from the noisy raw data (Direct Raw-data Reconstruction, DRR). The object is found at a displaced layer (z=0.6 cm). (c) Images using data from the fitted curve (by minimizing ∥g 2m(τ) - g 2c(τ)∥) to reconstruct the images (Smoothed Fitting Reconstruction, SFR). Image artifacts are greatly reduced. (d) Using noise information in the fitting process (by minimizing ) can further improve the image quality (Noise Fitting Reconstruction, NFR).
Fig. 7. Quantitative comparison of different image reconstruction schemes: (a) Volume weighted rCBF for the reconstructed object normalized to the simulation (Sim). Noise Fitting Reconstruction (NFR) gives the most accurate value compared to the simulation. (b) Distance from the center of the simulation to the center of the reconstructed object. NFR provides the best location accuracy (∼ 1 mm). (c) Contrast to noise ratio (CNR) of the reconstructed images. NFR gives the highest CNR over the three reconstruction schemes.
The noise in ϕs is estimated using Eq. (9) and 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 λ c. 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 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 mm3 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 (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 of the maximum change within the reconstructed object. are the mean values, σROI and σbg are the fluctuations, ωROI = VROI /(VROI + Vbg) and ωbg = Vbg/(VROI + Vbg) 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.
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 mm3 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 αDb = 4.5 × 10-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.
Fig. 8. System setup for the in-vivo rat brain CSD study. A grid-like pattern of source/detector fibers (3 sources, 8 detectors) was mounted on the back of a 35 mm camera body. The light was sent to and detected from the tissues through a relay lens avoiding contact with the tissue. Light from a CW, long coherence length laser operating at 800 nm was connected to the non-contact probe through optical switches in order to time-share the source positions. The output of the APDs was fed into a custom built correlator board which calculated the intensity auto-correlation function g 2(τ). The whole system was automated and controlled by a desktop computer and a full frame was acquired every ∼ 6.5 seconds.
Fig. 9. Cortical Spreading Depresstion (CSD); During experiments, a rat was fixed on a stereotaxic frame with the scalp retracted and the skull intact. CSD was induced by placing KCl solution on the rat brain through a small hole drilled on the skull. Periodic activations and deactivations of the neurons then spread out radially on the cortex as shown on the sketch.
Fig. 10. Reconstructed 3D rCBF images at different layers (from top to bottom) of the rat brain during CSD. CSD is induced at the top, slightly to the left of the midline (∼ x = 0 mm). Images (from left to right) are shown about every 20 seconds from immediately before KCl was applied until the end of the first CSD peak. rCBF responses are mainly localized to the cortex and spread across the cortex tangentially from the point where KCl was applied. No significant activity is visible in the top panel, which corresponds to the skull, and in the bottom panel, which penetrates below the cortex, indicating that the activity is localized in the cortex. Images are oriented as shown in Fig. 8(b).
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).
Fig. 11. rCBF changes on the cortex of the rat brain during CSD. Images (from left to right, from top to bottom) are shown roughly every 20 seconds from immediately before KCl was applied until the end of the second CSD peak. The panel titles indicate the corresponding time point. A strong increase in blood flow appears from the top and proceeds to the bottom of the image. After the peak, there is a sustained decrease in blood flow which covers most of the image area. Three regions of interest (ROI) were selected and the time series of rCBF changes from these ROIs are plotted in Fig. 12(a). A movie showing the rCBF changes at different brain layers during CSD is also provided as supporting media (1.5MB). [Media 1]
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 (ϕs) is increased in this regime. The interplay of these two effects is balanced by using a data set from n = 1, i.e. both the condition number and the measurement noise in ϕs are optimized. The development of this model was a key theoretical aspect of the present work.
Fig. 12. (a) Time series of rCBF changes from three ROIs illustrated in Fig. 11. From these curves, the propagation of the CSD waves can be clearly identified. (b) Dependence of maximum rCBF on depth in the second region of interest. The maximal change is localized at a depth of 1 mm below the skull This 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 (see text). There is no significant change at the surface (skull) or in deeper regions.
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. 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 αDb, 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 (ΔHbO2 ∼ +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 μa and μ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 (CMRO2) 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 Ĝ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
The noise of the normalized intensity autocorrelation function (g 2(τ) - 1) can be obtained by normalizing the variance of Ŝ(τ), as
Using Eqs. (13), (14), and (15), the variance of Ŝ(τ) can be written as
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 Ŝ(τ) can be derived for exponentially decayed (g 2(τ) - 1) = β e -2Γr after tedious math
Consequently, the noise (σ(τ)) of the normalized intensity autocorrelation function (g 2(τ)-1) can be expressed as
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 g1(τ). However, when the delay time τ is small g1(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. A10, 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]
