## 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.^{1}

Photoacoustic (PA) tomography (PAT) can quantify hemoglobin concentrations *in vivo* based on endogenous contrast with both fine spatial resolution and high sensitivity.^{2} 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\text{\hspace{0.17em}\hspace{0.17em}}{\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.^{3} In the optical ballistic regime, the lateral resolution of PAT is determined by optical focusing,^{4} 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,^{5} 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.^{6} Moreover, PAT inherently exploits optical absorption contrast, and thus it has a much higher sensitivity to optical absorption than other optical microscopy modalities^{7} 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^{3} 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 $:

## Eq. (1)

$${\mu}_{a}(\overline{r},\lambda )={\u03f5}_{\mathrm{HbR}}(\lambda )[\mathrm{HbR}](\overline{r})+{\u03f5}_{{\mathrm{HbO}}_{2}}(\lambda )[{\mathrm{HbO}}_{2}](\overline{r}),$$^{2}PAT images, however, are spatial mappings of the absorbed optical energy density $A(\overline{r},\lambda )={\mu}_{a}(\overline{r},\lambda )F(\overline{r},\lambda )$, where $F(\overline{r},\lambda )$ represents the local optical fluence in units of $J\text{\hspace{0.17em}\hspace{0.17em}}{\mathrm{m}}^{-2}$. To obtain the intrinsic quantity ${\mu}_{a}(\overline{r},\lambda )$ from $A(\overline{r},\lambda )$, we need to quantify $F(\overline{r},\lambda )$, which is usually unknown due to light transport in the scattering tissue.

$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^{19}

We applied this idea for the acoustic-resolution PA microscopy (AR-PAM).^{6} 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\text{\hspace{0.17em}\hspace{0.17em}}\mathrm{\mu m}$ for the 50 MHz transducer.

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\text{\hspace{0.17em}\hspace{0.17em}}\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:^{20}

## Eq. (3)

$$F(z,\lambda )={F}_{I}(\lambda ){\int}_{\mathrm{\Omega}}\mathrm{exp}[-{\mu}_{a}(\lambda )r]\frac{d\mathrm{\Omega}}{4\pi},$$^{21}(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}\mathrm{exp}[-{\mu}_{a}(\lambda )z/u]du,\text{\hspace{0.17em}}$$^{22}(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}\mathrm{exp}\left[\frac{-{\mu}_{a}(\lambda )ct}{u}\right]du\le ,$$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}}}.$$## Eq. (7)

$$\sum _{i=1}^{n}\int {|{\overline{S}}_{\mathrm{exp}}(\omega ,{\lambda}_{i})-{\overline{S}}_{\mathrm{model}}(\omega ,{\lambda}_{i})|}^{2}d\omega .$$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})={\u03f5}_{{\mathrm{HbO}}_{2}}({\lambda}_{i})[{\mathrm{HbO}}_{2}]+{\u03f5}_{\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.

## 2.2.

### Experimental Setups

In a phantom study, we used AR-PAM to image a fully oxygenated bovine blood phantom with four wavelengths (560, 565, 570, and 575 nm), where the absorption coefficients monochromatically increase [${\mu}_{a}(560)<{\mu}_{a}(565)<{\mu}_{a}(570)<{\mu}_{a}(575)$]. The blood phantom was in a 1 mm diameter tube, which was located 1.5-mm deep in the optical scattering medium (10% gelatin, 1% intralipid, and 5% ${\mathrm{CuCl}}_{2}$) and parallel to its surface [Fig. 1(b)].

In an *in vivo* experiment, we imaged an 8-mm by 8-mm region of the back of a nude mouse with two optical wavelengths (571 nm and 564 nm). Next, an optical phantom layer ($\sim 1.5\text{\hspace{0.17em}\hspace{0.17em}}\mathrm{mm}$ of 10% gelatin, 1% intralipid, and 5% ${\mathrm{CuCl}}_{2}$) was used to cover the back of the nude mouse, and the same region was imaged again.

## 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.

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\text{\hspace{0.17em}\hspace{0.17em}}\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\text{\hspace{0.17em}\hspace{0.17em}}\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.

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\text{\hspace{0.17em}\hspace{0.17em}}\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.

## Acknowledgments

We would like to thank Yan Liu, Dr. Xiao Xu, and Dr. Konstantin Maslov for useful discussions and Prof. James Ballard for editing the manuscript. This work was sponsored in part by the National Institutes of Health grants R01 EB008085, R01 EB000712, R01 CA134539, R01 CA157277, and U54 CA 136398. L. W. has a financial interest in Microphotoacoustics, Inc. and Endra, Inc., which, however, did not support this work.

## Appendices

## Appendix:

### Validation of Eq. (4) with Monte Carlo Simulations

The analytical solution of the fluence decay in blood [Eq. (4)] was validated by Monte Carlo simulations. The simulated object was assumed to be composed of two layers, while the top layer is the background tissue and the bottom layer is the blood. The optical properties of the background tissue are set to be ${\mu}_{a}=0.1\text{\hspace{0.17em}\hspace{0.17em}}{\mathrm{cm}}^{-1}$, ${\mu}_{s}=100\text{\hspace{0.17em}\hspace{0.17em}}{\mathrm{cm}}^{-1}$, and $g=0.9$, while those of the blood are set to be ${\mu}_{a}=200\text{\hspace{0.17em}\hspace{0.17em}}{\mathrm{cm}}^{-1}$, ${\mu}_{s}=300\text{\hspace{0.17em}\hspace{0.17em}}{\mathrm{cm}}^{-1}$, and $g=0.995$.^{21} Here ${\mu}_{a}$ is the absorption coefficient, ${\mu}_{s}$ is the scattering coefficient, and $g$ is the anisotropy factor. The thickness of the background tissue is 0.2 cm, which is greater than the transport mean free path (${l}_{t}^{\prime}=0.1\text{\hspace{0.17em}\hspace{0.17em}}\mathrm{cm}$). Therefore, the photons are almost completely diffused when they reach the blood layer.

In Fig. 4, the fluence decay profile in blood is plotted on a log scale, and it matches the analytical solution given by Eq. (4).