## 1.

## Introduction

Blood flow in an organ or tissues through capillaries represents an important factor for diagnosing pathological conditions, including heart failure,^{1} liver cirrhosis,^{2} and pancreatitis,^{3} and for evaluating the physiological condition, such as in renal transplants.^{4} Moreover, antiangiogenesis therapies for tumors have recently received wide attention in animal models used in preclinical research.^{5, 6} Treatment effects can be evaluated by longitudinal observation of tumor angiogenesis, which can be achieved by measuring the capillary flow in the tumor.^{7, 8} The flow rate in a considerable small region of interest (ROI) is commonly measured by monitoring the concentrations of exogenous indicators, which is also known as the contrast-specific measurement method.^{9, 10, 11} Recently, we have successfully developed contrast-specific methods for photoacoustic measurements of flow using gold nanorods as the contrast agent.^{12, 13, 14, 15} Although these methods have been shown to be effective, they are designed for measuring flows in large vessels, where the ROI is inside the vessel. In this work, we extend the method to vessels smaller than the photoacoustic sample volume using a high-frame-rate photoacoustic imaging system. One major advantage of such a 2-D photoacoustic flow measurement system is that it provides both anatomical and flow information.

Contrast-specific flow measurements have been developed using various imaging modalities, such as computed tomography,^{2, 3, 16} magnetic resonance imaging,^{17, 18} and ultrasound.^{19} Contrast agents are used as flow indicators, with their dilution monitored as a function of time. In general, contrast agents are materials that are injected into the human body to enhance the signal-to-noise ratio (SNR) of blood vessels in the ROI. Although concentration is the parameter of interest, most contrast-specific methods measure the signal intensity, implicitly assuming that the intensity is linearly proportional to the concentration.^{17, 19} The measured signal intensity over time is also known as the time-intensity curve (TIC). The flow rate is determined based on the indicator-dilution theory,^{20, 21, 22} which typically employs a mixing chamber that relates the two TICs obtained at the inflow and outflow of the mixing chamber. The mixing chamber is often modeled as a linear and time-invariant (LTI) system, and therefore the relation between the input and output of the mixing chamber can be described by a transfer function as a wash-out analysis.^{23} The transfer function is approximated as an exponential function with a time constant determined by the flow rate. Therefore, given both the input and output TICs, deconvolution can be applied to estimate the transfer function, from which the flow rate can be determined.^{24} The applicability of the deconvolution-based method is determined by the validity of the LTI system model and the accuracy of deconvolution.^{24}

Wash-in flow estimation methods have been proposed that overcome the difficulties associated with deconvolution-based methods.^{15, 19, 25} This method requires the concentration of the contrast agent to be reduced by “destroying” indicators in the ROI. In ultrasonic imaging, for example, a microbubble-based contrast agent can be destroyed by irradiation with acoustic pulses of sufficient energy,^{25, 26} with the subsequent temporal changes in the concentration determined by the replenishment rate of the contrast agent (i.e., new microbubbles flowing into the ROI). This replenishment process is related to the local flow rate.^{15} Because deconvolution is not required, the wash-in method is generally more stable and accurate.

We have previously developed wash-in flow measurements for photoacoustic imaging^{14, 15} using gold nanorods as the photoacoustic contrast agent. Gold nanorods are biomedically compatible materials that offer a strong photoacoustic response in the visible to near-infrared region.^{27} The wavelength at which a nanorod exhibits peak optical absorption is determined by its aspect ratio (defined as the length of the major axis divided by that of the minor axis). The shape of a nanorod can be altered by laser pulses with sufficient energy, which is known as the nanorod-to-nanosphere shape transition.^{28} The laser-induced shape transition of gold nanorods leads to a reduction of the aspect ratio (often being transformed into nanospheres) and causes a downshift of the absorption peak. In practice, the nanorods are not destroyed, but instead their shapes are changed, which effectively reduces the concentration of the nanorods with the original aspect ratio, and the associated shift in the peak absorption wavelength is monitored. The wash-in flow measurements used in ultrasonic imaging are still applicable.

The purpose of this work is to extend the previously proposed method from flow estimation from a single ROI to 2-D mean flow velocity measurements. Achieving this requires a 2-D photoacoustic imaging system with an adequate frame rate, the design and construction of which is also described here. The accuracy of flow rate measurements made by applying this system to flow phantoms is also reported.

This work is organized as follows. The principles of flow rate measurement based on the time-intensity method are described in Sec. 2. Section 3 describes the high-frame-rate imaging system and the experimental apparatus of the flow rate measurement. Section 4 describes the measurement results from flow phantoms. Section 5 discusses the results and the range of flow rates measurable by this system, and conclusions are drawn in Sec. 6.

## 2.

## Principles of Photoacoustic Wash-in Flow Measurement

The method used to estimate the flow rate consists of both the destruction (i.e., depletion) and the replenishment (i.e., wash-in) of the nanorods. Laser pulses with a sufficiently high energy are used to cause rod-to-sphere shape transitions, and also to monitor the reduced concentration of gold nanorods in the ROI. The process is illustrated in more detail in Fig. 1 . A steady state is first reached before the application of laser pulses $(t<{t}_{0})$ after the injection of the gold nanorods. A laser pulse transforms gold nanorods into nanorods with a smaller aspect ratio (possibly into nanospheres) at $t={t}_{0-}$ , and the detected photoacoustic-signal intensity is reduced $(t={t}_{0+})$ . After that, new nanorods (without shape transitions) flow into the ROI at a rate determined by the flow rate $({t}_{0}<t<{t}_{1})$ , and the concentration of gold nanorods with the original aspect ratio increases in the ROI. Another laser pulse irradiates after a pulse repetition interval (PRI) and causes another period of destruction and the following replenishment of the gold nanorods. After several periods of the destruction-replenishment process, the concentration of nanorods in the ROI will reach a steady state if the number of incoming new nanorods approximates that of the destroyed nanorods.

The TIC measured in this period of time can be used to calculate the flow rate by fitting it to a curve derived from the replenishment model. The model used was adapted and modified from our previous work.^{15} Assuming that the laser beam is much larger than the diameter of the capillaries inside a flow region and that the beam profile of the ultrasound transducer along the
$y$
axis is
$u\left(y\right)$
, the effective concentration measured by this transducer is

## Eq. 1

$${n}^{\prime}\left(t\right)=\frac{\int n(y;t)u\left(y\right)\mathrm{d}y}{\int u\left(y\right)\mathrm{d}y}=\int n(y;t)\widehat{u}\left(y\right)\mathrm{d}y,$$In the replenishment phase, predicting $n(y;t)$ for $t>{t}_{0}$ based on $n(y;{t}_{0})$ is required to estimate flow parameters, for which the initial condition is $n(y;{t}_{0})={n}_{0}-{n}_{0}l\left(y\right)$ , where $l\left(y\right)$ denotes the shape of the laser beam along the $y$ axis and is typically a low-pass function of $y$ . Assuming that the flow system is linear and shift invariant gives

## Eq. 2

$$n(y;t)=\int g(y-{y}_{0};t,{t}_{0})n({y}_{0};{t}_{0})\mathrm{d}{y}_{0}={n}_{0}-{n}_{0}g(y;t,{t}_{0})\otimes l\left(y\right),$$## Eq. 4

$$h\left[v(t-{t}_{0})\right]=\int [g(y;t,{t}_{0})\otimes l\left(y\right)]{u}^{\prime}\left(y\right)\mathrm{d}y.$$^{29}the density of label, i.e., nanoparticles, along the flow direction $x$ after a bolus injection can be represented as $F(x,\sigma )=(1\u2215\sigma \sqrt{2\pi})\mathrm{exp}[-{(x-\overline{x})}^{2}\u22152{\sigma}^{2}]$ , where $\sigma $ and $\overline{x}$ denote the random variable and the distance after the injection, respectively. In this case, the bolus injection, which is produced by the rapid laser destruction, is determined at $t={t}_{0}$ , and $\sigma =\sqrt{2D(t-{t}_{0})}$ at the central position of $\overline{x}=v(t-{t}_{0})$ . Therefore,

## Eq. 5

$$g(y-{y}_{0};t,{t}_{0})=\frac{1}{\sqrt{4\pi D(t-{t}_{0})}}\phantom{\rule{0.2em}{0ex}}\mathrm{exp}\{-\frac{{[y-{y}_{0}-v(t-{t}_{0})]}^{2}}{4D(t-{t}_{0})}\},$$^{29, 30}Assuming that both the laser and the ultrasound beam profile are Gaussian distributed in $y$ [i.e., ${u}^{\prime}\left(y\right)=\mathrm{exp}(-{y}^{2}\u2215{\sigma}_{u}^{2})\u2215\sqrt{\pi {\sigma}_{u}^{2}}$ and $l\left(y\right)={l}_{0}\phantom{\rule{0.2em}{0ex}}\mathrm{exp}(-{y}^{2}\u2215{\sigma}_{l}^{2})\u2215\sqrt{\pi {\sigma}_{l}^{2}}$ , where ${\sigma}_{u}$ and ${\sigma}_{l}$ are standard deviations of the ultrasound and the laser beam profile, respectively], Eq. 4 can be solved as

## Eq. 6

$$h\left[v(t-{t}_{0})\right]=\frac{{l}_{0}}{{\left\{\pi [4D(t-{t}_{0})+{\sigma}_{u}^{2}+{\sigma}_{l}^{2}]\right\}}^{1\u22152}}\phantom{\rule{0.2em}{0ex}}\mathrm{exp}\{-\frac{{\left[v(t-{t}_{0})\right]}^{2}}{4D(t-{t}_{0})+{\sigma}_{u}^{2}+{\sigma}_{l}^{2}}\}.$$## Eq. 8

$${n}^{\prime}\left(t\right)\cong {n}_{0}\{1-c\phantom{\rule{0.2em}{0ex}}\mathrm{exp}[-\lambda v(t-{t}_{0})]\}={n}_{0}\{1-c\phantom{\rule{0.2em}{0ex}}\mathrm{exp}[-\beta (t-{t}_{0})]\},$$## Eq. 9

$$\text{\hspace{0.17em}}{\phantom{\mid}1-{n}^{\prime}(t-{t}_{0})\mid}_{t={({t}_{0}+k\bullet \mathrm{PRI})}^{+}}\u2215{n}_{0}\cong c{\phantom{\mid}\mathrm{exp}[-\beta (t-{t}_{0})]\mid}_{t={({t}_{0}+k\bullet \mathrm{PRI})}^{+}},$$## Eq. 10

$$\text{\hspace{0.17em}}{\phantom{\mid}1-{n}^{\prime}(t-{t}_{0})\mid}_{t={[{t}_{0}+(k+1)\bullet \mathrm{PRI}]}^{-}}\u2215{n}_{0}\cong c{\phantom{\mid}\mathrm{exp}[-\beta (t-{t}_{0})]\mid}_{t={[{t}_{0}+(k+1)\bullet \mathrm{PRI}]}^{-}}.$$## Eq. 11

$$\text{\hspace{0.17em}}{\phantom{\mid}1-{n}^{\prime}(t-{t}_{0})\mid}_{t={[{t}_{0}+(k+1)\bullet \mathrm{PRI}]}^{-}}\u2215{n}_{0}\cong \mathrm{exp}(-\beta \bullet \mathrm{PRI})[{\phantom{\mid}1-{n}^{\prime}(t-{t}_{0})\mid}_{t={({t}_{0}+k\bullet \mathrm{PRI})}^{+}}\u2215{n}_{0}].$$On the other hand, in the destruction phases, it is noticeable that the replenishment during the destruction phase is negligible because the pulse duration is on the order of nanoseconds and is much smaller than the length of the PRI. In these periods, destruction mainly dominates the behavior of the concentration of nanorods. The concentration decreases asymptotically to a constant level and can be expressed as:

## Eq. 12

$$\text{\hspace{0.17em}}{\phantom{\mid}{n}^{\prime}(t-{t}_{0})\mid}_{t={({t}_{0}+k\bullet \mathrm{PRI})}^{-}}\cong ({n}_{0}-{n}_{\infty}){r}^{k}+{n}_{\infty}\phantom{\rule{1em}{0ex}}\text{for}\phantom{\rule{0.3em}{0ex}}k=0,1,2,\dots ,$$## Eq. 13

$$\text{\hspace{0.17em}}{\phantom{\mid}{n}^{\prime}(t-{t}_{0})\mid}_{t={({t}_{0}+k\bullet \mathrm{PRI})}^{+}}\cong r{\phantom{\mid}{n}^{\prime}(t-{t}_{0})\mid}_{t={({t}_{0}+k\bullet \mathrm{PRI})}^{-}}.$$Finally, we can obtain the concentration of the nanorods by combining the destruction phase in Eq. 13 with the replenishment phase in Eq. 11:

## Eq. 14

$$\text{\hspace{0.17em}}{\phantom{\mid}1-{n}^{\prime}(t-{t}_{0})\mid}_{t={[{t}_{0}+(k+1)\bullet \mathrm{PRI}]}^{-}}\u2215{n}_{0}\cong \mathrm{exp}(-\beta \bullet \mathrm{PRI})[1-r{\phantom{\mid}{n}^{\prime}(t-{t}_{0})\mid}_{t={({t}_{0}+k\bullet \mathrm{PRI})}^{-}}\u2215{n}_{0}].$$## Eq. 16

$${n}^{\u2033}\left(k\right)=\left(\frac{1-s}{1-rs}\right)+[1-\left(\frac{1-s}{1-rs}\right)]{\left(rs\right)}^{k}.$$Finally, the mean velocity can be obtained from the value of the fitted parameter:

According to Wei,^{19}$\lambda $ is inversely proportional to the width of the detection region (i.e., the width of laser beam $E=1\u2215\lambda $ ). Therefore, the flow rate can be calculated from the fitted parameter according to

## 3.

## Experimental Setup

## 3.1.

### High-Frame-Rate Photoacoustic Imaging System

The photoacoustic technique combines laser irradiation and acoustic detection. It offers advantages of both a higher optical contrast than conventional ultrasound and a better penetration depth as compared with all-optical imaging techniques such as microabsorption spectrometry and dark-field microscopy. The capabilities of photoacoustic imaging have been demonstrated in many applications, including functional imaging of rat brain,^{31} breast tumor detection,^{32} and molecular imaging.^{33} Nanoparticles have also been used as contrast agents and molecular probes.^{34, 35} Generally, when a tissue is irradiated by an incident laser pulse, absorption of the laser energy leads to a rapid temperature rise, thermal expansion, and the concomitant generation of broadband ultrasound waves.^{36} These waves can be used to estimate optical properties (i.e., the optical absorption) of the tissue.

According to the wash-in flow measurement model described in Sec. 2, flow rate measurements require the optical absorption of nanorods to be monitored over a considerable period of time. During the replenishment period, absorption changes have to be measured over a short time interval (i.e., requiring a high sampling rate). We have previously demonstrated that a sampling rate of
$15\phantom{\rule{0.3em}{0ex}}\mathrm{Hz}$
allows flow velocities from
$0.35\phantom{\rule{0.3em}{0ex}}\text{to}\phantom{\rule{0.3em}{0ex}}2.8\phantom{\rule{0.3em}{0ex}}\mathrm{mm}\u2215\mathrm{s}$
to be measured in a single sample volume.^{14}

The 2-D photoacoustic imaging system used in this study employed an array transducer instead of a single-crystal transducer to increase the frame rate by allowing data to be acquired simultaneously from multiple channels. This high-frame-rate photoacoustic imaging system consisted of a fiber laser system, an ultrasonic digital phased array system (DiPhAS), which is often used in research applications, and a custom-design photoacoustic probe, as shown in Fig. 2 . A Q-switched Nd:YAG laser (LS-2132U, Lotis TII, Belarus) operating at $1064\phantom{\rule{0.3em}{0ex}}\mathrm{nm}$ with a pulse duration of $8\phantom{\rule{0.3em}{0ex}}\mathrm{ns}$ was used for laser irradiation. The $5\text{-}\mathrm{mm}$ -diam beam from this laser was split and guided onto two multimode fiber light guides (LG-L30-6-H-1500-F-1, Taiwan Fiber Optics, Taiwan) using a half-reflectance beamsplitter (BS1-1064-50, CVI, NM). The output beams from the light guides were then focused by two cylindrical lenses positioned $25\phantom{\rule{0.3em}{0ex}}\mathrm{mm}$ from the light-guide surfaces, resulting in two irradiated zones of $30\times 0.8\phantom{\rule{0.3em}{0ex}}\mathrm{mm}$ . Both the light guides and the focusing lenses were mounted on a custom-design holder as part of the photoacoustic probe for confocal alignment. A motorized lens wheel (FW102, Thorlabs, New Jersey) placed between the laser and the beamsplitter was used to switch the laser energy.

Acoustic waves induced in the irradiated volume were detected by a 128-channel ultrasonic linear array (L6, Sound Technology, Pennsylvania) that was also confocally mounted between the two light guides in the photoacoustic probe, as shown in Fig. 2b. The transducer elements had a pitch of $0.3\phantom{\rule{0.3em}{0ex}}\mathrm{mm}$ , a $5\text{-}\mathrm{mm}$ elevational width, and a center frequency of $5\phantom{\rule{0.3em}{0ex}}\mathrm{MHz}$ with an 82% bandwidth. Reflecting foil (Ho Yan Tape, Taiwan) with a thickness of $9\phantom{\rule{0.3em}{0ex}}\mu \mathrm{m}$ was attached to the surface of the transducer to block backscattered laser irradiation from reaching the transducer surface. The DiPhAS was employed to amplify and digitize the detected rf array signals. It contained 64 transmitting and receiving channels and high-speed multiplexers to allow the use of a transducer array with up to 192 channels, although in this study the multiplexers were not used. The rf signals from the 64 transducer elements sent to the receiver were amplified by up to $80\phantom{\rule{0.3em}{0ex}}\mathrm{dB}$ and then digitized by analog-to-digital converters with $12\text{-}\text{bit}$ precision. The data were sampled at $40\phantom{\rule{0.3em}{0ex}}\text{Msamples}\u2215\mathrm{s}$ , and the on-board memory allowed 512 data samples per channel to be stored before they were transferred to a personal computer (PC). The DiPhAS allowed array data from 64 channels to be simultaneously acquired and transferred every $4\phantom{\rule{0.3em}{0ex}}\mathrm{ms}$ , giving a frame rate of up to $250\phantom{\rule{0.3em}{0ex}}\mathrm{Hz}$ . The raw rf data were transferred through a high-speed digital I/Q card (PCI-7300A, ADlink, Taiwan) to the PC, and dynamic focusing and image reconstruction were performed off-line. The DiPhAS and laser were synchronized using a programmable logic device (EPM3064A, Altera, California), with the laser triggered by the DiPhAS. The actual frame rate of the system was limited to $15\phantom{\rule{0.3em}{0ex}}\mathrm{Hz}$ due to the maximum PRF of the laser that we used.

In this study, 2-D photoacoustic images captured from the same flow region over the replenishment period were used to evaluate the concentration change of gold nanorods. ROIs were chosen to include the flow area of interest, such as capillaries or organs. The mean photoacoustic intensities in these ROIs were calculated at each time step as a data point in the TIC. Therefore, there is a tradeoff between the SNR and the size of the ROI (which determines the spatial resolution).

## 3.2.

### Flow Phantom

Figure 3 shows a schematic diagram of the flow measurement setup. A sample was made of chicken breast tissue, in which polyethylene tubing (Intramedic™ 427411, Sparks, Maryland) with a inner diameter of $580\phantom{\rule{0.3em}{0ex}}\mu \mathrm{m}$ was embedded at a depth of $2\phantom{\rule{0.3em}{0ex}}\mathrm{mm}$ . A volume of $1\phantom{\rule{0.3em}{0ex}}\mathrm{ml}$ of human blood and $55\phantom{\rule{0.3em}{0ex}}\mu \mathrm{l}$ of gold nanorods with an absorption peak at $985\phantom{\rule{0.3em}{0ex}}\mathrm{nm}$ at a concentration of $3.6\phantom{\rule{0.3em}{0ex}}\mathrm{nM}$ were mixed and injected into the tubing with a $1\text{-}\mathrm{ml}$ standard syringe. The absorption peak of the nanorods was measured using a spectrophotometer (V-570, Jasco, Japan). The syringe was controlled by an infusion pump (KDS 100, Montreal, Canada) to create mean flow velocities ranging from $0.125\phantom{\rule{0.3em}{0ex}}\text{to}\phantom{\rule{0.3em}{0ex}}2\phantom{\rule{0.3em}{0ex}}\mathrm{mm}\u2215\mathrm{s}$ . The theoretical flow velocity is the mean value of the entire flow region. It was calculated from the pumping rate of the infusion pump with the cross sectional area size of flow region (i.e., the tubing). Before these experiments, the pumping rate has been tested by measuring the weight of the pumped water. The photoacoustic probe was positioned $13\phantom{\rule{0.3em}{0ex}}\mathrm{mm}$ above the tube, and cross sectional images were captured. TICs were achieved by summing the photoacoustic intensities from within a $0.6\times 1.2\text{-}\mathrm{mm}$ ROI $(8\times 4\phantom{\rule{0.3em}{0ex}}\text{pixels}$ in the $z$ and $x$ axes, respectively) at each time point.

A pilot experiment indicated that applying laser irradiation at energy densities of $82.4\phantom{\rule{0.3em}{0ex}}\mathrm{mJ}\u2215{\mathrm{cm}}^{2}$ was suitable for the destruction of the nanorods in the tube. Figure 4 shows that the image intensity decayed within $4\phantom{\rule{0.3em}{0ex}}\mathrm{s}$ . Figure 5 shows the visible-to-infrared spectra of the nanorods before and after the applications of the destruction pulses as measured using the spectrophotometer, with the gray rectangle highlighting the large absorbance difference at $1064\phantom{\rule{0.3em}{0ex}}\mathrm{nm}$ (the wavelength of the laser irradiation). Figure 6 demonstrates the linearity of the relation between the concentration of the nanorods and the photoacoustic intensities before the flow experiments (for an energy density of $6.45\phantom{\rule{0.3em}{0ex}}\mathrm{mJ}\u2215{\mathrm{cm}}^{2}$ ), and thus the validity of monitoring nanorod concentration by measuring photoacoustic signals.

## 4.

## Experimental Results

Figure 7 shows focused cross sectional images of the flow phantom. In this figure, the positions of the chicken breast tissue and the tubing show consistent with their location. The tubing is at a depth of $2\phantom{\rule{0.3em}{0ex}}\mathrm{mm}$ from the chicken breast surface, a fact that the laser energy irradiates the tubing after propagating through light scattering and absorbing biological tissue. In this case, the destruction of the gold nanorods from the outer tissue sample becomes difficult, since the laser intensity decays along the propagation depth due to the light attenuation. Images captured during the measurement $(t=0\phantom{\rule{0.3em}{0ex}}\text{to}\phantom{\rule{0.3em}{0ex}}16.5\phantom{\rule{0.3em}{0ex}}\mathrm{s})$ are shown in Figs. 8a, 8b, 8c , displayed with a dynamic range of $-10\phantom{\rule{0.3em}{0ex}}\mathrm{dB}$ . These images display a depth of $2.5\phantom{\rule{0.3em}{0ex}}\mathrm{mm}$ and a width of $5\phantom{\rule{0.3em}{0ex}}\mathrm{mm}$ , and were normalized by subtracting the background image captured without gold nanorods. Therefore, these images show only the photoacoustic signals from the nanorods. The laser energy was set at $82.4\phantom{\rule{0.3em}{0ex}}\mathrm{mJ}\u2215{\mathrm{cm}}^{2}$ (i.e., at $t=0\phantom{\rule{0.3em}{0ex}}\mathrm{s}$ ). It is clear that the image intensity decreases with time, indicating that the gold nanorods inside the tube underwent shape changes due to the strong incident laser energy that reduced the optical absorption at the laser wavelength. A steady state was reached after continuous laser irradiation at the same energy level for $16\phantom{\rule{0.3em}{0ex}}\mathrm{s}$ , at which point the number of nanorods that underwent shape transition was roughly the same as the number of new nanorods flowing into the ROI. A square region with a size of $0.6\times 1.2\phantom{\rule{0.3em}{0ex}}\mathrm{mm}$ was selected as the ROI, as indicated by the dashed box in Fig. 8a. The mean image intensities within the ROI were summed to create the TIC. Figure 9 shows the normalized TICs measured at flow rates of 0.125, 0.25, 1, and $2\phantom{\rule{0.3em}{0ex}}\mathrm{mm}\u2215\mathrm{s}$ , which clearly indicate that the replenishment rate of the nanorods increases with the flow rate, as described by Eq. 17.

The TICs were then calculated by fitting, which was performed using MATLAB (The MathWorks, Natick, Massachusetts). The TICs were measured five times in each case, from which mean and standard deviation (STD) values were calculated. In Eq. 18, flow rate $v$ is proportional to the product of fitting parameter $\beta $ (i.e., rate constant) and the width of irradiated volume $E$ . Figure 10 shows the mean flow rate calculated for irradiated zones with widths ranging from $0.2\phantom{\rule{0.3em}{0ex}}\text{to}\phantom{\rule{0.3em}{0ex}}1.0\phantom{\rule{0.3em}{0ex}}\mathrm{mm}$ . The results calculated with a large $E$ (i.e., a broad irradiated zone along the $y$ axis) are overestimations, whereas a small $E$ results in underestimations of the flow velocity. It is clear from the figure that using $E=0.5\phantom{\rule{0.3em}{0ex}}\mathrm{mm}$ gives the best agreement with the theoretical flow rate. However, the width of the irradiated volume of this system is $0.8\phantom{\rule{0.3em}{0ex}}\mathrm{mm}$ (which was measured by using a knife-edge method).

A comparison between the experimental results using $E=0.8\phantom{\rule{0.3em}{0ex}}\mathrm{mm}$ and the actual flow velocities is shown in Fig. 11 . Because the ROI covers nearly the entire cross section of the tube, the measured results can be viewed as the mean flow velocity inside the sample volume. It is obvious that most of the flow velocities are overestimated. The linear regression curve between the measured velocities $\left(y\right)$ and the actual velocities $\left(x\right)$ is described by $y=1.68x-0.12$ , which indicates that the measured flow velocity is higher than the actual one. The correlation coefficient between the measured flow velocities and their linear regression fit was 0.83. Figure 11 also shows that the STDs increase with the mean values, with Fig. 12 plotting the normalized STD-to-mean ratios (the average value is 31.3%). This relationship may be attributable to the SNR of TICs. In high flow velocities, replenishment becomes more effective, which means that the steady state (in which the number of incoming new nanorods approximates to that of the destroyed nanorods) is easier to approach, as shown in Fig. 9. Therefore, fitting such low contrast curves without sufficient SNR leads to measurement errors. This problem can be solved if a larger ROI is chosen to increase the SNR.

## 5.

## Discussion

The results presented here show that the wash-in flow estimation method can be implemented utilizing images captured by a high-frame-rate photoacoustic imaging system. However, two issues need to be addressed before making *in-vivo* measurements. First, our experiments were conducted with a measured laser beam width. However, the width of the irradiated zone in the tissue is difficult to evaluate when the laser pulse propagates through a strongly scattering medium, and an incorrect width leads to measurement errors, as shown in Fig. 10. Therefore, the difference between the assumed and actual width is a critical problem in this method if absolute measurements are required. Nonetheless, relative measurements can be made even in the presence of errors in
$E$
. Second, the measurable range of the flow rate is determined by both the frame rate (i.e., sampling rate) and the width of the irradiated zone in the ROI. There is therefore a tradeoff between the elevational resolution (i.e., the resolution along the axis perpendicular to the cross sectional image) and the measurable flow rate. On the other hand, actual flow rates ranging from
$0.5\phantom{\rule{0.3em}{0ex}}\text{to}\phantom{\rule{0.3em}{0ex}}10\phantom{\rule{0.3em}{0ex}}\mathrm{mm}\u2215\mathrm{s}$
in capillaries smaller than
$200\phantom{\rule{0.3em}{0ex}}\mu \mathrm{m}$
have been verified using intravital microscopy.^{28} Therefore, such applications of this method require a system with a broad sampling range. In addition, preclinical investigations with small-animal models require a better spatial resolution for identifying their anatomic structures.

The rate constant is defined as the ratio of the flow velocity to the width of the irradiated zone, and represents the replenishment rate inside the irradiated volume (see Fig. 1). A large rate constant means that the concentration of new nanoparticles inside the irradiated zone rapidly reaches a steady state. Therefore, either increasing the flow velocity or decreasing the width of the irradiated zone increases the rate constant. Based on this principle, the flow velocity can be determined from the width of the irradiated zone and the sampling rate.

To investigate the relation between the rate constant $\beta $ and the sampling rate, TICs were calculated for various rate constants ranging from $0.2\phantom{\rule{0.3em}{0ex}}\text{to}\phantom{\rule{0.3em}{0ex}}10\u2215\mathrm{s}$ using Eq. 16. The simulated TICs were sampled at rates ranging from $1\phantom{\rule{0.3em}{0ex}}\text{to}\phantom{\rule{0.3em}{0ex}}50\phantom{\rule{0.3em}{0ex}}\text{samples}\u2215\mathrm{s}$ . These TICs were curve fitted for rate-constant estimation (as the last step in flow estimation). Comparisons between the measured rate constants with the predetermined values are presented in Fig. 13 . In this figure, most of the resulting ratios are more than unity, meaning that the measured rate constant is 10 to 20% higher than the determined one. This could be a possible reason of the overestimations shown in Fig. 11. On the other hand, it is also noticeable that the resulting ratios are less than unity when the sampling rate is insufficient, meaning that the measured flow rate is lower than the determined one. This may be due to a high-flow velocity in an irradiated zone, resulting in a large rate constant and thus a short replenishment period. With an irradiated zone of constant width and a fixed frame rate, underestimation is worse when the flow velocity is higher (i.e., $\text{higher}\beta $ ). In Fig. 13, the system frame rate of $15\phantom{\rule{0.3em}{0ex}}\mathrm{Hz}$ is indicated by the vertical dashed line. The flow velocity corresponding to rate constants of $0.2\u2215\mathrm{s}$ , $1\u2215\mathrm{s}$ , and $5\u2215\mathrm{s}$ are 0.16, 0.8, and $4\phantom{\rule{0.3em}{0ex}}\mathrm{mm}\u2215\mathrm{s}$ , respectively (assuming an irradiation width of $0.8\phantom{\rule{0.3em}{0ex}}\mathrm{mm}$ ).

In our previous work, we also reported a two-energy wash-in flow estimation method. This method contains a sequence of laser pulses to destroy gold nanorods, and a following sequence with lower laser energy to monitor the replenishment.^{15} However, lower laser energy propagated through turbid biological tissues cannot provide sufficient SNR for the replenishment measurement. Therefore, in strong scattering and absorbing tissue, it is much easier to achieve the flow estimation method in single energy than that in double energy.

Although color Doppler techniques have been routinely used in diagnostic ultrasound, there have been no reliable flow estimation methods available for photoacoustic imaging other than the proposed method in our study. It is therefore the purpose of this study to demonstrate the feasibility of our approach. Current color Doppler techniques estimate blood flow velocities by calculating the phase of the autocorrelation function of the baseband Doppler signals.^{37, 38} However, such phase-sensitive techniques are not applicable in photoacoustic imaging. In addition, the phase measurement is sensitive to noise and susceptible to tissue motion. Particularly for low velocity flows, design of the clutter filters becomes challenging, if not impossible.^{39} Note that the proposed method is less affected by noise because contrast agents (gold nanoparticles) are used. The estimation can also be performed more reliably because the data are collected over a longer period of time and signal averaging is generally applied.

## 6.

## Conclusions

We experimentally demonstrate the feasibility of nanorod-based flow measurement based on a high-frame-rate photoacoustic imaging system, which currently has a frame rate up to
$15\phantom{\rule{0.3em}{0ex}}\mathrm{Hz}$
. The linearity between the nanorod concentration and photoacoustic intensities is verified. The results of the flow experiments show that flow velocities ranging from
$0.125\phantom{\rule{0.3em}{0ex}}\text{to}\phantom{\rule{0.3em}{0ex}}2\phantom{\rule{0.3em}{0ex}}\mathrm{mm}\u2215\mathrm{s}$
can be measured, with an average normalized STD of 31.3%. The measurable flow rate of this method is limited by the PRF of the laser system. The feasibility of *in-vivo* studies is also been discussed. To apply this flow measurement method to preclinical research involving small-animal models, one future direction is to build a photoacoustic imaging system equipped with a transducer array operating at a higher frequency (to improve the imaging spatial resolution), and a laser with a higher PRF (to provide a high frame rate).

## Acknowledgments

The authors gratefully acknowledge Chung-Ren Chris Wang and Kuei Chen Pao for helpful comments and providing the gold nanorods used in this work. This work was financially supported by the National Science Council under grants NSC-94-2213-E-002-115 and NSC-94-2120-M-002-004, and by the National Taiwan University Nano Center for Science and Technology and National Health Research Institutes.