1. Introduction

Thermal Diffusion Flowmetry (TDF) (also called Heat Clearance Method or Thermal Clearance Method) is a longstanding technique for measuring blood flow or blood perfusion in living tissues [1–4]. The principles of the methods were first conceived by Gibbs as early as 1933 [1] and many variations and modifications have been proposed since. The common principle of all the variants is that the physical variables blood flow and thermal diffusion in a tissue are closely related and the former can be inferred by measurement of the latter. Typically, temperature transients and/or gradients are induced in a volume of interest and the temporal and/or spatial temperature variations which follow are recorded and analyzed. The analysis is usually based on a heat diffusion model which describes the conduction and convection of heat in the measured volume.

Photoacoustic (PA) imaging and spectroscopy is based on measurement of the acoustical waves which are generated due to the absorption of modulated light in a tested medium. PA is widely known for its excellent performance in vascular imaging [5–7]. The main reason for that is the relatively high absorption of hemoglobin in the visible and near-IR spectral regions. Photoacoustics enables high resolution imaging at depths where most purely optical imaging techniques are not functional due to excessive optical scattering [5]. One very useful attribute of PA imaging is its ability to determine the oxygen saturation level of the blood (sO₂) [5,6]. This is typically implemented by using multispectral PA excitation.

Recently it was demonstrated that the efficiency of the PA effect strongly depends on temperature [8,9] and is consequently affected by variations in flow [10]. In this paper we further investigate the PA dependence on temperature and flow and describe a photoacoustic version of TDF which is based on this dependence. In contrast with conventional TDF methods, in which the temperature is measured via a direct thermal contact with the tissue, we utilize the temperature dependence of the efficiency of PA generation for remote monitoring of the blood's temperature. Moreover, the same optics which delivers the PA optical excitation can be used for photothermal (PT) induction of the heat which is needed for implementing the thermal diffusion measurement. This approach offers a number of potential advantages over the traditional methods. One advantage is that with no thermal contact and since the light is absorbed primarily by hemoglobin, the technique selectively probes the blood and not the surrounding tissues. A second advantage is the optional high spatial resolution of the method, facilitated by optical as well as acoustical focusing, which enables measurement of the flow in a specific blood vessel whereas conventional TDF is typically limited to measurement of the global flow in a volume containing multiple vessels. Another very important advantage is the ability to use multispectral PA excitation or PT excitation for determining sO₂ simultaneously with the blood flow. In this way it will be possible to obtain the metabolic rate of oxygen which is an important physical parameter in many clinical situations [11–13].

While there were numerous studies of vascular PA imaging, only a few dealt with characterization of the blood flow. Most of them were based on the PA Doppler (PAD) effect [14–18]. PAD was demonstrated in-vivo in optical resolution mode, when the beam was focused to a size close to a single red blood cell, limiting the penetration depth to less than 1 mm [14,15]. Implementation of vascular PAD in the acoustical resolution mode, however, is rather challenging since the Doppler effect requires heterogeneity of the moving medium. When the spatial separation between the red blood cells is small relative to the acoustical wavelength, the blood seems homogeneous and the effect is suppressed [14,19]. The method described in this paper is not based on the Doppler effect but rather on heat transfer. Therefore, it does not require high spatial resolution and can be implemented at relatively deep tissues.

2. Theory

Thermal Diffusion Flowmetry (TDF) is based on the principle that the tissue temperature variations in a Volume Under Test (VUT) during some measurement time, $\tau =\tau (r,t)$ $\left\{0<t<t_{\text{meas}};r\in \text{VUT\hspace{0.17em}}\right\}$, in response to thermal excitation, are closely related to the flow in the volume and can be used for its evaluation. There are many ways to induce the thermal excitation and to monitor the resulting temperature variations. In this study we propose the use of the photothermal (PT) effect for thermal excitation and the photoacoustic (PA) effect for monitoring the temperature. The quantitative description of the method is divided into two parts. In the first part we will introduce the equations which govern the spatial and temporal behavior of $\tau (r,t)$ in the presence of PT excitation, thermal diffusion and flow-dependent heat convection. In the second part we will describe the effect of the photothermally generated temperature variations on the PA properties of the medium and how this can be utilized for PA monitoring of $\tau (r,t)$and consequently of the flow.

2.1. Photothermal excitation of the VUT

The temperature in a volume of tissue subjected to photothermal excitation is described by a heat equation [20]:

where τ is the deviation from the baseline temperature, α is the thermal diffusivity, ρ is the density, $C_p$ is the heat capacity at constant pressure, $\overrightarrow{u}$ is the velocity of the blood in the tissue and $H_{\text{PT}}$ is the PT heat energy deposited in the medium per unit time per unit volume. For simplicity, metabolic heat generation was neglected. Generally $H_{\text{PT}}(r,t)$ can be written as a product of a space-dependent function and a time-dependent function as $H_{\text{PT}}(r,t)=H_r^{\text{PT}}(r)H_t^{\text{PT}}(t)$. It is also convenient to use normalization where the area under $H_t^{\text{PT}}(t)$is 1.

The temporal variations of τ are determined by the three terms on the right hand side (RHS) of the equation. The first term describes thermal conduction. The second term describes the heat which is added or removed from a given location due to heat convection. This term exists only where $\overrightarrow{u}\ne 0$ and is closely related to the blood flow. The heat induced in the medium due to the absorbed light is described by the third term.

For an excitation pulse which satisfies the thermal confinement condition [6] (but not necessarily the stress confinement condition) $H_{\text{PT}}(r,t)$ can be expressed as $H_{\text{PT}}(r,t)=H_r^{\text{PT}}(r)\delta (t)$. In this case the problem can be formulated as an initial value problem:

with initial conditions

In general, the solution of Eq. (2) for the initial conditions in Eq. (3) and given boundary conditions can be found numerically. Here we concentrate on a simplified lumped solution which facilitates the estimation of $\overrightarrow{u}$ from measurement of the temperature.

2.2. Lumped model

To gain further physical insight regarding the temporal behavior of τ in response to impulse excitation we consider a simplified model, in which the VUT is considered a “lump” and its thermal state is described by a single parameter [21]. We define the lumped temperature of the system as a weighted average over the temperature in the volume:

where $\widehat{f}\left(r\right)$ is currently an arbitrary weight function of r, normalized such that its integral over the VUT is 1. It is introduced here for use in subsequent paragraphs. Using Eq. (2) we obtain

where $\overline{\tau }(r,t)\equiv \tau (r,t)/\tilde{\tau }(t)$.

This representation is particularly helpful in situations where the coefficient of $\tilde{\tau }(t)$ in the RHS of Eq. (5) is independent of time. In this case the decay of $\tilde{\tau }(t)$ following an impulse will be exponential with a time constant $\left(t_{\text{conduction}}^{-1}+t_{\text{convection}}^{-1}\right){}^{-1}$, where

An important observation that will be used to analyze the experimental results is that in the case of uniform flow $t_{\text{convection}}$ becomes inversely proportional to the velocity. In many practical situations a lumped model is used even if the characteristic decay times are not strictly constants. The suitability of the model depends on the thermal and geometrical properties of the VUT and its accuracy can be evaluated empirically. The usefulness of the lumped model in formalizing the thermal dynamics of a blood vessel phantom and its potential role in flow measurement is demonstrated in Section 4.

A simple but compelling description of the lumped model can be made in terms of an RC circuit (Fig. 1 ). In these terms the illuminated volume is conceived as a thermal “capacitor” $C_{\text{thermal}}$ which is charged by the PT source. The absorbed energy is dissipated through two parallel resistors: $R_{\text{conduction}}$ and $R_{\text{convection}}$. The temporal behavior of the thermal circuit is determined by the time constant $t_{\text{eff}}=\left(t_{\text{conduction}}^{-1}+t_{\text{convection}}^{-1}\right){}^{-1}$, where $t_{\text{conduction}}=C_{\text{thermal}}{R}_{\text{conduction}}$ and $t_{\text{convection}}=C_{\text{thermal}}{R}_{\text{convection}}$. In the case of zero flow (stationary blood) the convection resistance is infinite and all the energy is dissipated through $R_{\text{conduction}}$. As the flow increases, the convection resistance decreases and more energy is dissipated through $R_{\text{convection}}$, until, for very high flows, it completely shortens the capacitor and prohibits its charging. Moreover, in light of the definition of $t_{\text{conduction}}$ in Eq. (6) it can be seen that if the velocity is uniform the convection resistance is inversely proportional to the velocity and can be expressed as $R_{\text{convection}}={\tilde{R}}_{\text{convection}}/\left|\overrightarrow{u}\right|$.

 

Fig. 1 RC circuit description of the lumped model.

Download Full Size | PDF

2.3. Photoacoustic monitoring of the temperature variations

The generation and propagation of pressure waves in a fluid due to the photoacoustic effect, in the so-called thermal-confinement regime, can be described by the following inhomogeneous wave equation [22]:

where $p(r,t)$is the pressure field, c is the velocity of sound, $\Gamma \equiv \beta c^2/C_p$ is the Grüneisen coefficient, β is the volumetric thermal expansion coefficient, $C_p$ is the specific heat and $H_{\text{PA}}(r,t)$ is the PA heat energy deposited in the medium per unit time per unit volume. Here as well the heat rate, $H_{\text{PA}}(r,t)$, can be written as a product of a space-dependent function and a time-dependent function as $H_{\text{PA}}(r,t)=H_r^{\text{PA}}(r)H_t^{\text{PA}}(t)$. Again, the area under $H_t^{\text{PA}}(t)$ is normalized to 1.

To use PA for temperature monitoring we take advantage of the temperature dependence of Γ which, for limited temperature deviations, can be approximated by [8,9]

where T is the temperature, $T_0$ is the baseline temperature, τ is the deviation from $T_0$, ${\Gamma }_0$ is the Grüneisen coefficient at $T_0$ and b is the slope of the dependence curve. As described in Section 4, we experimentally found that in blood the relative change in Γ, $\left[\Delta \Gamma \left(T\right)/{\Gamma }_0\right]/\tau =b/{\Gamma }_0$, is approximately 3.5% per degree Celsius in the range $25°\le T\le 35°$. This is similar to the value measured in water and is attributed mainly to the increase in the thermal expansion coefficient [20,23].

In the present form Eq. (7) describes the situation where Γ is constant. In the presence of PT excitation Γ becomes a function of both space and time due to the spatial and temporal dependence of $\tau (r,t)$. Since Γ determines the fraction of deposited heat energy which is converted to pressure at each point in space at a given time, the Poisson solution to Eq. (7) in the case of an optical impulse at $t=t^{\prime }$ becomes

Here $A(t)$ denotes a surface where $\left|r-r^{\prime }\right|=ct$. Note that the Poisson solution at a given position,$\widehat{p}(r_{\text{0}},t-t^{\prime },t^{\prime })$, is no longer a function of $t-t^{\prime }$ alone but rather it depends on the time of excitation $t^{\prime }$ and $H_r^{\text{PA}}(r)$ is weighted according to the spatial dependence of Γ. It is also assumed that Γ changes very slowly with respect to the characteristic PA response time so it can be assumed constant during the response. For an excitation waveform other than an impulse the pressure field can be obtained by

2.4. Temperature monitoring using sinusoidal PA excitation

Consider the case where the PA excitation is sinusoidal, namely,

where ${\omega }_{\text{PA}}$ is the radial frequency of the modulation and $t_{\text{total}}$ is a normalization constant which takes into account the fact that in reality the excitation has a finite duration. Substituting in Eq. (10) and retaining only the term centered at $+{\omega }_{\text{PA}}$ yields

To proceed, we substitute the exact form of $\widehat{p}(r,t^{″},t-t^{″})$ from Eq. (9) and the expression for Γ from Eq. (8) and arrive after several mathematical steps to the following simple result (see Appendix A for details):

where we assumed that ${\Gamma }_0$ is spatially uniform, the observation point is at the origin, $A{e}^{j{\omega }_{\text{PA}}t}$ is the response obtained from Eq. (12) in the absence of photothermal excitation and

where

It can be seen that the PT excitation leads to Amplitude Modulation (AM) of the otherwise purely sinusoidal PA response. The carrier amplitude, A, can be found by turning off the PT excitation. With A and the relative change in the Grüneisen coefficient known, the temperature modulation term can be found from the measured PA response according to

For an impulse PT excitation, in the cases where the lumped model applies, $\tilde{\tau }(t)$ takes an exponential form according to Eqs. (5) and (6):

The Fourier transform of Eq. (17) yields the Modulation Frequency Response (MFR):

It can be seen that in the lumped model approximation the photothermal MFR takes the form of a complex Lorentzian. The applicability of this result to our experimental conditions is demonstrated in Section 4. By fitting the measured MFR to a Lorentzian, $t_{\text{eff}}$ was found and used for estimation of the velocity.

3. Experimental setup and measurements

The experimental setup (Fig. 2 ) comprised a pair of 830 nm Laser Diodes (LD's), each with a maximum peak power of 200 mW. The outputs of the LD's were polarization-combined and the resulting beam was coupled to a multimode fiber. The intensities of both LD's were sinusoidally modulated by direct modulation. One LD was modulated at a frequency of $f_{PA}=0.5\text{\hspace{0.17em}MHz}$ to generate a PA signal. To characterize the MFR the other LD was sinusoidally modulated at discrete frequencies in the range of $f_{PT}\in \left[0.1,20\right]\text{\hspace{0.17em}Hz}$. This choice of a frequency-domain characterization of the modulation response rather than a time-domain approach based on a thermal impulse was made to enable a higher signal to noise ratio. The unconnected end of the fiber was brought in close proximity to a vessel phantom comprising a transparent Tygon tube (1mm inner diameter) filled with defibrinated sheep blood. The size of the illuminated volume was roughly 0.75 mm³. The central part of the tube was immersed in water while one end was attached to a syringe pump and the other end to a container.

 

Fig. 2 The experimental setup.

Download Full Size | PDF

Detection of the PA signal was performed using a non-focused ultrasonic immersion transducer with a center frequency of 0.5 MHz (Olympus IR-0008-P) followed by a low-noise preamplifier (Olympus 5660C). The amplified signal was sampled using a digital acquisition card and recorded by a PC. The sampling rate was 1.25 MS/s and the acquisition time was 20s per data set leading to spectral resolution of 0.05 Hz. Each measurement of PA response was repeated twice and the average response was calculated. A Fourier transform of the PA response yielded a peak at $f_{PA}$ and modulation sidebands at $f_{PA}\pm f_{PT}$.

As a preliminary stage, the temperature dependence of the PA response in stationary blood was characterized. This was performed by a simultaneous measurement of the blood's temperature and its PA response. To measure the temperature, a thermocouple was inserted into the tube and brought in close proximity to the illuminated volume. The position of the thermocouple was determined as a compromise between the need to place it as close as possible to the VUT center and the attempt to avoid PA generation by its direct illumination. To minimize the inaccuracy, the readings were taken only after sufficient time when the transient component in the temperature signal decayed and the system reached a steady state. It should be noted that implementation of the PA-TDF method does not require a thermocouple nor does it use the value of the Grüneisen temperature dependence that was estimated from the thermocouple measurements.

To experimentally test the PA-TDF method, the PA modulation frequency response was measured for several flow rates by changing the PT modulation frequency in the range of $f_{PT}\in \left[0.1,20\right]\text{\hspace{0.17em}Hz}$. The amplitudes of the two PT sidebands were averaged for each value of $f_{\text{PT}}$ to generate the modulation frequency response (MFR). This was repeated for both stationary blood and flowing blood at 5 different flow rates between 0.05 and 1 ml/min. The corresponding blood velocities were between 1.06 and 21.2 mm/sec for the tube that was used.

For each blood velocity u, a normalized Lorentzian was fitted to the measured MFR according to Eq. (18). This procedure yielded the effective time constant as a function of velocity, $t_{\text{eff}}\left(u\right)$. The conduction characteristic time, $t_{\text{conduction}}^{}$, was obtained from the stationary measurement $t_{\text{conduction}}^{}=t_{\text{eff}}^{}\left(u=0\right)$. With $t_{\text{eff}}$ and $t_{\text{conduction}}^{}$ known, $t_{\text{convection}}^{}$ could be calculated according to $t_{\text{convection}}={\left(t_{\text{eff}}^{-1}-t_{\text{conduction}}^{-1}\right)}^{-1}$. As described in the next section, in our tested phantom $t_{\text{convection}}^{}$ was found to be inversely proportional to u. The parameter of proportion, having units of length, can be interpreted as a measure to the spatial dimension of the VUT in the direction of the flow.

4. Experimental results

The direct measurement of temperature using the intra-tube thermocouple enabled characterization of the temperature-dependence of the PA response. Temperature variations were induced in the measurement volume by sinusoidal PT modulation at a frequency of 0.1 Hz. To scan a wide temperature range the peak power of the photothermal-LD was tuned between 0 and 200 mW in increments of 25 mW. At each level of peak power the readings of the thermocouple, after a transient time, oscillated in correspondence with the PT optical excitation. We denote the amplitude of these temperature oscillations as ${\tau }_{\text{AC}}$ [C°] and the deviation of their mean value from the baseline temperature (the thermocouple output in the absence of photothermal excitation) as ${\tau }_{\text{DC}}$ [C°]. In parallel with the PT excitation, the volume was illuminated by the photoacoustic-LD whose power was sinusoidally modulated at a frequency of 0.5 MHz. The presence of PT excitation led to modulation of the generated PA signal. The recorded PA signal was spectrally analyzed and its spectral peaks were identified and recorded. An example of the PA response spectrum in the presence of PT excitation at frequency of 2 Hz is plotted in Fig. 3 . The modulation sidebands due to the PT excitation at $f_{\text{PA}}\pm f_{\text{PT}}$ are evident.

 

Fig. 3 The spectrum of a PA response with PT modulation at 2Hz.

Download Full Size | PDF

The dependence of the central PA peak on ${\tau }_{\text{DC}}$ normalized by its value without thermal excitation, is plotted in Fig. 4(a) . An estimate to this ratio can be found from Eq. (13) by approximating $\overline{\tilde{\tau }\left(t\right)}\approx {\tau }_{\text{DC}}$ (where $\overline{x\left(t\right)}$ denotes the time average of $x\left(t\right)$). It can then be shown that the graph in Fig. 4(a) should follow $1+\left(b/{\Gamma }_0\right){\tau }_{\text{DC}}$. The linear fit in Fig. 4(a) implies a slope of $b/{\Gamma }_0=0.034[1/°C]$. This result, which manifests a temperature dependence of 3.4% per degree Celsius of the Grüneisen coefficient of blood, is in close agreement with the value obtained for water which was around 3.85% [1/C°] [20,23]. The dependence of the modulation peak on ${\tau }_{\text{AC}}$ together with a linear fit is plotted in Fig. 4(b). In this case the anticipated slope is $\left(b/{\Gamma }_0\right)/2$ since the value is taken from a modulation sideband whose amplitude is half the amplitude of the modulating sinusoidal signal. The fitted slope, however, was 0.015 which is slightly lower than the expected value. The deviation is attributed to the very long response time of the thermocouple, which leads to underestimation of ${\tau }_{\text{AC}}$ even in the relatively low modulation frequency of 0.1 Hz. The thermal properties of the VUT, as well as the excitation frequency, affect the ratio between ${\tau }_{\text{AC}}$ and ${\tau }_{\text{DC}}$, and as will be seen in the next paragraphs, this ratio should increase with the increase of flow.

 

Fig. 4 Simultaneous measurements of temperature and PA response at different PT modulation power: (a) PA amplitude of carrier frequency, normalized by amplitude without PT modulation, vs. the average temperature τ_DC. Experimental data marked by squares, linear fit in dashed black. (b) The amplitude of 0.1 Hz modulation peak, normalized by the carrier amplitude without PT modulation, vs. τ_AC.

Download Full Size | PDF

Normalized modulation frequency responses for stationary blood ($q_0=0\text{\hspace{0.17em}ml}/\text{min}$) and 5 flow rates of $q_{1-5}\in \left[0.05,0.1,0.2,0.5,1\right]$ are shown in Fig. 5 . For the tube that was used these flow rates corresponded to the following blood velocities: $u_{0-5}\in \left[0,1.06,2.1,4.2,10.6,21.2\right]$ [mm/s]. All responses were normalized so that their value at $f_{\text{PT}}=0.1\text{\hspace{0.17em}Hz}$ is 1 and they are given in dB scale. In the lumped model approximation the MFR's follows a Lorentzian behavior. In accord with Eq. (18) the normalized responses were fitted to $\left[1+{\left(2\pi 0.1t_{\text{eff}}\right)}^2\right]{}^{1/2}/\left[1+{\left(2\pi f_{\text{PT}}{t}_{\text{eff}}\right)}^2\right]{}^{1/2}$. The suitability of the lumped model for the experimental phantom can be evaluated from the Lorentzian fits in Fig. 5. The agreement is generally good for the four lowest velocities but rather large deviations can be observed in the plots of the two highest velocities. These deviations are attributed to the distribution of the fluid velocity in the tube. As the flow rates increases, the velocity profile, which follows a parabolic shape in a laminar flow regime, becomes more and more peaked and the approximation of uniform velocity, which is implicitly assumed in the lumped model, becomes increasingly less accurate. Clearly, to describe the MFR for non-uniform flow the lumped model should be modified. It is shown below, however, that despite the discrepancy in the MFR, the estimated velocities for non-uniform flow fit well with the expected mean velocities.

 

Fig. 5 Normalized PT-PA modulation frequency responses for stationary (blue) and velocities of 1.06 (pink), 2.1 (green), 4.2 (red), 10.6 (black) and 21.2 (gray) mm/s. Solid lines with markers are experimental data. Dotted lines are Lorentzian fits.

Download Full Size | PDF

The best fit time-constants were, as expected, a decreasing function of the blood velocity: $t_{\text{eff}}^{\text{0-5}}\in \left[0.838,0.405,0.208,0.129,0.054,0.029\right]$ [s]. Using $t_{\text{convection}}={\left(t_{\text{eff}}^{-1}-t_{\text{conduction}}^{-1}\right)}^{-1}$ (where $t_{\text{conduction}}^{}=t_{\text{eff}}^{}\left(u=0\right)$) yielded $t_{\text{convection}}^{\text{0-5}}\in \left[\infty ,0.784,0.277,0.153,0.058,0.03\right]$ [s]. As can be seen in Fig. 6 , these results manifested a linear relation between $t_{\text{convection}}^{-1}$ and the fluid's velocity. The optimal proportion parameter, which we denote as $l_{\text{VUT}}^{}$, was found by minimizing the error between the “real velocity,” u, and the “estimated velocity,” $u_{\text{estimated}}^{}=l_{\text{VUT}}^{}/t_{\text{convection}}$. The error was defined as ${\displaystyle {\sum }_k\left[u_{}^k-l_{\text{VUT}}^{}/t_{\text{convection}}^k\right]{}^2}$ ($k=\mathrm{0..5}$). The physical interpretation of $l_{\text{VUT}}^{}$ (and the origin of its notation) is obtained by recalling that $t_{\text{convection}}^{}$ is the time it takes for the temperature to return to its steady state value, after a short heating pulse, in the absence of heat conduction. In this case the only heat-clearance route is flow and a steady state is reached as soon as the heated fluid leaves the VUT. Hence, $l_{\text{VUT}}\approx u\cdot t_{\text{convection}}$ can be interpreted as a length parameter describing the dimension of the VUT in the direction of the flow. The value of $l_{\text{VUT}}^{}$ which minimized the error between $u_{\text{estimated}}^{}$ and u (for all measured velocities) was found to be 0.62 mm,, which is in general agreement with the Full Width Half Max (FWHM) of the optical beam. It should be noted, however, that the VUT does not have sharp boundaries and analytical determination of $l_{\text{VUT}}^{}$ from the spatial distribution of the PT excitation requires additional analysis. Interestingly, it can be seen that the linear dependence extends even to velocities where the Lorentzian fit showed relatively large deviations from the measured frequency responses. The maximum error between the estimated and the real velocity was 0.3 mm/sec. Also plotted in Fig. 6 (green circles) are velocities that were estimated while ignoring thermal conduction, namely, using ${\tilde{u}}_{\text{estimated}}^{}=l_{\text{VUT}}^{}/t_{\text{eff}}$ with the best-fit $l_{\text{VUT}}^{}$ for this case. As can be seen in Fig. 6 this resulted in increased errors, in particular in the case of low velocities. The maximum error in this case was 0.75 mm/sec.

 

Fig. 6 Estimated velocity ($=l_{\text{VUT}}^{}/t_{\text{convection}}$) vs. real velocity: estimation based on the complete lumped model (red circles). Estimation based on a lumped model which ignores thermal conduction (green circles). Estimation based on only 2 PT modulation frequencies: 1 Hz and 15 Hz (blue x). Dotted black line indicates the 45° line. Inset: zoom-in on the range of low velocities.

Download Full Size | PDF

Another graph in Fig. 6 (blue crosses) shows the estimated velocity based on data from two PT modulation frequencies only (1 and 15 Hz). This can lead to a significant reduction in measurement time (<2 s) at the expense of increased estimation error (in particular in the high velocities).

The system’s sensitivity to temperature modulation can be estimated based on the measured value of $b/{\Gamma }_0$ and the photoacoustic SNR at the carrier frequency. Using Eq. (13) the ratio between the PA response at the modulation sideband $f_{\text{PA}}\pm f_{\text{PT}}$ and the PA response at the carrier frequency $f_{\text{PA}}$ can be approximated as

where it was assumed that $\left(b/{\Gamma }_0\right){\tau }_{\text{DC}}<<1$. Therefore, the minimum detectable temperature oscillation amplitude can be estimated by

with $SNR(f_{\text{PA}})$ defined as the ratio between $p(f_{\text{PA}})$to the noise power at $f_{\text{PA}}\pm f_{\text{PT}}$. In Fig. 3 for example, $SNR(f_{\text{PA}})\approx 74dB$, which leads to ${\tau }_{\text{AC\_min}}\approx 0.01$°C. Note, however, that this sensitivity was obtained when the PA excitation was at maximum power. Any decrease in the PA excitation power will result in an increase in the minimum detectable temperature oscillations.

Finally, the system’s ability to operate at lower input powers was tested. For PT modulation frequency of 0.1 Hz, PT power of 40 mW, PA power of 30 mW and stationary blood, a PT modulation peak of ~5dB above noise level was obtained. These levels give a rough lower bound to the required optical power. It corresponded to ${\tau }_{\text{DC}}$ of 0.9°C (out of which about 0.5°C was contributed by the PT excitation) and to ${\tau }_{\text{AC}}$ of about 0.1C°. As described in Section 3, the readings of the thermocouple possessed some inherent inaccuracy. By waiting sufficient time before taking a reading (until the transient component in the temperature signal decayed) and carefully positioning the thermocouple these inaccuracies were minimized. With these precautions it is estimated that the total increase in temperature did not exceed 2°C. Such temperature increments are comparable with the ones which are induced in other TDF methods, and are well below the threshold of thermal damage [24]. Moreover, the proposed measurement technique can potentially be performed using pulsed PA excitation, provided that the pulse repetition frequency is significantly higher than the photothermal modulation frequency. In this case it is expected that a much lower mean PA input power will be required.

5. Discussion

The experimental results described above provided a first proof of concept to Photoacoustic Thermal Diffusion Flowmetry. They demonstrated, albeit in-vitro, the potential of the method to become a viable tool in vascular characterization. The main advantage of the proposed method over PA Doppler flowmetry is its ability to operate in the acoustic-resolution regime (rather than in the optical-resolution photoacoustic-microscopy mode), which is challenging for PAD due to the homogeneity of the blood. The ability to operate in the acoustic-resolution regime leads to a much bigger measurement depth and simplifies the implementation of the measurement.

It was found that reasonable estimation of the blood velocity, in particular in the case of relatively high velocities, can potentially be obtained without a need for calibration. The only needed parameter is $l_{\text{VUT}}^{}$. While this parameter was obtained from a minimization process in our experiment, it may be possible to estimate it based on the dimensions of the optical beam and the estimation of the optical and thermal propagation in the tissue. For low velocities, where heat conduction plays a significant role, more accurate results can be obtained if the conduction time constant is used to correct the otherwise overestimated velocities. The conduction time can be found via a reference measurement of stationary blood. In-vivo, this can be implemented by a temporary occlusion of the vessel.

Another practical consideration is the duration of the measurement. It is of course desired to implement the measurement in the shortest time possible. While the experimental procedure that was used included measurements in multiple PT excitation frequencies, it was found that reasonable estimation of the velocity can be obtained from measurements at two frequencies only. There are three main factors that will affect the choice of the PT excitation frequencies in realistic scenarios: the measurement duration, the sideband amplitude and undesired spectral broadening of the PA peak at the carrier frequency. The latter effect is expected to be generated by PA-induced temperature variations, flow variations in the VUT or patient movement during the measurement. While in the present experiment the broadening of the carrier was not significant, it is expected to be more pronounced in in-vivo measurements and to limit the minimum useable PT frequency.

Our study was focused on the case of a single vessel. The same approach, however, can be applied to a tissue containing a large number of blood vessels with intricate structure. In this case the underlying model and equation can be replaced with a bio-heat model [3,4]. This case is more similar to conventional TDF techniques, which generally measure a global parameter such as the tissue perfusion. Hence, it is to be expected that the proposed measurement approach could also be implemented over a volume containing multiple vessels and make use of the excellent analytical tools that were developed to make TDF a viable quantitative measurement method of blood flow [2–4].

Another aspect that must be taken into consideration is the safety of the method. Here as well, the similarity to TDF can be used for assessment of the possibility for safe implementation. TDF implementations which heat the tissue for flow measurement are already commercially available [25]. Accordingly, it is reasonable to assume that the proposed method can be implemented under similar restrictions. Experimentally we observed that it was possible to obtain good SNR with heating of less than 2C°, which is well below the threshold for thermal damage [24], using sinusoidal modulation of both diodes. Although an additional source for the PA excitation is required, as mentioned in Section 4, the CW modulation of the PA source is not essential to the method and can be replaced by a pulsed source that produces considerably less heating. The pulse repetition rate (PRF) of such PA source needs to be high enough compared to the thermal clearance rates (i.e. >100 Hz). Spectrally, the PA signal will show discrete harmonics separated by the PRF while the PT modulation sidebands will be in the near vicinity of each PA harmonic. Further optimization of the system, including focusing of the optical beam and/or ultrasonic transducer, may also improve the sensitivity and allow detection at lower optical powers.

6. Conclusions

In this paper a new method for implementing Thermal Diffusion Flowmetry was proposed, analyzed and tested experimentally. The method is based on photothermal excitation of a VUT and photoacoustic measurement of the resulting temperature variations. The PA-TDF method enabled measurement of blood velocity in a phantom vessel in a velocity range of 0-22 mm/s with errors smaller than 0.3 mm/s for all velocities. The relation between the velocity and the measured PA signal was derived from the heat and the PA equations via the use of a simplified lumped model. As the method is based on light absorption its multi-spectral implementation would provide significant spectroscopic information. In this way, for example, it is expected to accomplish a simultaneous measurement of both blood velocity and its oxygenation level (sO₂).