1.

Introduction

Noninvasive flow mapping in blood vessels is a field of growing interest. The most common and mature techniques are Doppler ultrasound (DU) and laser Doppler velocimetry^{1, 2} (LDV). Other purely optical techniques such as optical coherence tomography (OCT) and laser speckle imaging (LSI) have also been utilized for blood-flow measurement.^{3, 4} All these techniques are based on scattering of waves from a moving medium and characterization of velocity-dependent features in the scattered waves. In particular, DU and LDV are based on characterization of the Doppler shift in the frequency of the scattered wave. Recently, it was proposed to use the photoacoustic Doppler (PAD) effect for in vivo measurement of flow.^{5, 6, 7, 8, 9, 10, 11} Photoacoustic (PA) imaging in general is based on measurement of the acoustical waves that are generated due to the absorption of modulated light in a tested medium. When the excited region comprises a moving fluid with light-absorbing particles, such as the red blood cells (RBCs) in blood, the generated PA waves exhibit a Doppler shift proportional to the flow velocity.

Vascularly related PA-based measurement techniques are typically based on the optical absorption of hemoglobin in the visible and near-IR spectral ranges, which is significantly higher than that of the plasma and the surrounding tissues. Therefore PA imaging of blood vessels typically provides superior contrast relative to scattering-based techniques such as ultrasonography. Another advantage of PA imaging is its ability to provide functional information, such as the level of hemoglobin oxygenation. This can be accomplished by taking advantage of the spectral differences in the absorption of oxygenated and deoxygenated hemoglobin. In addition to the advantages of the optical excitation, PA techniques also benefit from the advantages of ultrasonic detection, which enables relatively large penetration depth compared to purely optical techniques.¹²

In recent studies of PA flow measurement, four different techniques were described: the first was termed PAD and was based on a continuous sinusoidal modulation of a cw optical source and lock-in detection.^{5, 6} This technique enabled a direct and complete measurement of the Doppler spectrum but did not provide axial position information. The second technique employed the so called optical resolution mode of PA imaging and was based on a pulsed source and cross-correlation of responses of consecutive pulses.^{7, 8} While not providing the complete Doppler spectrum of the imaged medium, this approach enabled fast spatial characterization of velocity via the estimation of the spectrum's mean and variance. In a different study, partial spatial characterization of velocity was attained via direct estimation of the time shift between successive signals.⁹ Another study^{13, 14} demonstrated the implementation of ultrasound M-mode scanning approach for PA flow measurement. Here as well, axial information was attainable but the maximum measurable velocity was limited to ∼1 mm/s.

In this paper, we describe a method that enables simultaneous measurement of axial position and flow as well as complete characterization of the Doppler spectrum over a wide range of mean velocities (3.5 to 200 mm/s). Using the complete Doppler spectrum for estimation of the velocity and its distribution, rather than the spectrum's first two moments, is advantageous as it enables spectral filtering of noise and artifacts such as slow wall movements. It also ensures a more accurate estimation of velocity variability in the presence of turbulence.¹⁵ The inherent ability of the proposed method to measure high velocities, which is a result of the high pulse repetition frequency (PRF), does not trade off with the minimum measurable velocity and may become useful in several applications for mitigating aliasing, as described in Sec. 4. Therefore, this method may be applicable to both low velocities in the microvasculature, where ultrasound often experiences poor SNR, and for high velocities, where ultrasound may suffer from aliasing and the proposed technique is much less restricted.

An optical signal generation method, similar to the one that was used in this work, was reported in details in a recent publication.¹⁰ It was based on tone-burst excitation and was implemented by external electro-optic modulation of a cw laser source. The receiver end of our system, on the other hand, is significantly different from the one described in the previous report. While the previous implementation was based on heterodyne detection via an rf spectrum analyzer, in this work, the output from the transducer was directly digitized and time-gated spectral analysis (also known as short-time Fourier transform) was used to obtain simultaneous position and velocity information. To test the viability of the new method it was used to characterize flow in a vessel phantom with spatially varying lumen. By applying different levels of local external pressure it was possible to emulate flow inhomogeneity and irregularities characteristic of various vascular diseases. Such irregularities may appear in atherosclerosis in large vessels but also in peripheral and retinal microvasculature, for example, due to micorangiopathy.¹⁶ Another experimental test of the method was a quantitative estimation of flow velocities over a wide range of 3.5 to 203 mm/s. This range of velocities renders the proposed method applicable to large blood vessels as well as to microvasculature.

2.

Methods and Experimental Setup

When flowing light-absorbing particles are illuminated by intensity-modulated light, an acoustic wave is generated via the PA effect. Due to the motion of the particles the generated PA signal exhibits a Doppler shift [TeX:] $f_{\rm {PAD}}$ $f_{\mathrm{PAD}}$ , which can be expressed as.^{5, 10}

Similarly to other Doppler-based flow measurements,¹⁷ the PAD spectrum around [TeX:] $f_0 + f_{\rm {PAD}}$ $f_0+f_{\mathrm{PAD}}$ is broadened due to several mechanisms. These mechanisms include the finite dimensions and spatial inhomogeneities of the excitation beam which, respectively, truncates and modulates the response of each particle, the finite size of the ultrasound detector, which leads to a distribution of angles θ [and hence, according to Eq. (1), a distribution of Doppler frequencies] and the diversity in particles velocities in the measured volume. Assuming a laminar flow regime, all these mechanisms lead to broadening, which is linearly dependent on the mean particles velocity.¹⁰

According to the law of continuity, the flow velocity in a vessel with spatially varying lumen attains a maximum value at the point where the diameter of the vessel is at its minimum. At weak stenosis levels, the flow is laminar everywhere and the mean velocity is proportional to [TeX:] ${1 \mathord{\left/ {\vphantom {1 {D^2 }}} \right. \kern-\nulldelimiterspace} {D^2 }}$ $1/\phantom{1D^2}{D}^2$ , where D is the local diameter.¹⁸ In strong stenoses, the flow beyond the stricture may become turbulent and contain eddies similar to the ones created beyond a narrow gorge in a river (Fig. 1). The small eddies are irregular and change rapidly with time. They drift along the flow and add a rapidly varying component to the velocity. In the frequency domain, the presence of a turbulent component broadens the Doppler-shifted spectrum.¹⁹

Fig. 1

The flow through a stenotic tube.

To enable time-resolved spectral characterization of the measured volume, the intensity of the optical excitation was modulated by a tone-burst signal centered at the ultrasonic frequency range. The intensity of the optical excitation can be expressed as

Optically, the experimental setup (Fig. 2) comprised a fiber-coupled tunable laser source (TLS) in the range of 1510 to 1620 nm, a polarization controller, and a lithium-niobate electro-optic modulator (EOM) driven by an arbitrary waveform generator (AWG). The modulated light was amplified using an erbium-doped fiber amplifier (EDFA) and its output was delivered to the sample using a single-mode fiber with a cleaved end. Using a bias-T the EOM was biased to its minimum transmission operation point to reduce the optical power between pulses to minimum. The frequency of the sinusoidal modulation was set to 0.5 MHz, which yielded at the output of the EOM optical pulses with modulated power centered at 1 MHz. The external modulation scheme that was used enabled an accurate and flexible control over the parameters of the excitation such as center frequency, pulse rate, and duty cycle. A similar setup was previously proposed for enhancement of sensitivity and responsivity of PA measurements.²⁰ The specific excitation parameters in the current work were chosen to be [TeX:] $T_r$ $T_r$ = 3 μs, [TeX:] $T_c$ $T_c$ = 70 μs, and N = 14, 286 (the number of bursts in a 1-s sequence) dictating theoretical values for the axial resolution, 4.5 mm; the axial range, 10.5 cm; and a spectral resolution of 1 Hz.

Fig. 2

Experimental setup.

The experimental setup operated in the spectral range around 1550 nm, taking advantage of components available for fiber optic communications. While this setup is an excellent test bed for developing effective excitation waveforms and data analysis algorithms its spectral range has the disadvantage of relatively high optical absorption of water. Therefore, in vivo implementation of the measurement method requires shifting to the wavelength range between 500 and 1000 nm, where low water absorption and high hemoglobin absorption ensure superior contrast. Equivalent setups that can be implemented using direct modulation of laser diodes or cw sources and electro-optic modulators at the desired spectral range are currently under study.

The modulated light, at mean power of 100 mW, was directed through a circular aperture of 4mm diameter to a transparent C-flex tube (inner diameter; 2.4 mm). The tube was coiled to form a loop as illustrated in Fig. 2. The transducer was placed at one side of the loop, with a slight tilt with respect to the horizontal axis, and the gap between the transducer and the tube was filled with an ultrasonic gel to ensure good acoustical coupling. The illumination geometry and the transducer's position and orientation were designed to obtain a relatively small angle between the flow direction and the transducer beam. In this manner, the generated Doppler shift approached its maximum value. The acoustical detection was implemented using an immersion ultrasound transducer with a center frequency of 1 MHz (Olympus IR-0108-S-SU) followed by a low-noise preamplifier (Olympus 5660C). The output of the preamplifier was sampled by a data acquisition device (DAQ: Agilent U2331A), which transferred the digitized data to a PC. The DAQ allowed a maximum sampling rate of 3M samples/s, which limited the modulation frequency to 1.5 MHz and a maximum sequence length of 3M samples, which led to a maximum sequence duration of 1 s, and therefore frequency resolution of 1 Hz. The energy of each burst in this setup was quite low, however, sufficient to obtain a good SNR in the experimental conditions. The SNR can be improved using beam focusing, a higher PRF, and of course optical sources at higher power.

A syringe pump (Fusion-200, Chemyx) with a 50-cc syringe was used to generate flow rates of up to 55 ml/min (in both withdrawal and infusion modes) in the tube, giving a maximum mean velocity of 203 mm/s for the given tube diameter. The fluid used was a suspension of carbon particles with typical particles dimensions of 37 to 149 μm (Activated charcoal C3345, Sigma-Aldrich) in deuterium oxide (D₂O) with volume fraction of approximately 10%. D₂O was chosen as the solvent due to its relatively low absorption (around 1 cm⁻¹) in the experiment spectral range. The tested suspension viscosity is different from that of blood, since the D₂O has a 3 times lower viscosity than blood. The carbon particles should also affect the viscosity. Therefore, we must take into account that performing the same experiment with blood should yield a different Reynolds number and therefore a different extent of turbulence when the flow regime deviates from laminar.

To find the Doppler shift as a function of axial position, the recorded 1-s long periodic sequence was processed by a computerized time-gated spectral analysis procedure (also referred to as a short-time Fourier transform). The sequence was multiplied by a periodic gate sequence of the same length and period (see Fig. 3). The width of each gate was nine samples, which corresponded to 3 μs. The result was Fourier transformed [using a fast Fourier transform (FFT)] and its absolute value was recorded. This calculation was repeated for 70 temporally shifted versions of the gate sequence yielding a 2-D map of frequency and time (or position). The temporal shift between each two versions of the gate sequence was 1 μs. To improve the SNR, maps from eight successive independent measurements of the phantom were averaged to produce the presented experimental maps (Figs. 4 and 6 to 8 in Sec. 3). The time axis was converted to axial position by multiplication by the sound velocity in the matching gel, which was assumed to be approximately 1500 m/s.

Fig. 3

Data processing procedure: examples of two time gates and the resulting spectra in the 2-D map.

Fig. 4

(a) Two-dimensional maps of PSD versus axial position for withdrawal (right map) and infusion (left map) at a rate of 40 ml/min and (b) cross section at distance of 1.4 cm for 40 ml/min: infusion (black) and withdrawal (gray).

Fig. 6

Two-dimensional maps of PSD versus axial position for withdrawal at rate = 40 ml/min without stenosis at beam locations alongs the tube relative to the stenosis of (a) −4.5 mm, (b) −3 mm, (c) −1.5mm, (d) 0 mm (beam centered at stenosis wire), (e) +1.5 mm, and (f) +3 mm.

Fig. 7

Two-dimensional maps of frequency shift versus axial position for withdrawal at rate = 40 ml/min with weak stenosis at beam locations alongs the tube relative to the stenosis of (a) −4.5 mm, (b) −3 mm, (c) −1.5mm, (d) 0 mm (beam centered at stenosis wire), (e) +1.5 mm, and (f) +3 mm.

Fig. 8

Two-dimensional maps of frequency shift versus axial position for withdrawal at rate = 40 ml/min with tight stenosis at beam locations along the tube relative to the stenosis of (a) −4.5 mm, (b) −3 mm, (c) −1.5mm, (d) 0 mm (beam centered at stenosis wire), (e) +1.5 mm, and (f)+3 mm.

Two sets of experiments were performed. In the first set, undisturbed flows at 16 different flow rates between 1 and 55 ml/min in both infusion and withdrawal modes were produced in the tube and the PA responses were recorded. In the second set of experiments, stenosis was induced at the middle of the horizontal section of the tube using a thin metal wire (diameter ∼0.1 mm). The wire was coiled around the tube at different tightening levels: no stenosis; weak stenosis, which reduced the tube diameter by approximately 20%; and strong stenosis, which reduced the tube diameter by approximately 30%. The acoustical signal was recorded for each of these stenosis levels at six different displacements d of the illuminating beam relative to the stenosis. The displacements were 0, ±1.5, ±3, and −4.5 mm, where locations upstream (proximal) to the stenosis are considered negative and locations downstream (distal) to the stenosis are positive. The measurements were taken at withdrawal mode at a rate of 40 ml/min.

3.

Experimental Results

Examples of 2-D maps of PSD versus axial position for withdrawal and infusion at a rate of 40 ml/min (equivalent to a mean velocity of 147 mm/s) for an unperturbed tube are shown in Fig. 4. The color scale describes normalized power in decibels. Also shown are vertical cross sections of the maps at the positions where the observed Doppler shifts are at maximum. These positions correspond to the location of the tube.

The Doppler shift as well as the Doppler broadening at the axial position of the illuminated region is clearly seen. The maps also show a strong signal at the center frequency [TeX:] $f_0$ $f_0$ . This signal is attributed to the response of stationary absorbers, such as the tube and particles that were deposited on the tube inner sidewalls. This background signal, broadened according to the spectral resolution of the DAQ, limited the minimum measurable velocity to 3.5 mm/s (equivalent to a flow rate of 1 ml/min in our setup and a Doppler shift of ∼2 Hz). The maximum measurable velocity was limited to 203 mm/s by the maximum rate of the pump, however, this limitation can be easily alleviated by the use of a stronger pump.

The cross section plots at the location of maximum Doppler shift were used to estimate the mean Doppler shift. This was done by smoothing the response and averaging the frequencies at –1 dB from the peak, after filtering out a ±1 Hz window around [TeX:] $f_0$ $f_0$ to eliminate the stationary responses. The mean particles velocity was calculated from the estimated Doppler shift according to Eq. (1), assuming that the speed of sound in D₂O is [TeX:] $c_A = 1400$ $c_A=1400$  m/s and the effective angle between the flow and transducer is [TeX:] $\theta \approx 25\deg$ $\theta \approx 25\mathrm{deg}$ . The resulting estimated velocities, along with the theoretical velocities as calculated from the flow rate and the tube dimensions, as a function of flow rate are shown in Fig. 5. The velocity resolution is dictated by the spectral resolution (1 Hz) to 1.5 mm/s. Note that the spectral broadening can also be used for velocity estimation, when the flow regime is laminar, as was demonstrated in our previous paper,¹⁰ however, due to the multiple factors contributing to the broadening, this method requires prior calibration.

Fig. 5

Velocity versus flow-rate estimated from Doppler shift (blue with markers at measurement points) and theoretic, based on rate and tube dimensions (dotted red). Inset includes zoom-in on the low velocities section. (Color online only.)

Two-dimensional maps of PSD versus axial position for the three levels of stenosis and for the six values of d are shown in Figs. 6 through 8. The variation in the Doppler shift and the spectral broadening as a function of the displacement from the stenosis is clearly seen from the color maps. To further quantify this variation, the Doppler shift and the spectral broadening at the position of maximum shift for each of the maps was estimated. The shift was found as described earlier and the broadening was determined as the FWHM. The Doppler shift and the spectral broadening as a function of the displacement of the illumination point from the stenosis for the three stenosis levels are shown in Fig. 9. As in the first part of the experiment, a frequency window of [TeX:] $f_0\, \pm 1\,{\rm Hz}$ $f_0\pm 1\mathrm{Hz}$ was filtered out to exclude responses from static absorbers, such as the tube and the metal wire at the stenosis.

Fig. 9

(a) Doppler shift and (b) Doppler spectral broadening versus beam distance relative to stenosis location for three levels of stensois: no stenosis (dashed blue), weak stenosis (dotted green), and tight stenosis (solid red). Calculation was performed on map cross sections at the location of maximum shift. Beam locations proximal to the stenosis are negative and locations distal to the stenosis are positive. (Color online only.)

We see that the Doppler shift and the spectral broadening depend on the stenosis level and the illumination point and their behavior is consistent with the expected flow profile around a stenosis. At the locations proximal to the stenosis, −4.5 mm ⩽ d ⩽ 3 mm, the spectral parameters approach their unperturbed values for both stenosis levels. At the stenosis (− 1.5 mm ⩽ d ⩽ 1.5 mm) the Doppler shift and the spectral broadening increase by a factor of approximately 2 for tight stenosis and ∼1.7 for weak stenosis. Distal to the stenosis (d ⩾ 3 mm) the Doppler shift starts to decrease, however, the spectral broadening remains the same and even slightly increases for the tight stenois. This behavior is attributed to the presence of turbulence, which gives rise to enhanced spectral broadening. The dependence of the spectral parameters on the relative distance to the stenosis is similar to DU spectrum obtained in stenotic vessels.¹⁸ Proximal to the stenosis, a narrow spectrum is observed. At the stenosis, the spectrum remains pretty narrow, but the Doppler shift increases significantly due the high velocity. Distal to the stenosis, a broad spectrum is observed due to turbulent flow. Due to the random nature of turbulent flow, this broadening is an indicator of its presence but it is difficult to use it for flow quantification.

The Doppler shift change at the stenosis, relative to the locations proximal to the stenosis, enables estimation of the stenosis degree. The shift change corresponded to ∼23 and ∼30% decrease in the tube diameter for the weaker and tighter stenosis, respectively, assuming that the Doppler shift ratio is proportional to the inverse square of the tube diameters. Though it was difficult to accurately estimate the exact tube diameter at each test, these results approximately agree with the experimental conditions.

4.

Discussion

PAD imaging of flow has several potential advantages, with respect to DU, in both small and large vessels. In regard to microvasculature it was already pointed out^{5, 6} that PAD has better SNR than DU since the former can be tuned to specifically excite hemoglobin where the latter suffers from low scattering from the RBCs and relatively high background noise from the surrounding tissues. In the case of large vessels and in particular in the presence of stenosis, flow velocities can become very high. In these circumstances, the Doppler shift may exceed PRF/2, resulting in aliasing. Clearly, one simple step to alleviate this problem is to increase the PRF, however, the PRF also determines the measurement range and can not be arbitrarily increased. Typically, an aliasing problem is more severe in deep vessels, such as the aorta, but in some cases stenoses appear also in peripheral vessels (which can potentially be examined using PAD) such as the superficial femoral artery and the popliteal artery.^{21, 22} Interestingly, DU is much more prone to aliasing than PAD. The Doppler shift in DU is twice that of PAD and the PRF in DU must be half the PRF in PAD to impose an identical range limit. This leads to the PAD being 4 times less sensitive to aliasing than DU.

Clearly, the reduced Doppler shift in PAD is a limitation in low velocities where the higher the shift the easier it is to separate between moving and static responses. However, the improved SNR is a significant advantage of PA in microvasculature, as mentioned in previous PAD studies, which might compensate for this drawback. In these cases, the full spectral analysis, which enables filtering of individual frequencies, may help to filter out the background noise. In comparison with flow measurement techniques based on optical scattering, PAD can reach higher penetration depths, while these techniques are limited to a few millimeters.⁶ Another important potential use of PAD is to improve contrast in PA imaging of blood vessels, especially in situations of low SNR. Via spectral filtering the signal from static or slowly moving tissues in the surrounding of the blood vessel can be removed, leading to cleaner images of the vessels inner volumes.

For in vivo measurements, the spectral, temporal, and spatial resolution of the proposed technique should be adapted for the desired application. As expected, there is an inherent trade-off between the spectral and temporal resolution. The temporal resolution in the proposed method is determined by the duration of the sequence that is Fourier transformed and processed to obtain a single frequency versus axial position map. Assuming that the bursts are coherent, the spectral resolution is simply the inverse of this duration. In the experimental demonstration of the method, the flow was temporally constant. Hence, long sequences of 1 s and the corresponding spectral resolution of 1 Hz could be obtained. However, for measurement of pulsatile flow the sequence should be limited to 10 to 20 ms to enable adequate temporal sampling of the varying flow, leading to a spectral resolution of 50 to 100 Hz. This places a limit on the velocity resolution and the minimum measurable velocity, which can be reduced by increasing the modulation frequency. Figure 10 plots the minimum measurable velocity versus the processed sequence duration for several modulation frequencies. We can see that while a modulation of 1 MHz may not be sufficient for pulsatile flow, increasing the modulation frequency to 10 MHz will enable a resolution of 7.5 mm/s, which is reasonable for large vessels where the typical velocities are in the scale of a few tens to a few hundreds of millimeters per second.

Fig. 10

Minimum measurable velocity versus measurement duration (log scale) for different modulation frequencies: 1 MHz (solid gray), 10 MHz (dashed pink), 20 MHz (dotted green), and 50 MHz (solid black). Red dotted vertical line marks the duration limit for pulsatile flow measurement and the red cross marks the conditions of the experimental setup presented in this paper. (Color online only.)

For measurement of low velocities of few millimeters per second, typical of microvasculature, even the resolution obtained for modulation of 50 MHz may be insufficient, therefore a longer sequence is required, preventing tracking of the temporal pulse. However, note that in microvasculature, the pulsatility is significantly damped and steady flow is often assumed due to the low Wormsley number in vessels with small diameters.²³

The axial resolution is limited by the duration of the signal received from a point source, which is dictated by the minimum between the bandwidth of the ultrasonic detector and the pulse bandwidth. In our case, the transducer bandwidth was the limiting factor. Hence, using high-center-frequency transducers and increasing the modulation frequency will improve the axial resolution, however, it will also increase the acoustic attenuation and decrease the penetration depth. In the experimental setup already described, the transducer bandwidth yielded a rather poor resolution of 4.5 mm. To obtain a resolution below 100 μm, the bandwidth should be increased beyond 10 MHz. The lateral resolution can be improved as well by focusing the optical beam.

5.

Conclusions

Tone-burst excitation and time-gated spectral analysis PAD were used for the first time to characterize flow with and without spatial variations in a laboratory vessel model. Accurate and flexible control over the parameters of the tone-burst modulated excitation, which enabled the simultaneous measurement of full spectral content and axial position, was achieved via an EOM driven by an AWG. The PA responses were detected by an ultrasonic transducer and were digitized by a DAQ device. The digitized data was subjected to time-gated spectral analysis, which yielded 2-D maps of Doppler spectrum versus axial position.

Using the setup, velocities from 3.5 to 203 mm/s were accurately estimated from the Doppler shift. The introduction of stenoses in the measured vessel produced spatial flow variations that could be identified in the 2-D maps. The characteristic increase in the Doppler shift and spectral broadening at the stenosis were clearly observed as well as the asymmetry between spectra taken from two different sides of the stenosis. The same spectral attributes are hallmarks of stenotic flow in DU. The increase in the Doppler shift enabled quantification of the stenosis strength in addition to its localization. The spatial resolution achieved in the experiment was limited by the bandwidth of the DAQ and the ultrasonic transducer to 4.5 mm, but is expected to improve significantly with the implementation of a faster DAQ and a transducer with higher central frequency and bandwidth.

The results shown in this paper demonstrate the feasibility of the PAD method using tone-burst excitation for characterization of flow and identification of flow spatial irregularities. The method is very promising and possesses several potential advantages over other scattering-based optical or acoustical techniques, such as additional biological information, superior contrast in blood vessels, and reduced aliasing limit in high velocities. The burst excitation and the signal analysis method, which yield the full Doppler spectrum, enable noise reduction by means of spectral filtering and should be more immune to potential estimation errors relative to calculation of the mean and variance of the velocity alone. Adaptation of the setup to a spectral range that is more suitable for in vivo measurements, as well as resolution improvement, are currently under research.