Search tips
Search criteria 


Logo of optexpOSAOptics InfoBaseThis ArticleThis JournalSearchAuthor Information
Optics Express
Opt Express. 2010 February 1; 18(3): 2477–2494.
Published online 2010 January 22. doi:  10.1364/OE.18.002477
PMCID: PMC2837842

Quantitative cerebral blood flow with Optical Coherence Tomography


Absolute measurements of cerebral blood flow (CBF) are an important endpoint in studies of cerebral pathophysiology. Currently no accepted method exists for in vivo longitudinal monitoring of CBF with high resolution in rats and mice. Using three-dimensional Doppler Optical Coherence Tomography and cranial window preparations, we present methods and algorithms for regional CBF measurements in the rat cortex. Towards this end, we develop and validate a quantitative statistical model to describe the effect of static tissue on velocity sensitivity. This model is used to design scanning protocols and algorithms for sensitive 3D flow measurements and angiography of the cortex. We also introduce a method of absolute flow calculation that does not require explicit knowledge of vessel angles. We show that OCT estimates of absolute CBF values in rats agree with prior measures by autoradiography, suggesting that Doppler OCT can perform absolute flow measurements in animal models.

OCIS codes: (110.4500) Optical coherence tomography, (170.3880) Medical and biological imaging, (170.5380) Physiology, (170.1470) Blood or tissue constituent monitoring

1. Introduction

In the normal brain cerebral blood flow (CBF) is regulated to satisfy metabolic demand determined by brain activity. CBF is required for the supply of oxygen and nutrients and the removal of metabolites and heat. Alterations in CBF may contribute to the pathogenesis of stroke and Alzheimer's disease [1]. Therefore quantitative and localized measurements of CBF may significantly advance understanding of major cerebrovascular diseases. Rat and mouse models are of major importance in preclinical cerebrovascular research. However, a method for in vivo determination of absolute CBF values with high resolution in small animals is currently lacking. The gold standard for determination of regional blood flow is autoradiography. In this technique, a tracer, such as iodo[14C]antipyrine, is typically administered over a short period of time, followed by cardiac arrest and quick freezing. Autoradiography of frozen sections provides a representation of radioactivity levels in the different tissue layers, which can be converted to blood flow after incorporation of the time course of arterial blood radioactivity into an operational equation. Although autoradiographic methods provide three-dimensional spatial information, they contain no information about the temporal evolution of CBF changes [2]. Therefore, studies of disease progression or response to treatment using intra-animal comparisons cannot be performed, limiting the utility of this method.

Methods of flow measurement based on magnetic resonance imaging [3] and positron emission tomography [4] provide spatial maps of CBF but are limited in their temporal and spatial resolution. Two-photon microscopy can calculate both transverse velocity and linear density of red blood cells (RBCs) [5], which may be combined to calculate RBC flux. However, the determination of CBF requires knowledge of the flow profile as well as vessel orientation, which can be cumbersome to obtain. Scanning laser–Doppler can be used to obtain spatially-resolved relative CBF images by moving a beam across the field of interest [6]. Laser speckle imaging [7] has recently been applied for high spatiotemporal resolution relative measurements of CBF [8]. Laser speckle imaging uses speckle correlation times, which are inversely proportional to mean velocities, to obtain relative measurements of flow. However, it is difficult to relate correlation times (laser speckle) or the width of the power spectral density (laser Doppler) to absolute velocities and flow of red blood cells because the number of moving particles that the light interacts with and their orientations are unknown [9]. Furthermore, neither laser Doppler nor laser speckle offers depth specificity, since the wavelength, tissue optical properties, and the optical geometry determine the penetration depth. Therefore, a simple method that provides real-time three-dimensional, spatially-resolved, and absolute CBF measurements would aid in experimental studies of functional cerebral activation and cerebral pathophysiology.

Doppler Optical Coherence Tomography (OCT) overcomes many of the limitations of the aforementioned optical techniques by using the path delay of light to gate out multiply scattered light [1013]. Due to the “coherence gate” [14], OCT enables measurement of velocity averaged over well-defined volumes, determined by the confocal gate in the transverse direction and by the coherence length in the axial direction. Thus, Doppler OCT avoids the limitations associated with laser Doppler and laser speckle imaging.

Recently, Wang et al. demonstrated a method of OCT termed optical angiography (OAG) to visualize the three-dimensional cerebral circulation of adult living mice through the intact cranium [15]. OAG extends on an earlier observation by Ren et al. [16] that high-pass filtering of the complex OCT signal removes signal contribution from static scatterers and improves visualization of blood vessels. While Ren et al. used a real-valued high-pass filter, OAG uses a complex filter (isolating either positive or negative temporal frequencies) that can visualize either positive or negative velocities. Tao et al. noted that the minimum detectable velocity in this method decreases linearly with increased spatial oversampling [17]. Recently, measurements of velocity projections in the mouse cerebral cortex were performed by combining phase-resolved measurements with OAG [18]. Phase variance methods were recently applied to visualize major vasculature in the mouse brain [19]. Depth-resolved Doppler OCT measurements of relative flow changes in response to functional activation were also recently performed in the rat [20].

The utility of these OCT methods could be extended by 1) performing quantitative flow measurements, and 2) achieving comprehensive visualization of the cortical vasculature down to the capillary level. In this study, we develop a quantitative statistical model for static tissue and validate the model in a scattering phantom. Using this model to guide our signal processing, we perform quantitative, absolute flow measurements in the rat cortex. We measure flow values comparable to those obtained in rats under similar anesthetic conditions using autoradiography.

2. Experimental methods

2.1 System description

A spectral/Fourier domain OCT microscope for in vivo imaging of the rat cerebral cortex was designed [Fig. 1(a) ]. A superluminescent diode (model EX8705-2411, Exalos Inc.) with a center wavelength of 856 nm, bandwidth of 54 nm and a power of 5 mW was used as the light source. The axial (depth) resolution was approximately 6.4 µm in air (4.8 µm in tissue) after spectral shaping. The power on the sample was 1.5 mW, enabling a maximum sensitivity of 99 dB. A home-built spectrometer using a 2048 pixel, 12-bit line scan camera (Aviiva SM2 camera, e2v semiconductors) enabled an imaging speed of 22,000 axial scans per second and an axial imaging range of 2.5 mm in air. Quantitative flow measurements (Fig. 5 and Fig. 7 ) were performed using a fast 3D OCT raster scan (protocol 1, 160 images with 840 axial scans per image). A slow 3D OCT scan (protocol 2, 450 images with 32760 axial scans per image) was performed to acquire the OCT angiograms (Fig. 6 ). The protocols are summarized in
Fig. 1(d). Unless otherwise stated, the transverse resolution was 12.0 µm and the scanned area was 800 µm. For Fig. 6(c) only, the transverse resolution was 3.6 µm and the scanned area was 240 µm. The transverse resolution is defined as the 1/e2 diameter of the intensity profile in the focal plane.

Fig. 1
Schematic of the OCT system. (a) The optical path incorporates two steering mirrors (M1 and M2) in the beam path to achieve precise alignment of the OCT beam. (b) Zoom of the focal plane of the objective lens showing a point scatterer emitting a spherical ...
Fig. 5
Procedure for calculating flow without explicit calculation of vessel angle. (a) Blood vessels can be oriented at any angle [var phi] relative to the incident light. (b) The area of the vessel in the en face (xy) transverse plane is inversely proportional ...
Fig. 7
Absolute values of flow measurements in the rat cortex in ascending venules. (a) Locations for flow measurements ( Media 1) and (b) absolute flow measurements at the designated locations. Standard errors over 10 repeated measurements are shown.
Fig. 6
(a) OCT maximum intensity projection (MIP) angiogram of the rat somatosensory cortex showing comprehensive visualization of the vasculature at 12.0 µm transverse resolution. (b) zoom of arterial anasthmosis, showing visualization of surface vessels ...

A standard microscope design was used, where the galvanometer mirrors (Cambridge Technology Inc., Model 6230M) were relay imaged onto the back focal plane of a water immersion microscope objective (Olympus 20X, 0.95 NA). In order to minimize Doppler shifts from scanning the beam, two mirrors (M1 and M2) were used for alignment. The mirrors were adjusted to center the beam on all optical elements, while minimizing the phase difference between adjacent axial scans when scanning across a scattering phantom in the x and y directions. This alignment procedure ensured the accuracy and precision of flow measurements by removing residual Doppler shifts associated with transverse scanning and reducing the phase noise. The quantitative justification of this alignment procedure is discussed below.

2.2 Animal preparation

Sprague–Dawley rats (250 g – 350 g) were prepared under isoflurane anesthesia (1.5–2.5% in 25% Oxygen, 75% air). After tracheotomy and femoral arterial and venous catheterization, the head was fixed in a stereotaxic frame, the scalp was retracted, and a ~5 mm×5 mm area of skull over the somatosensory cortex was thinned to translucency. The IVth ventricle was opened to relieve intra-cerebral pressure, and the thinned skull area was removed along with the dura. 1.5% agarose (Fluka) in artificial cerebrospinal fluid [21] was placed onto the exposed cortex and covered with a small glass coverslip. The window was then sealed relative to the skull with dental acrylic to minimize brain motion. Following surgery, isoflurane was discontinued, and replaced with intravenous delivery of alpha-chloralose (30–40 mg/kg/h). Throughout data collection, the animals were mechanically ventilated with air and oxygen, and systemic blood pressure was monitored and maintained within the normal range (100 mm Hg±20 mm Hg). Blood gas values were measured when necessary to ensure good systemic physiology (PaCO2 between 35 and 45 mmHg, PaO2 between 90 and 120 mmHg, and pH between 7.35 and 7.45). All animal procedures were reviewed and approved by the Subcommittee on Research Animal Care at Massachusetts General Hospital, where these experiments were performed.

2.3 Noise and error sources

The effect of shot noise on differential phase measurements has been described in detail by others [22, 23]. In shot noise-limited conditions, the minimum measurable phase shift in the presence of photon noise (a zero-mean, complex Gaussian random process) is given by SNR-0.5 [24]. In addition to shot noise, there are additional factors that cause phase noise when imaging scattering tissue in vivo with actively controlled galvanometers. These include physiologic motion (respiration, heart rate, etc.) and galvanometer jitter. The effects of these additional sources of noise can be divided into axial and transverse components. A shift in axial position causes both a phase shift and decorrelation. Phase shifts, caused by axial position shifts on the order of λ/4, may be corrected by motion correction algorithms [25] since every axial position experiences the same phase shift. Decorrelation in the axial direction (caused by axial position shifts on the order of the axial resolution), leads to a different phase shift at axial positions separated by more than the axial resolution. Similarly, a shift in transverse position causes decorrelation [26], also leading to a different phase shift at axial positions separated by more than the axial resolution. However, small amounts of positional noise predominantly cause axial phase shifts because an axial position shift of λ/4 leads to a phase shift of π, whereas position shifts must be on the order of the resolution (typically greater than several λ in both the axial and transverse directions) in order to cause significant decorrelation. It should be noted that algorithms exist to correct for axial phase shifts, while decorrelation caused by transverse motion is typically not corrected.

As shown in Fig. 1(c), galvanometer motion can affect phase measurements if the direction of the OCT scan direction does not lie in the plane of the wavefront of the incident beam. The OCT scan direction is defined as the direction of the three-dimensional displacement vector of the Gaussian beam during beam scanning, and is shown as a red arrow in Fig. 1(c). In general, the OCT scan direction forms a non-zero angle θ relative to the wavefront of the incident beam. The effects of this misalignment are twofold – 1) a projection of galvanometer angular “jitter” onto the axial direction, and 2) a Doppler shift as the OCT beam is scanned across static tissue, leading to a Doppler frequency shift from static tissue. The transverse displacement of the beam (vxT) is related to the incremental change in optical angle caused by the galvanometer (d[var phi]) by vxT =Ld[var phi] where L= (f1/f2)fOL. As shown in
Fig. 1(c), any displacement along the axis of the scan is related to a displacement along the z axis by dz =sin(θ)vxT. Any z displacement causes a phase shift given by dΦ=(4πn/λ)dz, where n is the refractive index and λ is the wavelength. Therefore, the standard deviation of galvanometer jitter (σ[var phi]) can be related to the standard deviation of phase fluctuations caused by the corresponding z displacement (σΦ): σΦ=(4πnL/λ)sin(θ)σ[var phi]. Therefore, by minimizing the angle θ through careful optical alignment it is possible to minimize axial phase shifts caused by galvanometer jitter. In our setup we accommodated two steering mirrors in the beam path (M1 and M2), as shown in Fig. 1(a). The phase shifts along the x and y axes were monitored in real time, and the mirrors were adjusted to minimize the scanning-induced phase shift as the beam was scanned across a scattering phantom along each axis sequentially. Because both galvanometer mirrors cannot be simultaneously imaged at the objective back focal plane using our microscope design, it was impossible to achieve θ=0° at all transverse positions. However, this procedure achieved acceptably low residual Doppler shift levels for quantitative flow measurements. For qualitative angiography, correction of phase shifts due to bulk axial motion was performed as described in Section 2.7.

2.4 Bandwidth of static tissue and high-pass filtering

In order to understand the frequency sensitivity required for capillary visualization, we characterize the bandwidths of static and dynamic components of scattering from stationary tissue. For this purpose, we view the spatial speckle pattern generated by tissue as a complex white Gaussian random process filtered by the transverse point spread function [27].

As both temporal and spatial relationships are important, we characterize the speckle pattern temporally and spatially by autocorrelation functions (Rx(x) and Rt(t)) and power spectral densities (Px(fx) and Pt(ft)). The variables fx and ft correspond to spatial and temporal frequency, respectively. The temporal and spatial autocorrelation functions are given by the Fourier transforms of the temporal and spatial power spectral densities.

equation m1

The spatial and temporal autocorrelation functions have magnitudes with even symmetry about the origin and phases with odd symmetry about the origin.

equation m2

Because the speckle pattern is generated by filtering of white Gaussian noise, its autocorrelation and power spectral density are determined solely by the power in the white Gaussian noise process and parameters of the transverse point spread function h(x). In order to calculate h(x), the transverse point spread function, it is necessary to calculate the electric field profile. Assuming Gaussian beam propagation, the incident electric field is given by:

equation m3

For convenience we have normalized Ei(x,y,z) so that the norm squared integrated over any xy plane is equal to 1. This ensures that the power in the Gaussian beam (proportional to the norm squared of the electric field) is independent of wi, the width of the beam waist.

equation m4

As shown in Fig. 1(b), a point scatterer located at a position (x=x0,y=0,z=z0) will be excited by the incident field and emit a spherical wave which is eventually coupled back to the single mode optical fiber. Therefore, the amplitude and phase of the emitted light is determined by the amplitude and phase of the incident electric field. In order to calculate the field amplitude and phase coupled back to the single mode fiber it is necessary to perform an overlap integral of the fiber mode and the backscattered field. While it is possible to perform this integration in any plane, it is most convenient to perform the integration in the plane of the point scatterer. The field of the point scatter in the plane z=z0 can be represented by

equation m5

K is an arbitrary complex constant. The complex conjugate is necessary because the direction of propagation is reversed for the backscattered field. Therefore,

equation m6

For convenience we assume that the point scatterer is located in the focal plane (z0 = 0) and define the point spread function as:

equation m7

Therefore, the point spread function is essentially given by the incident electric field squared.

Then, the autocorrelation function and power spectral density of the speckle pattern may be expressed as follows:

equation m8

In the above expression σ2 represents the variance (power) of white noise Gaussian random process and * denotes convolution. The power (variance) of the speckle pattern, Rx(0), varies inversely as the beam width cubed. The fact that the OCT beam is scanned (x=vxt) imposes a relationship between spatial and temporal autocorrelation functions and power spectral densities.

equation m9

The interpretation of this equation is that the temporal power spectral density of stationary tissue can be either stretched or compressed by increasing or decreasing the velocity, respectively. The variance of the temporal speckle pattern, or the integral of the power spectral density over all frequencies, is independent of vx.

In order to verify these relationships we experimentally measured the spatial and temporal autocorrelation functions on a scattering phantom using different transverse scanning velocities. An OCT image of an Ocean Optics WS-1 diffuse reflectance standard is shown in Fig. 2(a) . A comparison of the measured and theoretical autocorrelation functions is shown in Fig. 2(b). The theoretical autocorrelation function is obtained from the measured Gaussian beam width and the theory presented above. Figure 2(c) shows measured spatial autocorrelation functions at different scanning speeds, demonstrating good agreement.
Figure 2(d) shows temporal autocorrelation functions at different scanning speeds, demonstrating an approximately twofold decrease in width as the speed is doubled.
Figure 2(e), 2(f) shows the spatial and temporal power spectral densities at different scanning speeds. The power spectral densities were obtained by Gaussian fitting of the autocorrelation functions followed by Fourier transformation. All measurements are normalized to the peak value of the slowest velocity (vx = 0.55 mm/s) curve to allow easy comparison.

Fig. 2
Slower scanning speeds improve frequency sensitivity by decreasing the width of the power spectral density of static tissue. (a) Image of scattering phantom used for autocorrelation and power spectral density measurements. (b) Measured and theoretical ...

While Fig. 2 shows the autocorrelation and power spectral density for stationary tissue, similar results hold for tissue moving in the axial direction. For tissue generating a Doppler frequency of f0, assuming that the dominant source of decorrelation is transverse scanning, the magnitude of the autocorrelation remains unchanged while the phase is given by 2πf0t. This leads to a shift of the power spectral density such that the maximum value is located at f0. However, it is important to note that for moving red blood cells there will be additional sources of decorrelation leading to increased power spectral density width including axial speckle decorrelation caused by RBC motion in the axial direction, transverse speckle decorrelation caused by RBC motion in the transverse direction, and a random RBC velocity distribution.

In the rat somatosensory cortex, capillaries exhibit significant fluctuations in speed and flux, with speed varying between 0 and 1 mm / s and red blood cell flux varying between 0 and 100 s−1 [5]. Given the illumination and detection geometry, only velocities with components along the z axis generate Doppler shifts. Therefore, in order to map capillaries as well as large blood vessels, it is necessary to measure velocities much lower than 1 mm / s. For instance, a capillary velocity of 0.5 mm / s at an angle of 85 degrees to the probe beam generates a Doppler shift of approximately 150 Hz. In addition to Doppler frequency shifts caused by axial motion, broadening of the power spectral density is caused by speckle decorrelation. The effect of red blood cell motion (e.g. in capillaries) is to cause a speckle decorrelation, which leads to a broadening of the power spectral density. As a simple order-of-magnitude calculation of this effect, assuming a transverse motion of ~0.5 mm /s and referring to Fig. 2, the decorrelation time is on the order of 20 ms [Fig. 2(d)], and the bandwidth generated is on the order of 50 Hz [Fig. 2(f)]. Broadening due to decorrelation from motion in the axial direction is also possible. The effect of power spectral density broadening is distinct from Doppler shifts, since it results in an increase in the width of the power spectral density as opposed to a shift in the center frequency of the power spectral density. In general red blood cell flow is accompanied by both frequency shifts and broadening.

The above discussion established estimates for the bandwidths of dynamic and static tissue under different conditions. The requirements for flow imaging and angiography are different; consequently they require different scanning protocols. Firstly, for flow quantification it is not necessary to measure capillary flow, since virtually all flow in the cortex can be traced to a surface artery or vein. In addition, repeated flow imaging may be advantageous to reduce physiologic flow noise. Therefore, a transverse scanning speed of ~21 mm / s was chosen for flow measurements, leading to a 3D OCT protocol of 160 images with 840 axial scans per image (protocol 1). This corresponds to an oversampling factor (defined as the number of A scans per FWHM of the autocorrelation function) of 10.5. Due to the short scanning time, this protocol was repeated 10 times in succession, reversing the direction of fast axis scanning for alternate 3D OCT scans. For angiography, it is desirable to reduce the temporal bandwidth of stationary tissue significantly to visualize capillary flows. A transverse scanning speed of 0.55 mm / s was chosen for angiography, leading to a 3D OCT protocol of 450 images with 32760 axial scans per image (protocol 2). This corresponds to an oversampling factor of 410. At this speed the bandwidth of stationary tissue was less than 50 Hz, as shown in Fig. 2(f).

A digital filter was designed to remove stationary scattering components from the OCT image [16]. The digital filtering procedure assumes that the contribution from stationary and moving scatterers can be separated by a linear operation. (This is not the case if there is spatially correlated or multiple scattering, where the contribution of a stationary scatterer and a moving scatter in combination cannot be obtained by superposition of their individual contributions.) The ability of angiographic methods to visualize slow flow is determined by the bandwidths and magnitudes of stationary and dynamic scattering as well as the cutoff frequency of the high-pass filter. Because the amount of stationary scattering in the brain exceeded the amount of dynamic scattering, a relatively high cutoff frequency was chosen for the filters in Fig. 3 . A filter with a lower cutoff frequency will lead to a higher background from stationary tissue, but will enable imaging of slower flows, while a filter with a higher cutoff frequency will lead to a lower background from stationary tissue, but may also miss slow flows. An additional criterion in choosing the filter is that the filter should preserve the spatial resolution of the flow image. For protocol 1, the temporal extent of the filter was limited to ~0.2 milliseconds, corresponding to a spatial extent of ~4 µm. For protocol 2, the temporal extent of the filter was limited to ~3 milliseconds, corresponding to a spatial extent of ~1.7 µm and resulting is negligible loss of spatial resolution.

Fig. 3
. High-pass filtering was performed along the transverse direction to separate moving and stationary scatterers. (a) Real-valued filter kernel used for high-pass filtering and (b) Fourier transform of the filter kernel, or frequency response, for protocol ...

An additional consideration is the effect of high-pass filtering on power spectral density and frequency estimation. Figure 4 shows the theoretical temporal power spectral densities before and after digital high-pass filtering for a range of Doppler shifts from 0 to 6 kHz. The width of the temporal power spectral density was assumed to be the same as that of static tissue for this figure. In other words, the only source of power spectral density broadening was assumed to be transverse scanning (Note that this assumption is not satisfied at high Doppler frequencies, corresponding to high axial velocities, due to additional decorrelation caused by axial motion). As can be seen by comparing power spectral densities before and after high-pass filtering in Fig. 4(a) (protocol 1), high-pass filtering leads to two effects: a signal loss, and a distortion of the power spectral density. Because of this distortion, when the mean frequency is measured after high-pass filtering, a significant error results below 2 kHz [Fig. 4(b)]. By comparison, Fig. 4(c), 4(d) shows that for protocol 2, there is negligible distortion and error in the mean frequency caused by high-pass filtering. Therefore, high-pass filtering leads to a distortion of the power spectral density, which can be mitigated by slower scanning speeds.

Fig. 4
. Digital high-pass filtering results in a Doppler mean frequency bias. Comparison of power spectral densities for protocol 1 (a-b) and protocol 2 (c-d) before and after digital high-pass filtering. The input center frequencies (0.0, 1.5, 3.0, 4.5, 6.0 ...

The preceding analysis does not account for continuous movement of the sample beam during the camera integration time. However, we note that the camera integration can be viewed as an analog low-pass filter, which can be analyzed in a similar manner to the high-pass filter in the preceding paragraph. This low-pass filter leads to two effects which have been described previously: a signal loss [28] and an error in Doppler shift estimation [29]. Both of these effects may be understood using the formalism we have described: At higher scanning speeds, the width of the temporal autocorrelation of static tissue is decreased, leading to increased temporal power spectral density width. If the temporal power spectral density width increases, static tissue experiences more attenuation from a temporal low-pass filter, leading to the signal loss described by Yun et. al [28]. Additionally, if the temporal power spectral density is wide, low-pass filtering causes distortion and a frequency estimation error. This distortion should be worse at higher Doppler shift frequencies, as confirmed by Walther et. al [29].

2.5 Velocity estimation

While angiograms can be used to visualize perfusion, quantitative interpretation of OCT angiograms is challenging. The brightness in an angiogram is related only to the amount of scattering from moving sources. Therefore, a decrease in intensity of the angiogram may be due to a decrease in scattering or a decrease in velocity, or a combination of both. Measurements of velocity and flow would therefore improve the utility of OCT for quantitative neuroscience research.

Doppler OCT measurements estimate either the phase slope of the autocorrelation function in the time domain or the center frequency of the power spectral density in the Fourier domain [30]. Our procedure for velocity estimation is an extension of the Kasai autocorrelation algorithm [31]. Velocity estimation was performed after high-pass filtering [16]. The Kasai autocorrelation, shown below, has been used frequently in phase-resolved Doppler OCT measurements [32]

equation m10

In the above equation Ai is the ith complex axial scan. This method can be viewed as a series of two steps. The first step is estimating the value of the autocorrelation function at a single time lag [31].

equation m11

The summation can be performed spatially in the axial and/or transverse direction over a window. Larger windows can be used to reduce noise at the expense of spatial resolution. The second step in the Kasai method is estimation of the phase derivative of the autocorrelation function at zero time lag, as shown below.

equation m12

However, when sampling is dense, the autocorrelation function generally extends over more than two axial scans. Therefore, an extension of the standard approach may be desirable in the case of dense sampling. In this case, it is possible to calculate the complex autocorrelation function over a range of delays (T, 2T, etc.) and use phase information over a range of lags, as opposed to just a lag of one exposure time,

equation m13

In this equation the autocorrelation is estimated at multiple delays. In our procedure, we first estimate the autocorrelation function as described and estimate the Doppler frequency ft,est as the frequency corresponding to the maximum value of the fft of the autocorrelation function estimate. This is equivalent to maximizing the power spectral density estimate:

equation m14

A window of N = 7 was used, and a 256 point fft was used. The estimation was only performed where the estimate of Rt(0) exceeded the noise threshold. The frequency estimate could then be related to the z projection of velocity using the following expression:

equation m15

2.6 Flow calculation

The topography of blood supply to the cortex consists of pial arteries on the surface of the brain which branch to form diving arterioles supplying well-defined cortical capillary beds. The capillary beds are drained by ascending venules, with multiple venules draining a capillary bed supplied by a single arteriole. Ascending venules then converge to form pial veins. There are pathways for lateral and collateral flow in the cortical parenchyma, but, in general, flow in surface vessels should account for all flow in the parenchyma below [33]. By measuring ascending venules at the locations where velocity measurements have the highest signal to noise ratio, it is possible to obtain precise measurements of blood flow. Figure 5(c) shows an en face slice through the Doppler OCT volume at a cortical depth of approximately 50 µm. A pial artery (A) generates both positive and negative Doppler shifts as it undulates relative to the probe beam. Small arterioles (red) and venules (blue) are visible. Total blood flow measurements obtained by venous and arterial flow should theoretically be identical. Veins have less pulsatility, venules are more numerous than arterioles, and measured venous blood velocities fall within the Doppler detection range of our system. Therefore, we chose to measure total venous flow by measuring ascending venules directly before they joined pial vessels. Venules were identified based on the OCT angiogram as well as the en face Doppler images.

A three-dimensional OCT data set consisting of 160 images with 840 axial scans per image over an 800 µm square area was acquired in approximately 7.5 seconds. The three-dimensional OCT data set, after Doppler processing, yields a three-dimensional map of the z-projection of velocity, vz(x,y,z). We assume that the vector velocity field, v(x,y,z), has zero divergence and therefore represents a conservative field. Therefore, the surface integral of v(x,y,z) representing velocity flux, or flow in a vessel, is independent of the surface chosen, provided that the surface bisects the vessel.

equation m16

Typically a surface in the xz or yz plane is chosen for this integral. This is the most obvious choice of plane as OCT images are typically acquired as cross-sections. For instance, in ophthalmology, a dual circumpapillary scanning protocol has been proposed, where two cross-sectional images are used to calculate the angle of the vessel relative to the optic axis [34]. However, this choice of plane has the disadvantage that it requires explicit calculation of the angle of the vessel relative to the optic axis. Figure 5(a) shows a vessel inclined at an angle [var phi] with respect to the optic axis. Referring to this figure, it is possible to calculate the flow by integrating over the yz plane:

equation m17

For this work, we obtain flow from a 3D OCT scan. We choose to perform the surface integration over the en face (xy) plane, which is perpendicular to the measured velocity z component. Figure 5(b) shows the projection of this vessel in the en face plane. Figure 5(c) shows an en face cut through the 3D Doppler OCT volume. Figure 5(d) shows a zoom of the vessel marked by a white square. The area of this vessel in the xy plane is proportional to 1/|cos([var phi])|, as shown in Fig. 5(b), whereas the z projection of the velocity is proportional to cos([var phi]). Therefore, these two effects cancel when the integration is performed in the en face plane. Mathematically, by choosing the surface normal in the –z direction,

equation m18

This procedure does not require explicit calculation of the vessel angle, since the surface normal is parallel to the component of velocity that is measured. Areas for integration were determined by manually choosing a depth and an area in the en face plane for each vessel. The area for integration was chosen to be larger than the vessel cross-section.

2.7 Data processing

All data processing, was performed using compiled Matlab 7.4 code operating on a network cluster consisting of either Dual Opteron 246 (4 GB) or Dual Opteron 248 (16 GB) nodes. Data sets were divided into blocks of approximately 8400 axial scans for processing. Data processing for the angiograms included filtering to remove fixed frequency noise from the CCD line scan camera, spectral shaping, λ-k interpolation using cubic splines, dispersion compensation, Fourier transformation, correction of phase shifts due to bulk axial motion, and high-pass filtering. To display angiograms, a maximum intensity projection along the axial direction was performed. This yielded a two-dimensional en face (xy) image showing vasculature. Data processing for the flow measurements included filtering to remove fixed frequency noise from the CCD camera, spectral shaping, λ-k interpolation using cubic splines, dispersion compensation, Fourier transformation, high-pass filtering, and fitting to determine the phase slope. A series of ten 3D OCT data sets was acquired for the Doppler measurements. To ensure accurate phase measurements, the fast axis of the raster scan was reversed for alternate 3D OCT measurements. Averaging of Doppler OCT data sets acquired with opposite scan directions reduced fixed pattern speckle noise in the Doppler images [35] and helped to compensate for any residual Doppler shift caused by misalignments (See Section 2.3 above).

3. Results

Figure 6(a) shows an OCT angiogram of the rat cortex, obtained by the procedure described in Section 2.4. The transverse resolution was 12.0 µm. Figure 6(b) shows a zoom of an arterial anasthmosis. Figure 6(c) shows an angiogram obtained with a 3.6 µm transverse resolution. The scanning velocity was chosen to yield the same temporal speckle statistics for both resolutions, and the high-pass filter used for both images is shown in Fig. 3(c), 3(d).

The results of absolute flow measurements in the rat cortex are shown in Fig. 7. Venules are marked and numbered on the OCT angiogram [Fig. 7(a)]. Flow measurements at each location along with standard errors over 10 repeated measurements are shown in Fig. 7(b). Quantitative flow measurements were performed in the somatosensory cortex in 4 rats. Because flow measurements are performed by scanning an optical beam over the cortex, the most natural units for the Doppler OCT flow measurements are volume per unit cortical surface area per unit time. Our measurements yielded 1.07 +/− 0.16 x 10−3 mL/mm2/min.

4. Discussion

In this study we present methods for quantification of CBF in the rat cortex. Total retinal blood flow in humans was recently measured for the first time with Doppler OCT [34, 36]. These methods relied on calculation of the vessel angle with respect to the probe beam from two OCT cross-sectional images acquired at different transverse locations. Here we perform three-dimensional imaging of the rat cortex and calculate absolute blood flow by integration of the two dimensional flow profile in the en face (xy) plane. This method does not require explicit calculation of the vessel angle with respect to the probe beam. By limiting our flow measurements to ascending venules it is possible to achieve quantitative flow measurements across the cortex.

Visualization of the capillary network using flow-based contrast alone presents a challenge for several reasons. Flow is frequently less than 1 mm / s and capillaries are predominantly perpendicular to the brain surface, reducing the Doppler shift measured by OCT. Additionally, both static and dynamic scattering can occur within the same voxel. We experimentally and theoretically investigated the effects of the scanning velocity on the bandwidth of speckle fluctuations caused by transverse scanning. Based on our results we optimized scan protocols and processing methods to improve sensitivity to slow flow by drastically reducing the transverse scanning speed to < 1 mm/s. Because the high-pass filtering process enables separation of moving (RBC) and stationary scatterers (vessel wall, smooth muscle, and cells) within a resolution voxel, it is possible to visualize capillaries with sizes at or below the optical resolution. Therefore, the resulting OCT angiograms presented here (Fig. 6) show surprising spatial detail given the transverse resolution of ~12 µm.

4.1 Comparison with autoradiography

Measurements from 4 rats yielded a flow value of 1.07 +/− 0.16 x 10−3 mL/mm2/min in the somatosensory cortex under alpha-chloralose anesthesia. Assuming a cortical thickness of 1.5 to 2.0 mm and a brain tissue density of 1.05 g/mL [37], our results correspond to a blood flow measurement in the somatosensory cortex of 0.51-0.68 mL/g/min. A recent iodo[14C]antipyrine autoradiographic study measured resting CBF in the somatosensory cortex of ~0.6 mL/g/min in similar rats and anesthesia conditions [38], which is consistent with our results. This comparison should not be regarded as proof of agreement between flow measurements by OCT and autoradiography. Rather this comparison suggests that our method yields flow values within a known physiologic range. As heterogeneity of blood flow over the cortex is possible, we expect that our measurement variance can be further reduced by quantifying flow over a larger area (3 mm x 3 mm) of the cortical surface.

4.2 Image acquisition time

3D OCT scanning protocols developed typically require less than 10 seconds. The requirement for fast imaging is dictated by sample motion or data processing and storage limitations. Additionally, conventional wisdom states that axial scans should be spaced less than approximately half the transverse resolution in order to satisfy the Nyquist sampling criteria, which gives the minimum number of voxels required to sample a given volume. For high sensitivity flow imaging, however, a transverse axial scan spacing much denser than the Nyquist criteria is necessary. The time-frequency uncertainty principle states that in order to resolve two frequencies separated by Δf it is necessary to observe for at least Δt ~1/Δf. For the confocal flying spot implementation of OCT, the beam is typically moving in the transverse direction across the sample tissue over time. The spatial bandwidth of the speckle pattern is determined by the transverse resolution (wi), and the relationship between space and time is determined by the velocity of transverse scanning. The coordinate relationships are summarized in the equation below.

equation m19

Therefore slower scan speeds and longer dwell times are required in order to improve velocity sensitivity. Previous work in the brain has used scanning densities of 1000 axial scans over 2.5 mm [39]. Figure 2 shows that scanning at slower speeds leads to a narrower temporal bandwidth of stationary tissue. Therefore, for this study a scan velocity of 0.55 mm/s was used to visualize capillary flow. The high quality angiogram in Fig. 6(a) required 12 minutes of data acquisition. While the acquisition time is comparable to the time required for a two photon imaging stack, for some applications better time resolution is required. The speed limitation of the method presented here arises from the need to scan an optical beam in two dimensions while dwelling on each pixel for a minimum time Δt. An alternative approach is to use a line illumination configuration where a single B scan can be acquired simultaneously [40]. The line could then be scanned and high-pass filtering could be performed along the axis of the scan. Another related approach is to rapidly raster scan a flying spot and perform high-pass filtering along the slow axis of the raster scan [41].

6. Conclusion

In conclusion, we demonstrated Doppler OCT for absolute measurements of CBF and flow-based capillary angiography in the rat cortex. A model for the effect of static tissue on velocity sensitivity was derived and validated. The model was then used to design scanning protocols and algorithms for sensitive 3D Doppler measurements and angiography of the cerebral cortex. Finally, a method of calculating absolute flow that does not require vessel angles was introduced. These measurements do not require exogenous contrast and can be directly compared longitudinally in a single animal over time or cross-sectionally across groups of animals. This capability will greatly aid studies of cerebral pathophysiology and aid in the testing of pharmacological agents in animal models.


We would like to acknowledge scientific contributions and helpful advice from Chao Zhou at MIT, and James Jiang and Alex Cable from Thorlabs. Dr. Iwona Gorczynska’s present address is the Institute of Physics, Nicolaus Copernicus University, Toruń, Poland. We acknowledge support by the National Institutes of Health ( R01-NS057476, P01NS055104, P50NS010828, K99NS067050 and R01-CA075289-13), the Air Force Office of Scientific Research FA9550-07-1-0014, and the Medical Free Electron Laser Program FA9550-07-1-0101.

References and links

1. Girouard H., Iadecola C., “Neurovascular coupling in the normal brain and in hypertension, stroke, and Alzheimer disease,” J. Appl. Physiol. 100(1), 328–335 (2006). doi: 10.1152/japplphysiol.00966.2005. [PubMed] [Cross Ref]
2. Sakurada O., Kennedy C., Jehle J., Brown J. D., Carbin G. L., Sokoloff L., “Measurement of local cerebral blood flow with iodo [14C] antipyrine,” Am. J. Physiol. 234(1), H59–H66 (1978). [PubMed]
3. Calamante F., Thomas D. L., Pell G. S., Wiersma J., Turner R., “Measuring cerebral blood flow using magnetic resonance imaging techniques,” J. Cereb. Blood Flow Metab. 19(7), 701–735 (1999). doi: 10.1097/00004647-199907000-00001. [PubMed] [Cross Ref]
4. Heiss W. D., Graf R., Wienhard K., Löttgen J., Saito R., Fujita T., Rosner G., Wagner R., “Dynamic penumbra demonstrated by sequential multitracer PET after middle cerebral artery occlusion in cats,” J. Cereb. Blood Flow Metab. 14(6), 892–902 (1994). doi: 10.1038/jcbfm.1994.120. [PubMed] [Cross Ref]
5. Kleinfeld D., Mitra P. P., Helmchen F., Denk W., “Fluctuations and stimulus-induced changes in blood flow observed in individual capillaries in layers 2 through 4 of rat neocortex,” Proc. Natl. Acad. Sci. U.S.A. 95(26), 15741–15746 (1998). doi: 10.1073/pnas.95.26.15741. [PubMed] [Cross Ref]
6. Lauritzen M., Fabricius M., “Real time laser-Doppler perfusion imaging of cortical spreading depression in rat neocortex,” Neuroreport 6(9), 1271–1273 (1995). doi: 10.1097/00001756-199506090-00010. [PubMed] [Cross Ref]
7. Fercher A. F., Briers J. D., “Flow visualization by means of single-exposure speckle photography,” Opt. Commun. 37(5), 326–330 (1981). doi: 10.1016/0030-4018(81)90428-4. [Cross Ref]
8. Dunn A. K., Bolay H., Moskowitz M. A., Boas D. A., “Dynamic imaging of cerebral blood flow using laser speckle,” J. Cereb. Blood Flow Metab. 21(3), 195–201 (2001). doi: 10.1097/00004647-200103000-00002. [PubMed] [Cross Ref]
9. Bonner R., Nossal R., “Model for laser Doppler measurements of blood flow in tissue,” Appl. Opt. 20(12), 2097–2107 (1981). doi: 10.1364/AO.20.002097. [PubMed] [Cross Ref]
10. Chen Z., Milner T. E., Dave D., Nelson J. S., “Optical Doppler tomographic imaging of fluid flow velocity in highly scattering media,” Opt. Lett. 22(1), 64–66 (1997). doi: 10.1364/OL.22.000064. [PubMed] [Cross Ref]
11. Izatt J. A., Kulkarni M. D., Yazdanfar S., Barton J. K., Welch A. J., “In vivo bidirectional color Doppler flow imaging of picoliter blood volumes using optical coherence tomography,” Opt. Lett. 22(18), 1439–1441 (1997). doi: 10.1364/OL.22.001439. [PubMed] [Cross Ref]
12. Leitgeb R. A., Schmetterer L., Drexler W., Fercher A. F., Zawadzki R. J., Bajraszewski T., “Real-time assessment of retinal blood flow with ultrafast acquisition by color Doppler Fourier domain optical coherence tomography,” Opt. Express 11(23), 3116–3121 (2003). doi: 10.1364/OE.11.003116. [PubMed] [Cross Ref]
13. White B. R., Pierce M. C., Nassif N., Cense B., Park B. H., Tearney G. J., Bouma B. E., Chen T. C., de Boer J. F., “In vivo dynamic human retinal blood flow imaging using ultra-high-speed spectral domain optical coherence tomography,” Opt. Express 11(25), 3490–3497 (2003). doi: 10.1364/OE.11.003490. [PubMed] [Cross Ref]
14. Izatt J. A., Hee M. R., Owen G. M., Swanson E. A., Fujimoto J. G., “Optical coherence microscopy in scattering media,” Opt. Lett. 19(8), 590–592 (1994). doi: 10.1364/OL.19.000590. [PubMed] [Cross Ref]
15. Wang R. K., Jacques S. L., Ma Z., Hurst S., Hanson S. R., Gruber A., “Three dimensional optical angiography,” Opt. Express 15(7), 4083–4097 (2007). doi: 10.1364/OE.15.004083. [PubMed] [Cross Ref]
16. Ren H., Sun T., MacDonald D. J., Cobb M. J., Li X., “Real-time in vivo blood-flow imaging by moving-scatterer-sensitive spectral-domain optical Doppler tomography,” Opt. Lett. 31(7), 927–929 (2006). doi: 10.1364/OL.31.000927. [PubMed] [Cross Ref]
17. Tao Y. K., Davis A. M., Izatt J. A., “Single-pass volumetric bidirectional blood flow imaging spectral domain optical coherence tomography using a modified Hilbert transform,” Opt. Express 16(16), 12350–12361 (2008). doi: 10.1364/OE.16.012350. [PubMed] [Cross Ref]
18. Wang R. K., An L., “Doppler optical micro-angiography for volumetric imaging of vascular perfusion in vivo,” Opt. Express 17(11), 8926–8940 (2009). doi: 10.1364/OE.17.008926. [PubMed] [Cross Ref]
19. Vakoc B. J., Lanning R. M., Tyrrell J. A., Padera T. P., Bartlett L. A., Stylianopoulos T., Munn L. L., Tearney G. J., Fukumura D., Jain R. K., Bouma B. E., “Three-dimensional microscopy of the tumor microenvironment in vivo using optical frequency domain imaging,” Nat. Med. 15(10), 1219–1223 (2009). doi: 10.1038/nm.1971. [PubMed] [Cross Ref]
20. Srinivasan V. J., Sakadzić S., Gorczynska I., Ruvinskaya S., Wu W., Fujimoto J. G., Boas D. A., “Depth-resolved microscopy of cortical hemodynamics with optical coherence tomography,” Opt. Lett. 34(20), 3086–3088 (2009). doi: 10.1364/OL.34.003086. [PubMed] [Cross Ref]
21. Kleinfeld D., Delaney K. R., “Distributed representation of vibrissa movement in the upper layers of somatosensory cortex revealed with voltage-sensitive dyes,” J. Comp. Neurol. 375(1), 89–108 (1996). doi: 10.1002/(SICI)1096-9861(19961104)375:1<89::AID-CNE6>3.0.CO;2-K. [PubMed] [Cross Ref]
22. Choma M. A., Ellerbee A. K., Yang C. H., Creazzo T. L., Izatt J. A., “Spectral-domain phase microscopy,” Opt. Lett. 30(10), 1162–1164 (2005). doi: 10.1364/OL.30.001162. [PubMed] [Cross Ref]
23. Park B., Pierce M. C., Cense B., Yun S. H., Mujat M., Tearney G., Bouma B., de Boer J., “Real-time fiber-based multi-functional spectral-domain optical coherence tomography at 1.3 microm,” Opt. Express 13(11), 3931–3944 (2005). doi: 10.1364/OPEX.13.003931. [PubMed] [Cross Ref]
24. Yazdanfar S., Yang C., Sarunic M., Izatt J., “Frequency estimation precision in Doppler optical coherence tomography using the Cramer-Rao lower bound,” Opt. Express 13(2), 410–416 (2005). doi: 10.1364/OPEX.13.000410. [PubMed] [Cross Ref]
25. Makita S., Hong Y., Yamanari M., Yatagai T., Yasuno Y., “Optical coherence angiography,” Opt. Express 14(17), 7821–7840 (2006). doi: 10.1364/OE.14.007821. [PubMed] [Cross Ref]
26. Vakoc B. J., Tearney G. J., Bouma B. E., “Statistical properties of phase-decorrelation in phase-resolved Doppler optical coherence tomography,” IEEE Trans. Med. Imaging 28(6), 814–821 (2009). doi: 10.1109/TMI.2009.2012891. [PubMed] [Cross Ref]
27. Schmitt J. M., Xiang S. H., Yung K. M., “Speckle in optical coherence tomography,” J. Biomed. Opt. 4(1), 95–105 (1999). doi: 10.1117/1.429925. [PubMed] [Cross Ref]
28. Yun S. H., Tearney G., de Boer J., Bouma B., “Motion artifacts in optical coherence tomography with frequency-domain ranging,” Opt. Express 12(13), 2977–2998 (2004). doi: 10.1364/OPEX.12.002977. [PMC free article] [PubMed] [Cross Ref]
29. Walther J., Koch E., “Transverse motion as a source of noise and reduced correlation of the Doppler phase shift in spectral domain OCT,” Opt. Express 17(22), 19698–19713 (2009). doi: 10.1364/OE.17.019698. [PubMed] [Cross Ref]
30. Szkulmowski M., Szkulmowska A., Bajraszewski T., Kowalczyk A., Wojtkowski M., “Flow velocity estimation using joint Spectral and Time domain Optical Coherence Tomography,” Opt. Express 16(9), 6008–6025 (2008). doi: 10.1364/OE.16.006008. [PubMed] [Cross Ref]
31. Kasai C., Namekawa K., Koyano A., Omoto R., “Real-Time Two-Dimensional Blood-Flow Imaging Using an Auto-Correlation Technique,” IEEE Trans. Sonics Ultrason. 32, 458–464 (1985).
32. Yang V. X. D., Gordon M. L., Mok A., Zhao Y. H., Chen Z. P., Cobbold R. S. C., Wilson B. C., Vitkin I. A., “Improved phase-resolved optical Doppler tomography using the Kasai velocity estimator and histogram segmentation,” Opt. Commun. 208(4-6), 209–214 (2002). doi: 10.1016/S0030-4018(02)01501-8. [Cross Ref]
33. Nishimura N., Schaffer C. B., Friedman B., Lyden P. D., Kleinfeld D., “Penetrating arterioles are a bottleneck in the perfusion of neocortex,” Proc. Natl. Acad. Sci. U.S.A. 104(1), 365–370 (2007). doi: 10.1073/pnas.0609551104. [PubMed] [Cross Ref]
34. Wang Y., Bower B. A., Izatt J. A., Tan O., Huang D., “In vivo total retinal blood flow measurement by Fourier domain Doppler optical coherence tomography,” J. Biomed. Opt. 12(4), 041215 (2007). doi: 10.1117/1.2772871. [PubMed] [Cross Ref]
35. Wang R. K., Ma Z., “Real-time flow imaging by removing texture pattern artifacts in spectral-domain optical Doppler tomography,” Opt. Lett. 31(20), 3001–3003 (2006). doi: 10.1364/OL.31.003001. [PubMed] [Cross Ref]
36. Michaely R., Bachmann A. H., Villiger M. L., Blatter C., Lasser T., Leitgeb R. A., “Vectorial reconstruction of retinal blood flow in three dimensions measured with high resolution resonant Doppler Fourier domain optical coherence tomography,” J. Biomed. Opt. 12(4), 041213 (2007). doi: 10.1117/1.2771553. [PubMed] [Cross Ref]
37. Pellicer A., Gayá F., Madero R., Quero J., Cabañas F., “Noninvasive continuous monitoring of the effects of head position on brain hemodynamics in ventilated infants,” Pediatrics 109(3), 434–440 (2002). doi: 10.1542/peds.109.3.434. [PubMed] [Cross Ref]
38. Nakao Y., Itoh Y., Kuang T. Y., Cook M., Jehle J., Sokoloff L., “Effects of anesthesia on functional activation of cerebral blood flow and metabolism,” Proc. Natl. Acad. Sci. U.S.A. 98(13), 7593–7598 (2001). doi: 10.1073/pnas.121179898. [PubMed] [Cross Ref]
39. Wang R. K. K., Hurst S., “Mapping of cerebro-vascular blood perfusion in mice with skin and skull intact by Optical Micro-AngioGraphy at 1.3 mum wavelength,” Opt. Express 15(18), 11402–11412 (2007). doi: 10.1364/OE.15.011402. [PubMed] [Cross Ref]
40. Yasuno Y., Endo T., Makita S., Aoki G., Itoh M., Yatagai T., “Three-dimensional line-field Fourier domain optical coherence tomography for in vivo dermatological investigation,” J. Biomed. Opt. 11(1), 014014 (2006). doi: 10.1117/1.2166628. [PubMed] [Cross Ref]
41. Srinivasan V. J., Jiang J. Y., Yaseen M. A., Radhakrishnan H., Wu W., Barry S., Cable A. E., Boas D. A., “Rapid volumetric angiography of cortical microvasculature with optical coherence tomography,” Opt. Lett. 35(1), 43–45 (2010). doi: 10.1364/OL.35.000043. [PubMed] [Cross Ref]

Articles from Optics Express are provided here courtesy of Optical Society of America