## Abstract

We introduce an integration of dynamic light scattering (DLS) and optical coherence tomography (OCT) for high-resolution 3D imaging of heterogeneous diffusion and flow. DLS analyzes fluctuations in light scattered by particles to measure diffusion or flow of the particles, and OCT uses coherence gating to collect light only scattered from a small volume for high-resolution structural imaging. Therefore, the integration of DLS and OCT enables high-resolution 3D imaging of diffusion and flow. We derived a theory under the assumption that static and moving particles are mixed within the OCT resolution volume and the moving particles can exhibit either diffusive or translational motion. Based on this theory, we developed a fitting algorithm to estimate dynamic parameters including the axial and transverse velocities and the diffusion coefficient. We validated DLS-OCT measurements of diffusion and flow through numerical simulations and phantom experiments. As an example application, we performed DLS-OCT imaging of the living animal brain, resulting in 3D maps of the absolute and axial velocities, the diffusion coefficient, and the coefficient of determination.

© 2012 OSA

## 1. Introduction

Dynamic light scattering (DLS) is widely used to quantify the dynamics of scattering particles by analyzing the autocorrelation function of light scattered from the particles [1–3]. Optical coherence tomography (OCT) uses coherence gating to enable the measurement and analysis of light back-scattered from small volumes (typically on the order of a few wavelengths), facilitating high-resolution 3D imaging of the sample structure [4–6]. Therefore, an integration of these two techniques would provide high-resolution 3D imaging of particle dynamics (e.g., a 3D map of the diffusion coefficient quantifying the diffusive motion of particles at every voxel).

This paper describes the integration of DLS and OCT leading to a novel technique for high-resolution 3D imaging of particle dynamics, named DLS-OCT. For this integration, we propose and validate a theory that considers three concerns simultaneously. First, as OCT provides a very small probing volume, the theory considers the effect of such a wavelength-scale probing volume on the autocorrelation function of scattered light. Second, the theory considers the non-ergodic effect, since static and moving particles can be mixed within the OCT resolution volume in a highly heterogeneous sample. Finally, as both diffusion and flow can appear in a sample though spatially separated (e.g., the brain where both translational blood flow and diffusive intracellular organelle motions appear), the theory should be able to distinguish diffusive motions from translational flow and measure either the diffusion coefficient or the velocity according to the distinction. Although each of these concerns on DLS have been studied separately in the literature: the effect of the finite sample volume [7–9], the non-ergodic effect [10,11], and the mixture of diffusion and flow [12]; no theory has considered all of these simultaneously and has related these to OCT. Therefore, we derived the field autocorrelation function directly from the phase-resolved OCT signal under the assumption that static and moving particles are mixed in the resolution volume of OCT and the moving particles can exhibit either diffusive or translational motion.

This paper describes the derivation of the DLS-OCT theory, from the OCT signal when a single particle moves, to the field autocorrelation function when many particles exhibit heterogeneous dynamics. We propose a fitting algorithm to accurately estimate the diffusion coefficient and velocity from the autocorrelation function. The DLS-OCT theory has been validated through numerical simulations starting with the numerical position data of particles. Further, measurements of the axial and transverse velocities and the diffusion coefficient have been validated with phantom experiments. Based on these validations, we present an example application where DLS-OCT imaging was performed on the living rodent brain cortex, resulting in a velocity map showing blood flow and a diffusion map showing non-translational motions at the boundary of vessels.

## 2. Materials and methods

#### 2.1 Spectral-domain optical coherence tomography system

We used a spectral-domain optical coherence tomography (SD-OCT) system (Thorlabs, Inc) and optimized it for DLS-OCT imaging. The details of the OCT system are described in our previous publication [13]. SD-OCT generally has several advantages over time-domain OCT, one of which is to provide phase information that is very sensitive to the motion of particles [14,15]. The system employs a large-bandwidth near-infrared light source (1310 ± 170 nm in wavelength) for a large imaging depth and high spatial resolution (an axial resolution of 3.5 μm in biological tissue). The transverse resolution depends on which objective lens is used. Our 10X objective lens was used for the phantom and animal experiments described in this paper, leading to the transverse resolution of 3.5 μm in biological tissue. The scanning speed is 47,000 A-scans/s.

#### 2.2 Animal preparation

A detailed description of the animal preparation is given in our previous publication [13]. All experimental procedures were approved by the Massachusetts General Hospital Subcommitee on Research Animal Care. Sprague Dawley rats (250-300 g) were anesthetized with isofluorane and ventilated with a mixture of air and oxygen during surgical procedures. A craniotomy was performed, and a 3 mm × 3 mm area over the somatosensory cortex was exposed. The dura was carefully removed, and then the brain surface was covered with agarose gel and a glass cover slip. Physiological signs such as heart rate, body temperature and blood pressure were continuously monitored during surgery and during the experiment.

## 3. Results

#### 3.1 DLS-OCT theory

Our derivation starts from the SD-OCT signal when only one particle exists in the sample. SD-OCT generally measures the spectrum of the interference between the reference beam and light reflected from the sample. With reasonable approximations and the assumption of Gaussian-shaped source light spectrum, the Fourier transform of the interference spectrum results in the signal arising from a voxel (i.e., the complex-valued field reflectivity) where a single particle moves as *z*_{1}(*t*) along the z axis with respect to the center of the voxel:

*R*

_{01}is proportional to the field reflectivity of the particle. The other quantities are defined in Table 1 . The spectral width of the light source Δ

*q*determines the axial length of the coherence gating of SD-OCT (i.e. the Gaussian function in the magnitude of Eq. (1)). The transverse resolution is generally determined by the objective lens used in the focusing optics. Therefore, the signal from a particle moving in the three-dimensional space is

Now, we derive the field autocorrelation function of the OCT signal (Eq. (2)) when a single particle moves in the voxel. The movement of a particle can be defined by the self-part of the Van Hove space-time correlation function [16], that is, the probability that the particle which was initially at position **r _{i}** at time

*t*is found at position

_{i}**r**at time

_{f}*t*, ${P}_{1}({r}_{f},{t}_{f}|{r}_{i},{t}_{i})$. The probability of a particle exhibiting both translational and diffusive motions with the velocity

_{f}*v*and the diffusion coefficient

*D*can be formulated by

*E*[ ] means the average over random initial positions of the particle and ${v}_{t}^{2}={v}_{x}^{2}+{v}_{y}^{2}$ is the squared transverse speed. Here, we assumed that the initial position of the particle is evenly distributed through the whole space. In the regime of this study, where the maximum time lag

*τ*

_{max}= 25 Δ

*t*limits the maximum measureable velocity to 1/

*h*/

*τ*

_{max}~5 mm/s and the corresponding maximum measurable diffusion coefficient to (

*h*/

*q*)

^{2}

*v*

_{max}

^{2}

*τ*

_{max}~10 μm

^{2}/s, the factors that couple the diffusive movement and the small voxel size, $1/(1+4{h}_{t}^{2}D\tau )$ and $1/(1+4{h}^{2}D\tau )$, are very close to 1 (the minimum was 0.9969). Therefore, Eq. (4) becomes

*N*is the number of particles. The field autocorrelation function of this signal is the same as Eq. (5),

When particles exhibit heterogeneous dynamics in the voxel, those dynamics can be categorized into three groups (static (S), diffusing or flowing (F), and entering or exiting (E); see Fig. 1(a)
). We now assume an isotropic voxel (*h _{t}* =

*h*) for simplicity. The OCT signal from the voxel is

*N*

_{Ω}is the number of Ω-type particles in the voxel (Ω = S, F, or E). Here, the

*r*(

_{Ej}*t*) and

*z*(

_{Ej}*t*) will randomly vary by the definition of the E-particles. Again, we assumed ${\left|{R}_{S0j}\right|}^{2}\cong {\left|{R}_{F0j}\right|}^{2}\cong {\left|{R}_{E0j}\right|}^{2}=1$ for simplicity, leading to

*w*(

*t*) is complex-valued white noise. In practice, the third term may include stochastic fluctuations in the OCT signal due to imperfections of the OCT system.

The field autocorrelation function of this OCT signal has several groups of coupled terms. Note that static and moving particles are mixed in this case (i.e., non-ergodic) and thus the system is not stationary; that is, $\u3008{R}^{*}(t)R(t+\tau )\u3009\ne {R}^{*}(0)R(\tau )$.

*δ*(

*τ*) is the delta function,

*r*is the

_{jFk0}*j*-axis component of the initial position of the

*k*-th F-particle (

*j*=

*x*,

*y*,

*z*). In deriving Eq. (11), we assumed that the

*m*-th S-particle and

*k*-th F-particle are uncorrelated such that $P[({r}_{Fk},t+\tau )\cap ({r}_{Sm},t)]=P({r}_{Fk},t+\tau )P({r}_{Sm},t)$, leading to

*q*/Δ

*q*)

^{2}>> 1, which is satisfied in this study as (

*q*/Δ

*q*)

^{2}> 300. Consequently, the field autocorrelation function of mixed dynamics is given by

*M*(

_{S}**r**),

*M*(

_{F}**r**),

*v*(

_{t}**r**),

*v*(

_{z}**r**), and

*D*(

**r**) are the parameters of particle dynamics to be estimated for each position.

*M*is the composition ratio of static particles,

_{S}*M*is that of flowing/diffusing particles,

_{F}*v*is the transverse component of the flow velocity,

_{t}*v*is the axial component, and

_{z}*D*is the diffusion coefficient.

The general behavior of this autocorrelation function is illustrated in Fig. 1(b). For a given voxel, the constant term (*M _{S}*) is the center of rotation, and

*M*is the initial amplitude of rotation. The speed of rotation is determined by the axial velocity-dependent phase term (${e}^{iq{v}_{z}\tau}$), which is similar to the one used in Doppler OCT. The decay of the amplitude of rotation is governed by the diffusion-oriented decay term (${e}^{-{q}^{2}D\tau}$) that appears in the traditional DLS analysis for diffusion measurements. Interestingly, the decay is also affected by the velocity-dependent term (${e}^{-{h}_{t}^{2}{v}_{t}^{2}{\tau}^{2}-{h}^{2}{v}_{z}^{2}{\tau}^{2}}$), which made the decay differ from a simple exponential decay. This term, which is similar to that introduced in previous DLS studies [7–9,12], enabled us to estimate the transverse component of the flow velocity in addition to the axial component as provided by the phase term. We assumed that the composition ratios of the static and flowing/diffusing particles (i.e.,

_{F}*M*and

_{S}*M*) do not vary during a short correlation time (e.g., 0.5 ms in this study).

_{F}#### 3.2 Fitting algorithm

We developed and optimized an algorithm for estimating the dynamic parameters through fitting the theory to the autocorrelation data. The fitting has to determine five independent coefficients (*M _{S}*,

*M*,

_{F}*v*,

_{t}*v*, and

_{z}*D*) such that they minimize the sum of squared residuals (i.e., maximizing the coefficient of determination,

*R*

^{2}):

As shown in Fig. 2
, we found it inaccurate to use a simple fitting algorithm containing a single simplex search method [17] to find a global minimum in the five-dimensional parameter space. For this reason, we developed an advanced fitting algorithm. First, we reduced the number of parameters by determining *D*, *v _{z}* and

*v*directly in the least square manner given

_{t}*M*and

_{S}*M*:

_{F}*M*-

_{S}*M*)

_{F}*δ*(

*τ*). The radius of rotation |

*g*(

*τ*)-

*M*| was weighted in calculating

_{S}*v*to reduce the contribution of noisy movements of small-radius data points. Also,

_{z}*g*

_{1}(

*τ*) was weighted in calculating

*v*and

_{t}*D*to minimize distortion of the residuals that is caused by applying the logarithm.

*D*was forced to have a positive value.

Even when fitting for the two remaining parameters (*M _{S}* and

*M*), searching for the minimum was very sensitive to the initial guess (i.e., there was more than one local minimum). Therefore, our algorithm guesses and tests three different initial guess of

_{F}*M*and

_{S}*M*.

_{F}For the initial guesses of *M _{S}* and

*M*, an initial estimation of 1-

_{F}*M*-

_{S}*M*=

_{F}*M*was helpful.

_{E}*M*was simply estimated from the first three non-zero timelag data under the assumption that the first four data points can be approximated to a second-order polynomial. That is,

_{E}*τ*=

_{k}*k*Δ

*τ*is the k-th timelag.

*M*was forced to have a value between zero and one.

_{E}For the first initial guess of *M _{S}*, the center of rotation was used. We derived an equation to find the center of rotation (

*A*+

*iB*) in the complex plane for a given set of autocorrelation data points

*g*(

*τ*) =

_{k}*a*+

_{k}*ib*such that the center minimizes the variation in the radius of rotation (i.e., the distance from the center). This equation (Eq. (18)) was used to choose an initial guess of

_{k}*M*as the center of rotation of the latter half of the autocorrelation function. We used only the latter half because the radius of rotation varied largely during the first half of the timelag when flow and/or diffusion were large. The real part of the center of rotation was chosen as the initial guess for

_{S}*M*.

_{S}*k*>

*N*/2 and

_{τ}*N*is the number of autocorrelation data points (i.e., the number of timelag points).

_{τ}*M*was determined by

_{F}*M*= 1-

_{F}*M*-

_{S}*M*. Both

_{E}*M*and

_{S}*M*were forced to have a value between zero and one.

_{F}For the second initial guess of *M _{S}* and

*M*, a mesh of

_{F}*M*and

_{S}*M*was tested. For each pair of

_{F}*M*and

_{S}*M*in the mesh grid,

_{F}*D*,

*v*and

_{z}*v*were determined (Eq. (16)) and then

_{t}*R*

^{2}was calculated (Eq. (15)). One pair of

*M*and

_{S}*M*was chosen as the second initial guess where it maximizes

_{F}*R*

^{2}. As the third initial guess, we chose a set of typical values that empirically turned out to result in good fitting,

*M*=

_{S}*M*= (1-

_{F}*M*)/2.

_{E}For each initial guess, the simplex search method was used to find *M _{S}* and

*M*minimizing the sum of squared residuals, where

_{F}*D*,

*v*and

_{z}*v*were directly determined from

_{t}*M*and

_{S}*M*during each iteration (Eq. (16)). This process resulted in three sets of

_{F}*M*,

_{S}*M*,

_{F}*D*,

*v*and

_{z}*v*. We chose the one with the maximum

_{t}*R*

^{2}. A diagram summarizing our fitting algorithm is presented in Fig. 3 .

We tested if the algorithm results in better estimation of the parameters than the simple algorithm shown in Fig. 2. Autocorrelation data were generated for various values of *M _{S}*,

*M*,

_{F}*D*,

*v*= (

*v*

_{z}^{2}+

*v*

_{t}^{2})

^{1/2}and

*v*/

_{z}*v*(the cosine of the flow angle) according to Eq. (13). The estimated parameters agreed with the given values much better than with the simple algorithm (Fig. 4 ).

#### 3.3 Validation through numerical simulation

Our DLS-OCT theory has been validated using numerical simulations. As summarized in Fig. 5
, we numerically generated two-dimensional (transverse and axial axes) position data of 100 particles for 100 time steps (Δ*t* = 20 μs) for various parameter values (Table 2
). Position data were used to generate SD-OCT signals, and autocorrelation data were obtained from the signals. The autocorrelation data were fit using the algorithm described in the above section, leading to the estimation of the dynamic parameters.

First, we validated the estimation of the flow velocity for various parameters while fixing the diffusion coefficient to zero. As can be seen in Fig. 6(a)
, the autocorrelation function behaved as our theory predicted across several combinations of the parameters. Note that the center of rotation moved when the number of static particles changed. In this study, we mainly investigated the diffusion coefficient and velocities (*D*, *v*, and *v _{z}*) multiplied by

*M*, because a small

_{F}*M*in practice can make the estimation result in unreasonably large diffusion coefficient or velocities. Since

_{F}*M*is dimensionless,

_{F}*M*and

_{F}D*M*have the units of μm

_{F}v^{2}/s and mm/s, respectively. As a result, the absolute velocity and the axial velocity were estimated close to the true values (Fig. 6(b)). In addition, the diffusion coefficient was estimated negligible (0.20 ± 0.27 μm

^{2}/s). These results verify that the field autocorrelation function of the OCT signal numerically obtained from the position data of flowing particles exhibits the behavior that our theory predicts and thus the flow velocity can be accurately estimated by fitting the model of Eq. (14) to data.

Next, we validated the estimation of the diffusion coefficient. The position data of diffusing particles was generated with the fractional Brownian motion synthesis algorithm [18]. Our model was fit to the numerical autocorrelation function data (Fig. 7(a)
), and the *M _{F}*-weighted diffusion coefficient was estimated close to the true values (Fig. 7(b)). In this simulation, the absolute velocity was estimated negligible (0.24 ± 0.47 mm/s).

Results from these numerical simulations confirm that the DLS-OCT theory enables us to effectively distinguish diffusion from flow for a given autocorrelation function and accurately estimate the diffusion coefficient or the absolute and axial velocities according to the distinction. The distinction is mainly attributed to the difference in the waveforms of decay of the *M _{F}*-term (Fig. 7(c)) as predicted in Eq. (14).

#### 3.4 Validation through phantom experiment

DLS-OCT measurement of the velocity and diffusion coefficient was further validated through phantom experiments. We used the SD-OCT system described in the Materials and methods. In order to verify measurement of the flow velocity, a piezoelectrically actuated static sample was used to simulate axial movements of particles while transverse movements were implemented by galvanometric lateral scanning of OCT beam. For each combination of the axial and transverse velocities, we acquired the OCT signal and obtained 3D autocorrelation function data, *g*(*x*,*z*,*τ*). The autocorrelation function was averaged over neighboring 3 × 3 voxels and then was used to estimate the dynamic parameters, leading to 2D maps of *M _{S}*(

*x*,

*z*),

*M*(

_{F}*x*,

*z*),

*D*(

*x*,

*z*),

*v*(

_{z}*x*,

*z*),

*v*(

_{t}*x*,

*z*), and

*R*

^{2}(

*x*,

*z*), where

*R*

^{2}is the coefficient of determination (Eq. (15)). The values of a bad-fitting voxel (

*R*

^{2}<0.5) were replaced with the mean value of the neighboring good-fitting voxels. As noise can result in a small

*M*and large

_{F}*v*or

*D*, the velocity and diffusion maps were weighted by

*M*while

_{F}*M*smaller than 0.1 was forced to zero. Then, each map of the dynamic parameters was convolved with a 2D Gaussian kernel (10 μm in diameter). As a result, the absolute and axial velocities were reliably measured across various true velocities and flow angles (Fig. 8(a) ). Microsphere samples of 0.1 and 1 μm in diameter were used for validation of the diffusion measurement, where the measured diffusion coefficient agreed with the theoretical values given by the Einstein-Stokes equation (Fig. 8(b)). For this phantom experiment, monodisperse polystyrene microspheres in 2.5% solids (w/v) aqueous suspension (Polysciences, Inc.) were used.

_{F}#### 3.5 DLS-OCT imaging of the rodent brain

Based on the validated DLS-OCT measurement of the velocity and diffusion, we performed DLS-OCT imaging of the living rodent cortex. A-scans were repeated 100 times at a fixed position, and the position was moved in a raster manner to scan the volume of 600 μm × 600 μm × 300 μm (*x* × *y* × *z*) of the cortical surface. This FOV consisted of 400 × 400 positions, leading to the scanning step size of 1.5 μm and a total scanning time of ~10 min. We chose 100 A-scans per position because our fitting algorithm worked well when the autocorrelation function data had ≥25 time points, and the measurement time of ~2 ms corresponding to 100 A-scans was sufficiently shorter than the characteristic time constants of the primary sources of motion artifacts (i.e., cardiac and respiratory motions) as revealed in our previous study [13]. The 4D (space and time) complex-valued field reflectivity of the sample was obtained, and then was used to produce the 4D (space and timelag) autocorrelation function data. The autocorrelation function was averaged over 3 × 3 × 3 neighboring voxels, since the theory considers the autocorrelation function being the average over not only the measurement time but also various initial positions of particles. The spatial averaging of the autocorrelation function is based on the assumption that the initial position of the particles will vary across the neighboring voxels whereas their dynamics will be similar. For this assumption to be as valid as possible, we chose the scanning step size slightly smaller than the conventional Nyquist sampling frequency. Analysis of the 4D autocorrelation function data led to 3D maps of the transverse and axial velocities, the diffusion coefficient, and the coefficient of determination (*R*^{2}). This processing took ~3 hours with our 200-core parallel processing cluster system.

As shown in Fig. 9(a) , the absolute velocity map clearly revealed the structure of the vascular anatomy and cerebral blood flow. The axial velocity map showed the axial component of the flow velocity, which looked very similar to conventional Doppler OCT images. The flow direction determined by the axial and transverse velocities agreed with the structural direction of vessels. This agreement between the structural and flow direction can be more clearly seen in the cross-sectional maps (Fig. 9(a), right).

When the diffusion map is overlaid with the absolute velocity map (Fig. 9(b)), we observed high-diffusion and low-velocity voxels at the boundary of vessels. Interestingly, the vessel boundaries also exhibited a characteristic low coefficient of determination (Fig. 9(c)). Single-plane images clearly showed the result that the vessel boundaries exhibit high diffusion, low velocity, and low *R*^{2}. At the circular cross-sections of the vessels (Fig. 9(c), blue arrows), high-velocity and high-*R*^{2} blood flow are surrounded by the high-diffusion, low-velocity, and low-*R*^{2} dynamics of the vessel boundaries. In particular, the low *R*^{2} means that the motion was neither translational nor diffusive; we hypothesize that it might be oscillatory due to the interaction between blood flow and the tension of vessel walls. As can be seen in the examples of the autocorrelation functions (Fig. 9(d)), the composition ratio of static particles (*M _{S}*; the center of rotation in the complex plane) increased as the voxel is located far from the center of the vessel. When the voxel is located at the vessel boundary, it exhibited the characteristic dynamics with a low

*R*

^{2}(i.e., bad fitting). As the present paper focuses on describing this new technique, a detailed interpretation of the brain imaging data will be described elsewhere.

## 4. Discussion

DLS-OCT has several advantages over conventional techniques. As for the application to blood flow imaging, laser Doppler flowmetry is used to measure blood flow at a fixed point [19,20], and its imaging corollary provides the 2D map of blood flow [21]. Doppler OCT enables 3D imaging of the axial flow velocity with microscopic resolution [22]. Compared to these techniques, DLS-OCT can simultaneously and independently measure the axial and transverse components of the flow velocity. Further, DLS-OCT can distinguish diffusion from flow, leading to a separate 3D map of the diffusion coefficient. In addition, DLS-OCT may offer additional contrast mechanisms (e.g., those based on *R*^{2} and *M _{F}*). The

*R*

^{2}map quantitatively images the degree of how much the motion is close to translational or diffusive ones; and for example it clearly revealed characteristic dynamics of the vessel boundaries. The

*M*map will quantify the fraction of moving particles in each voxel, which is similar to the mobile fraction suggested in non-ergodic DLS studies [10,11].

_{F}Each term in the DLS-OCT model we derived is similar to those predicted in several DLS studies. The velocity-dependent decay term is similar to that predicted in the study of the effect of the finite sample volume on DLS analysis [7–9], although the autocorrelation function in our theory was directly derived from the OCT signal whereas the finite size of the sample volume in the literature resulted from the illuminating and collecting optics. The linear decomposition of the static component and the flowing/diffusing component is similar to the decomposition shown in the study of DLS in non-ergodic media [10,11]. The combined diffusion-oriented and flow-oriented decay terms were similarly introduced in the study of DLS where diffusion and flow are mixed [12]. Therefore, our theory can be understood as a mathematical combination of the phase-resolved OCT signal, the finite sample volume DLS model, the non-ergodic DLS model, and the model for the mixture of diffusion and flow.

Our DLS-OCT theory assumed that the composition ratios of the static and diffusing/flowing particles (*M _{S}* and

*M*) do not vary during the measurement time. The validity of this assumption will depend on the measurement time and the magnitude of dynamics within a sample. For example, the assumption was reasonable when imaging blood flow of the brain with the measurement time of ~2 ms, because the blood flow velocity is typically 1-5 mm/s and thus leads to 2-10 μm displacement during the measurement time, smaller than or approximately comparable to the resolution volume. When we used a longer measurement time (200 ms), the autocorrelation function deviated from our model. In contrast, the fitting result was not reliable if only very short correlation times were used. Therefore, the correlation and measurement times for DLS-OCT imaging should be chosen carefully, taking into account both the scale of target dynamics and the fitting performance.

_{F}This study used the DLS-OCT theory for distinguishing whether the motion is translational or diffusive, not for measuring the mixture of the two motions. Nevertheless, the model seems to be able to measure the mixture of translational and diffusive motions. It was reported, however, that the diffusion is estimated inaccurately when mixed with large flow [12]. We also observed this tendency in a phantom experiment, where the diffusion coefficient was overestimated by ≈1.5 times when it was mixed with flow whose velocity is larger than >1 mm/s. In this study, the relative magnitudes of flow and diffusion can be estimated based on how much they contribute to the decay in the autocorrelation function. With the parameters used in this study, *v* = 0.1 and 1 mm/s approximately correspond to *D* = 0.004 and 0.4 μm^{2}/s, respectively. The coupling between diffusion and large flow will not be a critical problem in the studies investigating a sample where diffusion and flow are spatially separated. Meanwhile, considerable caution would be required if one wants to apply the present technique to the measurement of diffusion that is mixed with large flow within the resolution volume. On the other hand, the voxels with large blood flow in vessels often exhibited high diffusion, which can be either a real diffusive motion or decorrelation of the OCT signal that can be quantified by the exponential decay with the diffusion coefficient. This high diffusion in vessels might be attributed to the mixture of non-translational motion of red blood cells including rotation, turbulence in blood flow, and the variance in the velocity distribution within the resolution volume. Since this interpretation has not been yet validated, and since the present study used the DLS-OCT theory for determining whether the motion is translational or not, we overlaid the diffusion map with the velocity map as in Figs. 9(B) and 9(C). This overlay means that our analysis gave priority to flow over diffusion so that diffusion is of interest only at the voxels with low flow. Therefore, the coupling between diffusion and large flow was not an important concern at those vessel boundaries.

OCT uses coherence gating to collect light only scattered from the resolution volume and is known to effectively exclude multiply scattered photons [23]. However, strong multiple scattering can give rise to a distortion in the OCT signal, which often causes an undesired shadow of large vessels [13]. This multiple scattering also affected DLS-OCT estimation of dynamics at the voxels located beneath the surface vessels. This effect can be seen clearly by comparing the merged images in Fig. 9, where Fig. 9(b) presents the maximum projection over the whole depth while Fig. 9(c) only presents a single plane at the depth of 120 μm. High velocity and diffusion appeared in the shadow of the large surface vessels as shown in the single-plane image. This multiple scattering effect also can be found in the cross-section image of the velocity map of Fig. 9(a), where the absolute velocity was estimated large in the vessel shadow but the axial velocity was not. Although this multiple scattering would not cause a serious problem as one may generally build an *en face* image by including surface vessels as in Fig. 9(b), in the future it will be interesting to derive a DLS-OCT model that considers the effect of multiple scattering.

Future theoretical efforts should consider the possibility of measuring a mixture of translational and diffusive motions, and the effect of multiple scattering. Although there can be various modified versions of the DLS-OCT theory, the one described here will work well for 3D imaging of dynamics in a highly heterogeneous sample where static and moving particles can be mixed within the micrometer-scale resolution volume and the moving particles can exhibit either diffusion or flow, as long as it is used in the single-scattering regime.

## 5. Conclusion

We demonstrated a technique based on the integration of DLS and OCT for high-resolution 3D imaging of heterogeneous diffusion and flow. A theory for this purpose was proposed and validated with numerical simulations and phantom experiments. The DLS-OCT theory enabled us to simultaneously measure the axial and transverse velocities and the diffusion coefficient with the micrometer-scale resolution. We tested the utility of this technique by applying it to 3D *in vivo* imaging of translational blood flow and non-translational motion of the vessel boundary in the living animal brain cortex.

## Acknowledgments

This study was supported by the NIH (R01-NS057476, R01-EB000790, P01NS055104) and the AFOSR (MFEL FA9550-07-1-0101).

## References and links

**1. **N. A. Clark, J. H. Lunacek, and G. B. Benedek, “A study of Brownian motion using light scattering,” Am. J. Phys. **38**(5), 575–585 (1970). [CrossRef]

**2. **D. J. Durian, D. A. Weitz, and D. J. Pine, “Multiple light-scattering probes of foam structure and dynamics,” Science **252**(5006), 686–688 (1991). [CrossRef] [PubMed]

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

**4. **D. Huang, E. A. Swanson, C. P. Lin, J. S. Schuman, W. G. Stinson, W. Chang, M. R. Hee, T. Flotte, K. Gregory, C. A. Puliafito, and J. Fujimoto, “Optical coherence tomography,” Science **254**(5035), 1178–1181 (1991). [CrossRef] [PubMed]

**5. **G. J. Tearney, M. E. Brezinski, B. E. Bouma, S. A. Boppart, C. Pitris, J. F. Southern, and J. G. Fujimoto, “In vivo endoscopic optical biopsy with optical coherence tomography,” Science **276**(5321), 2037–2039 (1997). [CrossRef] [PubMed]

**6. **Y. Chen, L. N. Vuong, J. Liu, J. Ho, V. J. Srinivasan, I. Gorczynska, A. J. Witkin, J. S. Duker, J. Schuman, and J. G. Fujimoto, “Three-dimensional ultrahigh resolution optical coherence tomography imaging of age-related macular degeneration,” Opt. Express **17**(5), 4046–4060 (2009). [CrossRef] [PubMed]

**7. **R. V. Edwards, J. C. Angus, and M. J. French, “Spectral analysis of the signal from the laser Doppler flowmeter: time-independent systems,” J. Appl. Phys. **42**(2), 837–850 (1971). [CrossRef]

**8. **D. P. Chowdhury, C. M. Sorensen, T. W. Taylor, J. F. Merklin, and T. W. Lester, “Application of photon correlation spectroscopy to flowing Brownian motion systems,” Appl. Opt. **23**(22), 4149–4154 (1984). [CrossRef] [PubMed]

**9. **T. W. Taylor and C. M. Sorensen, “Gaussian beam effects on the photon correlation spectrum from a flowing Brownian motion system,” Appl. Opt. **25**(14), 2421–2426 (1986). [CrossRef] [PubMed]

**10. **P. N. Pusey and W. Van Megen, “Dynamic light scattering by non-ergodic media,” Physica A **157**(2), 705–741 (1989). [CrossRef]

**11. **J. G. H. Joosten, E. T. F. Geladé, and P. N. Pusey, “Dynamic light scattering by nonergodic media: Brownian particles trapped in polyacrylamide gels,” Phys. Rev. A **42**(4), 2161–2175 (1990). [CrossRef] [PubMed]

**12. **A. B. Leung, K. I. Suh, and R. R. Ansari, “Particle-size and velocity measurements in flowing conditions using dynamic light scattering,” Appl. Opt. **45**(10), 2186–2190 (2006). [CrossRef] [PubMed]

**13. **J. Lee, V. Srinivasan, H. Radhakrishnan, and D. A. Boas, “Motion correction for phase-resolved dynamic optical coherence tomography imaging of rodent cerebral cortex,” Opt. Express **19**(22), 21258–21270 (2011). [CrossRef] [PubMed]

**14. **J. F. de Boer, B. Cense, B. H. Park, M. C. Pierce, G. J. Tearney, and B. E. Bouma, “Improved signal-to-noise ratio in spectral-domain compared with time-domain optical coherence tomography,” Opt. Lett. **28**(21), 2067–2069 (2003). [CrossRef] [PubMed]

**15. **R. Leitgeb, C. Hitzenberger, and A. Fercher, “Performance of fourier domain vs. time domain optical coherence tomography,” Opt. Express **11**(8), 889–894 (2003). [CrossRef] [PubMed]

**16. **L. Van Hove, “Correlations in space and time and Born approximation scattering in systems of interacting particles,” Phys. Rev. **95**(1), 249–262 (1954). [CrossRef]

**17. **J. C. Lagarias, J. A. Reeds, M. H. Wright, and P. E. Wright, “Convergence properties of the Nelder-Mead simplex method in low dimensions,” Siam J Optimiz **9**(1), 112–147 (1998). [CrossRef]

**18. **P. Abry and F. Sellan, “The wavelet-based synthesis for fractional Brownian motion proposed by F. Sellan and Y. Meyer: remarks and fast implementation,” Appl. Comput. Harmon. Anal. **3**(4), 377–383 (1996). [CrossRef]

**19. **M. D. Stern, “In vivo evaluation of microcirculation by coherent light scattering,” Nature **254**(5495), 56–58 (1975). [CrossRef] [PubMed]

**20. **U. Dirnagl, B. Kaplan, M. Jacewicz, and W. Pulsinelli, “Continuous measurement of cerebral cortical blood flow by laser-Doppler flowmetry in a rat stroke model,” J. Cereb. Blood Flow Metab. **9**(5), 589–596 (1989). [CrossRef] [PubMed]

**21. **B. M. Ances, J. H. Greenberg, and J. A. Detre, “Laser Doppler imaging of activation-flow coupling in the rat somatosensory cortex,” Neuroimage **10**(6), 716–723 (1999). [CrossRef] [PubMed]

**22. **J. A. Izatt, M. D. Kulkarni, S. Yazdanfar, J. K. Barton, and A. J. Welch, “In vivo bidirectional color Doppler flow imaging of picoliter blood volumes using optical coherence tomography,” Opt. Lett. **22**(18), 1439–1441 (1997). [CrossRef] [PubMed]

**23. **B. Karamata, M. Laubscher, M. Leutenegger, S. Bourquin, T. Lasser, and P. Lambelet, “Multiple scattering in optical coherence tomography. I. Investigation and modeling,” J. Opt. Soc. Am. A **22**(7), 1369–1379 (2005). [CrossRef] [PubMed]