1.

Introduction

Oxygenated and deoxygenated hemoglobin concentrations ([${\mathrm{HbO}}_2$] and [HbR]) as well as the derived quantities—total hemoglobin concentration ([HbT]) and hemoglobin oxygen saturation (${\mathrm{sO}}_2=[{\mathrm{HbO}}_2]/[\mathrm{HbT}]$)—are important physiological parameters. Several techniques have been developed to quantify hemoglobin concentrations in vivo; however, they all have disadvantages. Diffuse optical tomography (DOT) has poor spatial resolution due to strong scattering in biological tissue; blood oxygen level dependent (BOLD) contrast magnetic resonance imaging (MRI) can monitor only [HbR]. Positron emission tomography (PET) measures only oxygen metabolism, and they not only suffer from poor spatial resolution but also require intravenous administration of radioactive isotopes.¹

Photoacoustic (PA) tomography (PAT) can quantify hemoglobin concentrations in vivo based on endogenous contrast with both fine spatial resolution and high sensitivity.² In PAT, the sample is typically illuminated by a pulsed laser. Following the absorption of optical energy, an initial pressure is generated via thermo-elastic expansion. The initial pressure then propagates as ultrasonic waves, which are detected by ultrasonic sensors. The strength of the initial pressure $P_0(\overline{r})$ in the unit of Pa at the location $\overline{r}$ in the biological tissue is proportional to the local absorbed optical energy density $A(\overline{r})$ in units of $J{\mathrm{m}}^{-3}$. From multiwavelength PA measurements, we can obtain the optical spectrum $A(\overline{r},\lambda )$ [i.e., $A(\overline{r})$ versus optical wavelength $\lambda$], which can be used to quantify hemoglobin concentrations in the same way as DOT.³ In the optical ballistic regime, the lateral resolution of PAT is determined by optical focusing,⁴ and thus it is comparable to that of other optical microscopy modalities. In the optical diffusive regime, however, the resolution of PAT is determined by ultrasonic waves,⁵ and PAT provides much better spatial resolution than DOT, in which the inverse algorithm is ill-posed. While DOT can only monitor ${\mathrm{sO}}_2$ which is volume-averaged over multiple blood vessels, PAT can pinpoint blood vessels and evaluate their individual ${\mathrm{sO}}_2$ levels.⁶ Moreover, PAT inherently exploits optical absorption contrast, and thus it has a much higher sensitivity to optical absorption than other optical microscopy modalities⁷ and DOT.

Nevertheless, using PAT to quantify hemoglobin concentrations conventionally requires knowledge of the local optical fluence.^8,9 In the quantification model, hemoglobin concentrations are derived³ from the optical absorption coefficient ${\mu }_a(\overline{r},\lambda )$ in the unit of ${\mathrm{m}}^{-1}$ by solving the following equation for multiple values of $\lambda$:

$F(\overline{r},\lambda )$ can be quantified in vivo either invasively^1,10 or noninvasively.^11–14 In the invasive method, an exogenous optical absorber or contrast agents with known absorption spectrum ${\mu }_a({\overline{r}}^{\prime },\lambda )$ is inserted at ${\overline{r}}^{\prime }$ near the target objects of interest. The PA amplitudes due to $A(\overline{r},\lambda )$ of the target objects of interest are normalized by the PA amplitude due to $A({\overline{r}}^{\prime },\lambda )$ of the calibration target. By assuming $F({\overline{r}}^{\prime },\lambda )\approx F(\overline{r},\lambda )$, we have ${\mu }_a(\overline{r},\lambda )={\mu }_a({\overline{r}}^{\prime },\lambda )A(\overline{r},\lambda )/A({\overline{r}}^{\prime },\lambda )$. This method physically compensates for the fluence attenuation; however, it is invasive. The noninvasive methods model light transport and acoustic wave propagation in biological tissues with the optical diffusion equation and the PA wave equation, respectively. The optical absorption coefficients are recovered by solving both equations. However, the inverse problems are usually ill-posed and computationally intensive.

The temporal profile of the PA signal has also been used to quantify optical absorption coefficients with reflection-mode PA imaging systems.^15–18 For example, if a pencil beam (with an infinitesimal beam width) incident perpendicularly on the blood vessel surface, then the energy deposition in the vessel decays exponentially along the beam propagation direction. If at wavelength $\lambda$ the reduced scattering coefficient ${\mu }_s^{\prime }(\lambda )$ is much smaller than the absorption coefficient ${\mu }_a(\lambda )$ of blood, fitting the received PA signals with Beer’s law yields ${\mu }_a(\lambda )$. Here, knowledge of the local optical fluence is not required, because ${\mu }_a(\lambda )$ is quantified from the relative temporal decay profile. However, various acoustic effects may distort the received PA signals. Acoustic attenuation in biological tissue has a power law dependence on the frequency, and therefore the shapes of the acoustic pulses change as they propagate. Also, as ultrasonic detectors have limited bandwidths, the detected PA signal is the convolution of the received acoustic pulse and the mechanical-electrical impulse response of the detector. As such, the temporal profiles no longer follow the exponential decay. Consequently, directly fitting the temporal profile for optical absorption coefficients can be inaccurate.

We propose a method to quantify optical absorption coefficients from the acoustic spectra of PA signals at multiple optical wavelengths for two main reasons. First, similar to fitting the PA temporal profiles for ${\mu }_a(\lambda )$, quantification of ${\mu }_a(\lambda )$ from acoustic spectra does not require the knowledge of the local optical fluence. Second, the response function of the ultrasonic detection system and the acoustic properties of the tissue are independent of the optical wavelength; therefore, their effects can be eliminated mathematically in the inverse problem of quantifying optical properties.

2.

Method

2.1.

Quantification of Optical Absorption Coefficients from the Acoustic Spectra of Multiwavelength PA Signals

Mathematically, the acoustic spectrum $S(\omega ,\lambda )$ of the received PA signal from an absorbing target object of interest can be written, based on the assumptions of linearity and shift-invariance, as¹⁹

Eq. (2)

$$S(\omega ,\lambda )=F_I(\lambda )O(\omega ,\lambda )a(\omega )H(\omega ).$$

Here, $\lambda$ is the optical wavelength. If the acoustic spectra are measured at multiple optical wavelengths, we have $\lambda ={\lambda }_i$ ($i=1,2,\cdots n$), where $n$ is the total number of optical wavelengths. $\omega$ is the acoustic angular frequency, $F_I(\lambda )$ is the optical fluence incident on the target object, $a(\omega )$ is the acoustic attenuation in the tissue between the target object and the ultrasonic transducer, and $H(\omega )$ is the system transfer function of the ultrasonic detection system. Since $a(\omega )$ and $H(\omega )$ are independent of the optical wavelength $\lambda$, and $F_I(\lambda )$ is independent of the acoustic frequency $\omega$, these factors can be eliminated, as will be shown later. $O(\omega ,\lambda )$ represents the original acoustic spectrum of the target object of interest (e.g., a blood vessel), and it equals to the Fourier transform of $O(t,\lambda )$, the time-domain PA signal from the target object of interest measured locally with an ideal ultrasonic detection system of infinite bandwidth at unit optical fluence. Here, $t$ denotes time. $O(t,\lambda )$ is related to the target object’s shape, size, optical properties, and the incident fluence distribution. If $O(\omega ,\lambda )$ is recovered, $O(t,\lambda )$ becomes known, and subsequently ${\mu }_a(\lambda )$ can be quantified through an appropriate model of $O(t,\lambda )$.

We applied this idea for the acoustic-resolution PA microscopy (AR-PAM).⁶ Figure 1(a) shows the AR-PAM system. A dye laser pumped by a Nd:YLF laser served as the irradiation source. The laser beam from the dye laser was delivered through an optical fiber and passes through a conical lens to provide a ring-shaped area of illumination. At each location, a focused ultrasonic transducer was employed to record the PA wave, which was converted into a one-dimensional (1-D) depth-resolved image (A-scan or A-line). A three-dimensional (3-D) image was achieved by raster scanning in the $x-y$ plane. The lateral resolution of the AR-PAM, determined by the focal diameter of the ultrasonic transducer, was $\sim 45\mathrm{\mu m}$ for the 50 MHz transducer.

Fig. 1

Schematic of the AR-PAM system and phantom experiment setups. (a) The schematic of the AR-PAM system. A dye laser pumped by a Nd:YLF laser was the irradiation source. The laser beam from the dye laser was delivered through an optical fiber and passes through a conical lens to provide a ring-shaped area of illumination. A focused ultrasonic transducer was employed to detect PA waves. (b) Use AR-PAM to image a blood-filled tube inserted 1.5 mm deep in the optical scattering medium. The fluence decay profile in the blood is modeled.

With AR-PAM, however, several assumptions can simplify the modeling of $O(t,\lambda )$. First, the surface of blood vessels with sufficiently large diameters (e.g., greater than 300 μm for the 45-μm lateral resolution) may be treated approximately as a flat surface locally. Second, when imaging blood vessels at depths greater than one transport mean free path ($l_t^{\prime }\sim 1\mathrm{mm}$ in biological tissue), we can assume that the light is completely diffused and can be considered as isotropic point sources [Fig. 1(b)]. The fluence in the blood vessel can be expressed as:²⁰

Eq. (3)

$$F(z,\lambda )=F_I(\lambda ){\int }_{\mathrm{\Omega }}\exp [-{\mu }_a(\lambda )r]\frac{d\mathrm{\Omega }}{4\pi },$$

where $z$ is the depth of the target layer in the blood vessel and $r$ is the distance between the isotropic point source and one point in the layer at the depth $z$. $d\mathrm{\Omega }=\sin \theta d\theta d\phi$ is the differential solid angle in the spherical coordinates, and the integration ranges of the polar angle $\theta$ and the azimuthal angle $\phi$ are [0, $\pi$] and [0, $2\pi$], respectively. $F_I(\lambda )$ is the incident fluence on the blood vessel, and $F_I(\lambda )\exp [-{\mu }_a(\lambda )r]/4\pi$ is the fluence at one point in the layer at the depth $z$ per unit solid angle. In blood, ${\mu }_s^{\prime }(\lambda )$ is much less than ${\mu }_a(\lambda )$ [e.g., $1.5{\mathrm{cm}}^{-1}\ll 241.0{\mathrm{cm}}^{-1}$ at 571 nm (21)], because the anisotropy factor ($g$) is so close to 1 in the optical spectral region we used²¹ (around 585 nm); therefore scattering is negligible.

Equation (3) can be further simplified to

Eq. (4)

$$F(z,\lambda )=F_I(\lambda ){\int }_0^1\exp [-{\mu }_a(\lambda )z/u]du,$$

where $u$ is a dimensionless scaling factor defined by $u=z/r$. Equation (4) was validated by Monte Carlo simulations.²² (see Appendix). The original PA signal from the target object induced by unit incident optical fluence is expressed as

Eq. (5)

$$O(t,\lambda )=\frac{\mathrm{\Gamma }{\mu }_a(\lambda )F(z,\lambda )}{F_I(\lambda )}=\mathrm{\Gamma }{\mu }_a(\lambda ){\int }_0^1\exp \left[\frac{-{\mu }_a(\lambda )ct}{u}\right]du\le ,$$

where $z$ is converted to time $t$ through the speed of sound $c:t=z/c$.

An algorithm, named the absorption-fitting algorithm, was developed to fit for the optical absorption coefficients as well as the relative incident fluence from the acoustic spectra at multiple wavelengths. We normalize the acoustic spectrum $S(\omega ,{\lambda }_i)$ ($i=1,2,\cdots n$) at each optical wavelength by the root-mean-square spectrum $S_{\mathrm{rms}}(\omega )=\sqrt{\sum _{i=1}^n{|S(\omega ,{\lambda }_i)|}^2/n}$ as

Eq. (6)

$$\overline{S}(\omega ,{\lambda }_i)=\frac{S(\omega ,{\lambda }_i)}{S_{\mathrm{rms}}(\omega )}=\frac{F_I({\lambda }_i)O(\omega ,{\lambda }_i)}{\sqrt{\frac{\sum _{i=1}^n{|F_I({\lambda }_i)O(\omega ,{\lambda }_i)|}^2}{n}}}.$$

The effects of both the system-dependent response $H(\omega )$ and the tissue related acoustic attenuation $a(\omega )$ are cancelled in the normalization process. The optical absorption coefficient and the relative incident fluence at each optical wavelengths ${\lambda }_i$ ($i=1,2,\cdots n$) are optimized by minimizing the following cost function:

Eq. (7)

$$\sum _{i=1}^n\int {|{\overline{S}}_{\exp }(\omega ,{\lambda }_i)-{\overline{S}}_{\mathrm{model}}(\omega ,{\lambda }_i)|}^{2}d\omega .$$

Here, ${\overline{S}}_{\exp }(\omega ,{\lambda }_i)$ is calculated from the experimentally measured acoustic spectra using Eq. (6); ${\overline{S}}_{\mathrm{model}}(\omega ,{\lambda }_i)$ is computed from our theoretical model containing the fitting parameters ${\mu }_a({\lambda }_i)$ and $F_I({\lambda }_i)$ using Eq. (5), Fourier transformation, and Eq. (6). The Fourier transformation, converting $O(t,\lambda )$ to $O(\omega ,\lambda )$, was implemented numerically. The minimization process is accomplished with a nonlinear optimization algorithm—the Levenberg—Marquardt algorithm (LMA)—available as a standard optimization package in Matlab.

In special circumstances, when the compositions of the absorbing medium are known, we can directly fit the acoustic spectra for the concentrations of the optical absorbers. For example, in blood the major optical absorbers are ${\mathrm{HbO}}_2$, HbR, and water. In the optical spectral region we used (around 585 nm), the absorption of water is negligible in comparison to that of ${\mathrm{HbO}}_2$ and HbR. An algorithm, named the concentration-fitting, was developed to directly fit for the hemoglobin concentrations ([${\mathrm{HbO}}_2$]and [HbR]) in blood. In the cost function [Eq. (7)], we substitute the optical absorption coefficient at ${\lambda }_i$ by ${\mu }_a({\lambda }_i)={ϵ}_{{\mathrm{HbO}}_2}({\lambda }_i)[{\mathrm{HbO}}_2]+{ϵ}_{\mathrm{HbR}}({\lambda }_i)[\mathrm{HbR}]$. Then [${\mathrm{HbO}}_2$], [HbR], and $F_I({\lambda }_i)$ are optimized from Eq. (7) by the LMA.

Compared with the concentration-fitting algorithm, which requires the knowledge of the compositions of the absorbing medium, the absorption-fitting algorithm is universally applicable. However, when the number of optical wavelengths measured is larger than the types of absorbers in the absorbing medium, the concentration-fitting algorithm, which involves a reduced number of fitting parameters, is preferred.

3.

Results and Discussion

In the phantom study, PA A-line signals at four wavelengths were acquired from the samples [Fig. 2(a)] and were normalized to the peak of the PA signal measured at 570 nm. The corresponding acoustic spectra of the PA signals [$S(\omega ,\lambda )$] were calculated [Fig. 2(b)]. Light penetrated deeper at longer wavelengths, and the corresponding PA signals decayed more slowly in the time domain. As the wavelengths increased, the fluence decayed faster in the tube, and the produced PA signals were sharper with time. Therefore, the acoustic spectra acquired at longer optical wavelengths contained more high-frequency components.

Fig. 2

Quantification of the optical absorption coefficients and hemoglobin concentrations of oxygenated bovine blood phantom. (a) PA time domain A-scan signals of the oxygenated-bovine-blood-filled tube phantom buried in the scattering medium. (b) Acoustic spectra of the PA A-scan signals. (c) Comparison of the optical absorption coefficients quantified with the acoustic spectral method and measured by spectrophotometry. (d) Relative incident fluence quantified with the acoustic spectral method. (e) Comparison of the relative optical absorption coefficients quantified from the PA signal amplitudes with and without the fluence compensation and measured by spectrophotometry. (f) Comparison of the [${\mathrm{HbO}}_2$] and [HbR] quantified with the acoustic spectral method and spectrophotometry. Amp.: the amplitude method with and without the fluence compensation. Abs.-fitting: the absorption-fitting algorithm. Conc.-fitting: the concentration-fitting algorithm.

By fitting the acoustic spectra with the absorption-fitting algorithm, we found the absorption coefficients at four wavelengths [Fig. 2(c)] as well as the relative incident fluence [Fig. 2(d)]. In comparison to the gold standard—spectrophotometry, the absorption-fitting algorithm quantified the absorption coefficients with relative errors of 1.9%, 2.8%, 1.5%, and 1.3% at the optical wavelengths of 560, 565, 570, and 575 nm, respectively [Fig. 2(c)], where the error bars indicate the fitting standard errors. In Fig. 2(d), the relative incident fluence at each of the four wavelengths is normalized to the maximum incident fluence at 570 nm. Normalizing the PA signal amplitude at each wavelength by the corresponding relative incident fluence yields the relative optical absorption coefficient [Fig. 2(e)]. This normalization process is referred to as fluence compensation, because it calibrates the wavelength-dependent fluence attenuation. In contrast, quantifying the absorption coefficients from the PA signal amplitudes without the fluence compensation is inaccurate [Fig. 2(e)]. From the optical absorption coefficients at four wavelengths, we obtained the [${\mathrm{HbO}}_2$] and [HbR] based on Eq. (1) [Fig. 2(f)]. [${\mathrm{HbO}}_2$] and [HbR] [Fig. 2(f)] as well as the relative incident fluence [Fig. 2(d)] were also fitted from PA acoustic spectra using the concentration-fitting algorithm, and the optical absorption coefficients at the four wavelengths were obtained based on Eq. (1) [Fig. 2(c)]. In comparison to the spectrophotometry, the absorption coefficients were obtained with relative errors of 1.6%, 2.2%, 1.3%, and 1.1% at the optical wavelengths of 560, 565, 570, and 575 nm, respectively. The two methods achieved similar quantification accuracy in terms of error relative to the gold standard.

Figure 3(a) and 3(b) shows the PA maximum amplitude projection (MAP) images, in which each point corresponds to the maximum value of a Hilbert transformed PA A-scan. Since these blood vessels are shallow [$\sim 150\mathrm{\mu m}$ deep, as shown in Fig. 3(c)], the effect of wavelength-dependent fluence attenuation is negligible. We obtained the ${\mathrm{sO}}_2$ and [HbT] control images [Fig. 3(d) and 3(e)] based on Eq. (1) from the PA signal amplitudes. From the ${\mathrm{sO}}_2$ image, the arteries and the veins can be clearly identified [Fig. 3(d)]. Figure 3(f) and 3(g) shows the PA MAP images of the same region with the optical phantom layer, where the blood vessels are $\sim 1.6\mathrm{mm}$ deep [Fig. 3(h)]. In this case, if we ignored the wavelength-dependent fluence variations induced by the optical phantom layer, the ${\mathrm{sO}}_2$ quantified from the amplitudes of the PA signals [Fig. 3(i)] became inaccurate by approximately 32%. The average quantification error of the relative [HbT] [Fig. 3(j)], however, was only 5%, because the blood vessels in the FOV was located approximately at the same depth [Fig. 3(h)]; thus, the fluence attenuation effect was eliminated by the normalization process. Nevertheless, the amplitude method provides only relative quantifications.

Fig. 3

Quantification of the ${\mathrm{sO}}_2$ and [HbT] of blood vessels in the back of a nude mouse in vivo. MAP image acquired at (a) 571 nm and (b) 564 nm. (c) Depth-encoded image of blood vessels. (d) Control images of ${\mathrm{sO}}_2$ and (e) [HbT]. MAP image acquired with an optical phantom layer at (f) 571 nm and (g) 564 nm. (h) Depth-encoded image of blood vessels below the optical phantom layer. (i) ${\mathrm{sO}}_2$ and (j) [HbT] quantified from the amplitude method. Absolute optical absorption coefficients at (k) 571 nm and (l) 564 nm and (m) relative incident fluence quantified with the absorption-fitting algorithm. (n) ${\mathrm{sO}}_2$ and (o) [HbT] quantified from the absolute optical absorption coefficients. (p) [${\mathrm{HbO}}_2$], (q) [HbR], and (r) relative incident fluence quantified with the concentration-fitting algorithm. (s) ${\mathrm{sO}}_2$ and (t) [HbT] quantified from [${\mathrm{HbO}}_2$] and [HbR]. (u) ${\mathrm{sO}}_2$ and (v) [HbT] quantified from amplitude method with the fluence compensation. Ctrl.: the control image. Amp.: the amplitude method with and without the fluence compensation. A.S.: the acoustic spectral method. Abs.-fitting: the absorption-fitting algorithm. Conc.-fitting: the concentration-fitting algorithm.

Since the light was completely diffused by the optical phantom layer and the blood vessel diameters were much larger than the imaging resolution, the assumptions in our acoustic spectral method were valid. For each scanning position, the optical absorption coefficients [Fig. 3(k) and 3(l)] as well as the relative incident fluence [Fig. 3(m)] were acquired with the absorption-fitting algorithm. Then the ${\mathrm{sO}}_2$ and [HbT] images [Fig. 3(n) and 3(o)] were quantified from the absolute optical absorption coefficients. By comparing the quantification results to the control images, we found the absorption-fitting algorithm achieved average errors of 8% and 6% in quantifying ${\mathrm{sO}}_2$ and [HbT], respectively. Here, we normalized [HbT] by the maximum to make a fair comparison. We also acquired the [${\mathrm{HbO}}_2$] [Fig. 3(p)], the [HbR] [Fig. 3(q)], and the relative incident fluence [Fig. 3(r)] using the concentration-fitting algorithm. Then the ${\mathrm{sO}}_2$ and [HbT] images [Fig. 3(s) and 3(t)] are derived from [${\mathrm{HbO}}_2$] and [HbR] images. The concentration-fitting algorithm achieved average errors of 7% and 6% in quantifying ${\mathrm{sO}}_2$ and [HbT], respectively.

Normalizing the amplitudes of the MAP PA images at two wavelengths by the corresponding relative incident fluences yields the relative optical absorption coefficients. Specifically, the PA MAP image at 571 nm [Fig. 3(f)] is divided by 1, while the PA MAP image at 564 nm [Fig. 3(g)] is divided, point by point, by the relative incident fluence [Fig. 3(m)]. After the fluence compensation, the ${\mathrm{sO}}_2$ [Fig. 3(u)] and the relative [HbT] [Fig. 3(v)] were quantified from the calibrated PA amplitudes with average errors of 9% and 5%, respectively.

We should note that the acoustic spectral method also has limitations. First, the object geometrical effect has been neglected. The current quantification model treats the surface of blood vessels with sufficiently large diameters approximately as a flat surface locally. Therefore, it cannot be applied to vessels whose diameters are comparable or smaller than the lateral resolution. Second, the incident light on the vessel surface was assumed to be completely diffused. Thus, the quantification model can be applied only to the blood vessels at depths greater than one transport mean free path ($l_t^{\prime }\sim 1\mathrm{mm}$ in biological tissue). Third, the ultrasonic transducer impulse response was assumed to be spatially invariant, which requires the blood vessels to be located within the acoustic focal zone.

4.

Conclusions

In summary, we have proposed a method to quantify [HbT] and ${\mathrm{sO}}_2$ in vivo using acoustic spectra of PA signals from multiple optical wavelength measurements. Using AR-PAM, we first demonstrated this method with phantom experiments and then quantified the [HbT] and ${\mathrm{sO}}_2$ in a live mouse. The acoustic spectral method provides greater quantification accuracy than the conventional amplitude method in the optical diffusive regime.