## 1.

## Introduction

Accurate estimation of blood oxygenation is of great importance in making critical decisions associated with clinical situations. Blood oxygen level can be monitored noninvasively by using a near-infrared spectroscopy method. It is a very convenient method and uses translucent body parts (e.g., earlobes, fingertips, etc.) to determine the blood oxygen saturation.^{1} However, the major limitations of this technique are the (*i*) small penetration depth due to strong light scattering and (*ii*) the inability to distinguish between venous and arterial blood due to poor spatial resolution. The optoacoustic (OA) technique promises to provide a real-time, noninvasive blood oxygenation monitoring method by overcoming these drawbacks and is expected to be useful for the assessment of oxygenation in many sites of physiological interest, including the brain.^{2} Note that the brain is very sensitively dependent on oxygen for its normal function, and accurate quantification of the brain tissue oxygenation can play a crucial role for treating patients with traumatic brain injury or neurophysiological disorders.^{1}

In OA methods, a short pulsed laser is incident on a tissue medium, which absorbs the light energy and subsequently emits pressure waves due to the thermoelastic expansion. The pressure pulse is detected using an ultrasonic transducer, and the signal can be utilized to image, as well as to characterize, the medium.^{2} Using this modality, quantification of oxygenation of deep tissue regions can be made. This is because oxygen estimates are obtained by analyzing the detected pressure waves, whose scattering in tissue is two to three orders of magnitude weaker than optical scattering and thus can travel a longer distance due to lower attenuation.^{2} The spatial resolution of the measurement can be improved by choosing an appropriate ultrasonic transducer and one that facilitates to obtain localized information. OA tomography and microscopy instruments have been designed for small animal studies by exploiting these advantages.^{2, 3} Moreover, in recent years several novel research tools using OA techniques have been developed to assist in projects such as cancer cell visualization using gold nanoparticles,^{4} circulating melanoma cell detection,^{5, 6} port wine stain depth determination,^{7} single-cell imaging,^{8} and vascular imaging,^{9} etc.

A number of experimental studies have also been performed to investigate the potential of the OA technique to monitor blood oxygenation. For example, Wang
^{10} imaged small animal cerebral hemoglobin concentration and hemoglobin oxygenation noninvasively using a multiwavelength laser-based OA system. Esenaliev
^{11} carried out an *in vitro* experimental work and showed that OA signal amplitude increased monotonically with blood oxygenation when probed using a 1064-nm laser source. The same group also measured *in vivo* sheep blood oxygenation from the external jugular vein in the neck area employing three laser wavelengths.^{12} They demonstrated that the OA signal amplitude decreased for 700 nm and increased at 1064 nm when blood oxygenation increased approximately from 20 to 95%. However, to the best of our knowledge the effects of erythrocyte oxygenation on OA signals have never been studied at the red blood cell (RBC) level.

A theoretical framework has recently been developed to express the OA field generated by an ensemble of RBCs, approximated as fluid spheres with uniform size and identical physical properties.^{13} The resultant pressure field generated by many optical absorbers in this framework has been obtained by using the linear superposition principle for the waves originating from those absorbers. The field emitted by an individual spherical absorber was evaluated by solving the wave equation for the pressure field (generated due to the absorption of light) in the frequency domain and using the appropriate boundary conditions.^{14} This theoretical model has also been utilized to investigate the capability of the OA technique to differentiate samples with nonaggregated and aggregated RBCs.^{13} In this study, it was shown that spectral intensity increased ∼11 dB at 15.6 MHz at an aggregation level compared to that of nonaggregated blood.

The objective of the paper is to discuss the generalized version of the theoretical model reported in Ref. 13. This framework enables one to compute the OA field generated by a collection of spherical absorbers with heterogeneous sizes and different physical properties. The presented model has been explored here to examine how the OA field would vary when the proportion of oxygenated and deoxygenated RBCs [(RBCOs) and (RBCDs), respectively] varies as a function of optical irradiation wavelength. The oxygen saturation level was assumed to be 100% for RBCOs and 0% for RBCDs. The light absorption coefficients of these cells are generally different (though similar at the isosbestic points), and consequently, OA emission strengths would be different. Therefore, the OA signals are expected to differ for samples with different proportions of RBCOs and RBCDs, even when probed with the same laser source. A Monte Carlo algorithm was implemented to simulate spatially random distributions of RBCOs and RBCDs. The OA signals were simulated for those configurations and subsequently analyzed to obtain signal envelope and spectral features. Our results show that the simulated trends are in good agreement with published experimental results, suggesting the potential of this model to explain experimental observations. Thus, this model could be useful to develop insight about the generation of OA signals from mixtures of RBCOs and RBCDs and also to predict signal behaviors at different blood oxygen saturation levels.

The organization of the paper is as follows. The generalized theoretical model is derived in Sec. 2. The simulation method is described in Sec. 3. The simulation results and discussion of the results are presented in Sec. 4 and 5, respectively. The conclusions of this study are highlighted in Sec. 6.

## 2.

## Theory

## 2.1.

### Optoacoustic Field Calculation for a Single Spherical Absorber

Consider a fluid sphere with density ρ_{s}, speed of sound *v*
_{s}, and radius *a*. This sphere is suspended in another fluid medium with properties ρ_{f} and *v*
_{f} (the density and speed of sound, respectively) and is heated by optical radiation. The subscript s is used to indicate the absorbing medium and f refers to the surrounding fluid medium. If *H* is the amount of heat deposited by the incident radiation per unit time per unit volume, the wave equation for the pressure field inside the absorber can be presented as,^{14}

## 1

[TeX:] \documentclass[12pt]{minimal}\begin{document} \begin{eqnarray} \nabla ^2p-\frac{1}{v_{\rm s}^2}\frac{\partial ^2 p}{\partial t^2}=-\frac{\beta }{C_P}\frac{\partial H}{\partial t}, \end{eqnarray} \end{document} $$\begin{array}{c}\hfill {\nabla}^{2}p-\frac{1}{{v}_{\mathrm{s}}^{2}}\frac{{\partial}^{2}p}{\partial {t}^{2}}=-\frac{\beta}{{C}_{P}}\frac{\partial H}{\partial t},\end{array}$$*C*

_{P}are the thermal expansion coefficient and isobaric heat capacity per unit mass, respectively, for the fluid medium inside the absorber. In deriving Eq. 1, it has been assumed that no heat conduction takes place before the acoustic pulse is launched. If the optical radiation with intensity

*I*

_{0}varies sinusoidally with time and propagates along the

*x*-axis, then the heating function can be written as,

*H*(

*x*,

*t*) = μ

*I*

_{0}

*e*

^{−iωt}. Here, μ is the absorption coefficient of the fluid medium within the absorber and ω is the modulation frequency of the incident optical radiation. For such a heating function, the time-independent wave equation for the pressure field is given by,

^{14}

## 2

[TeX:] \documentclass[12pt]{minimal}\begin{document} \begin{eqnarray} \nabla ^2p+k^2 p=\left\lbrace \begin{array}{l@{\quad}l}\frac{i\omega \mu \beta I_0}{C_P} & \rm {inside\; the\; absorber,}\\[6pt] 0 & \rm {outside\; the\; absorber}, \end{array} \right.\ \end{eqnarray} \end{document} $$\begin{array}{c}\hfill {\nabla}^{2}p+{k}^{2}p=\left\{\begin{array}{cc}\hfill \frac{i\omega \mu \beta {I}_{0}}{{C}_{P}}& \hfill \mathrm{inside}\phantom{\rule{0.28em}{0ex}}\mathrm{the}\phantom{\rule{0.28em}{0ex}}\mathrm{absorber},\\ \hfill 0& \hfill \mathrm{outside}\phantom{\rule{0.28em}{0ex}}\mathrm{the}\phantom{\rule{0.28em}{0ex}}\mathrm{absorber},\end{array}\right.\phantom{\rule{0.33em}{0ex}}\end{array}$$*k*is the wave number of the pressure wave. The analytical solution of Eq. 2 for a spherical absorber can readily be obtained by demanding the continuity of the pressure field and the normal component of particle velocity at the boundary. Thus, the expression of the pressure field after implementing these boundary conditions at a large distance

*r*from the absorber becomes,

^{14}

## 3

[TeX:] \documentclass[12pt]{minimal}\begin{document} \begin{eqnarray} &&\hspace*{-12pt}p_{\rm f}^{{\rm single}}(r, k_f)\nonumber\\ &&=\frac{i\mu \beta I_0 v_{\rm s}a^2}{C_P r}\frac{[\sin \hat{q}-\hat{q}\,\cos \hat{q}]e^{ik_{\rm f}(r-a)}}{\hat{q}^2[(1-\hat{\rho })(\sin \hat{q}/\hat{q})-\cos \hat{q}+i\hat{\rho }\hat{v}\,\sin \hat{q}]},\nonumber\\ \end{eqnarray} \end{document} $$\begin{array}{ccc}& & {p}_{\mathrm{f}}^{\mathrm{single}}(r,{k}_{f})\hfill \\ & & =\frac{i\mu \beta {I}_{0}{v}_{\mathrm{s}}{a}^{2}}{{C}_{P}r}\frac{[\mathrm{sin}\widehat{q}-\widehat{q}\phantom{\rule{0.16em}{0ex}}\mathrm{cos}\widehat{q}]{e}^{i{k}_{\mathrm{f}}(r-a)}}{{\widehat{q}}^{2}[(1-\widehat{\rho})(\mathrm{sin}\widehat{q}/\widehat{q})-\mathrm{cos}\widehat{q}+i\widehat{\rho}\widehat{v}\phantom{\rule{0.16em}{0ex}}\mathrm{sin}\widehat{q}]},\hfill \end{array}$$*k*

_{f}= ω/

*v*

_{f}. The superscript single states that only one OA source is considered. The corresponding time-dependent pressure field for a δ function heating pulse can be derived by using the Fourier transformation and that yields,

^{14}

## 4

[TeX:] \documentclass[12pt]{minimal}\begin{document} \begin{eqnarray} &&\hspace*{-7pt}p_{\rm f}^{\rm single}(r, t) = \frac{i\mu \beta F v_{\rm s} a^2}{2\pi C_P r}\nonumber\\ &&\hspace*{2pt}\times\int _{-\infty }^{\infty }\!\! d\omega \frac{[\sin \hat{q}{-}\hat{q}\cos \hat{q}]}{\hat{q}^2[(1{-}\hat{\rho })(\sin \hat{q}/\hat{q}){-}\cos \hat{q}{+}i\hat{\rho }\hat{v}\,\sin \hat{q}]}e^{ik_{\rm f}(r-a-v_{\rm f}t)},\nonumber\\ \end{eqnarray} \end{document} $$\begin{array}{ccc}& & {p}_{\mathrm{f}}^{\mathrm{single}}(r,t)=\frac{i\mu \beta F{v}_{\mathrm{s}}{a}^{2}}{2\pi {C}_{P}r}\hfill \\ & & \phantom{\rule{2.0pt}{0ex}}\times {\int}_{-\infty}^{\infty}d\omega \frac{[\mathrm{sin}\widehat{q}-\widehat{q}\mathrm{cos}\widehat{q}]}{{\widehat{q}}^{2}[(1-\widehat{\rho})(\mathrm{sin}\widehat{q}/\widehat{q})-\mathrm{cos}\widehat{q}+i\widehat{\rho}\widehat{v}\phantom{\rule{0.16em}{0ex}}\mathrm{sin}\widehat{q}]}{e}^{i{k}_{\mathrm{f}}(r-a-{v}_{\mathrm{f}}t)},\hfill \end{array}$$*F*is the optical radiation fluence.

## 2.2.

### Optoacoustic Field Calculation for an Ensemble of Spherical Absorbers

If the illuminated region contains a collection of absorbers, then the resultant pressure field can be obtained by using the linear superposition principle. This is essentially based on the single-particle approach and provides a simple way to model the generated pressure field. The single-particle model works well for a turbid medium and has been extensively used to explain experimental results in light and ultrasound scattering problems.^{15, 16} If the observation point is far away from the illuminated region, then the pressure field thus in this framework can be expressed as

## 5

[TeX:] \documentclass[12pt]{minimal}\begin{document} \begin{eqnarray} p_{\rm f}^{\rm {ensemble}}(r, k_{\rm f}) & = & \sum _{n=1}^{N} \frac{i\mu _n\beta _n I_0 v_{{\rm s}n} a_n^2}{C_{Pn}}\frac{[\sin \hat{q}_n-\hat{q}_n\,\cos \hat{q}_n]}{\hat{q}_n^2[(1-\hat{\rho }_n)(\sin \hat{q}_n/\hat{q}_n)-\cos \hat{q}_n+i\hat{\rho }_n\hat{v}_n\,\sin \hat{q}_n]} \frac{e^{i{\rm k}_{\rm f}(|{\bf {r}}-{\bf {r}_n}|-a_n)}}{|{\bf {r}}-{\bf {r}_n}|}\nonumber \\ & \approx & \sum _{n=1}^{N} \frac{i\mu _n\beta _n I_0 v_{{\rm s}n} a_n^2}{C_{Pn}r}\frac{[\sin \hat{q}_n-\hat{q}_n\cos \hat{q}_n]e^{ik_f(r-a_n)}}{\hat{q}_n^2[(1-\hat{\rho }_n)(\sin \hat{q}_n/\hat{q}_n)-\cos \hat{q}_n+i\hat{\rho }_n\hat{v}_n\,\sin \hat{q}_n]} e^{-i {\bf k}_{\rm f} \cdot {\bf {r}_n}}.\end{eqnarray} \end{document} $$\begin{array}{ccc}\hfill {p}_{\mathrm{f}}^{\mathrm{ensemble}}(r,{k}_{\mathrm{f}})& =& \sum _{n=1}^{N}\frac{i{\mu}_{n}{\beta}_{n}{I}_{0}{v}_{\mathrm{s}n}{a}_{n}^{2}}{{C}_{Pn}}\frac{[\mathrm{sin}{\widehat{q}}_{n}-{\widehat{q}}_{n}\phantom{\rule{0.16em}{0ex}}\mathrm{cos}{\widehat{q}}_{n}]}{{\widehat{q}}_{n}^{2}[(1-{\widehat{\rho}}_{n})(\mathrm{sin}{\widehat{q}}_{n}/{\widehat{q}}_{n})-\mathrm{cos}{\widehat{q}}_{n}+i{\widehat{\rho}}_{n}{\widehat{v}}_{n}\phantom{\rule{0.16em}{0ex}}\mathrm{sin}{\widehat{q}}_{n}]}\frac{{e}^{i{\mathrm{k}}_{\mathrm{f}}\left(\right|\mathbf{r}-{\mathbf{r}}_{\mathbf{n}}|-{a}_{n})}}{|\mathbf{r}-{\mathbf{r}}_{\mathbf{n}}|}\hfill \\ & \approx & \sum _{n=1}^{N}\frac{i{\mu}_{n}{\beta}_{n}{I}_{0}{v}_{\mathrm{s}n}{a}_{n}^{2}}{{C}_{Pn}r}\frac{[\mathrm{sin}{\widehat{q}}_{n}-{\widehat{q}}_{n}\mathrm{cos}{\widehat{q}}_{n}]{e}^{i{k}_{f}(r-{a}_{n})}}{{\widehat{q}}_{n}^{2}[(1-{\widehat{\rho}}_{n})(\mathrm{sin}{\widehat{q}}_{n}/{\widehat{q}}_{n})-\mathrm{cos}{\widehat{q}}_{n}+i{\widehat{\rho}}_{n}{\widehat{v}}_{n}\phantom{\rule{0.16em}{0ex}}\mathrm{sin}{\widehat{q}}_{n}]}{e}^{-i{\mathbf{k}}_{\mathrm{f}}\xb7{\mathbf{r}}_{\mathbf{n}}}.\hfill \end{array}$$*n*is used to indicate the n'th absorber with position vector

**r**

_{n}and radius

*a*

_{n}. Other physical properties of the n'th absorber are designated as μ

_{n}, β

_{n}, and

*C*

_{Pn}. The direction of the observation is defined by the wave vector

**k**

_{f}. Figure 1 presents the geometry of the OA setup. The illuminated region contains

*N*OA sources, and therefore, their contributions are summed up in Eq. 6. In this derivation, it has been assumed that the attenuation of the light beam is negligible during its propagation through the medium. Therefore, the incident-light intensity is the same for all OA sources and independent of their spatial locations. The interactions of generated acoustic waves with particles have also been ignored. The time-dependent pressure field for a δ function heating pulse can be obtained as

## 6

[TeX:] \documentclass[12pt]{minimal}\begin{document} \begin{eqnarray} p_{\rm f}^{\rm {ensemble}}(r, t) \approx \int _{-\infty }^{\infty } d\omega \left[\sum _{n=1}^{N}\frac{i\mu _n\beta _n F v_{{\rm s}n} a_n^2}{2\pi C_{Pn}r} \frac{[\sin \hat{q}_n-\hat{q}_n\,\cos \hat{q}_n]e^{ik_{\rm f}(r-a_n-v_{\rm f}t)}}{\hat{q}_n^2[(1-\hat{\rho }_n)(\sin \hat{q}_n/\hat{q}_n)-\cos \hat{q}_n+i\hat{\rho }_n\hat{v}_n\,\sin \hat{q}_n]} e^{-i {\bf k}_{\rm f} \cdot {\bf {r}_n}}\right]. \end{eqnarray} \end{document} $$\begin{array}{c}\hfill {p}_{\mathrm{f}}^{\mathrm{ensemble}}(r,t)\approx {\int}_{-\infty}^{\infty}d\omega \left[\sum _{n=1}^{N}\frac{i{\mu}_{n}{\beta}_{n}F{v}_{\mathrm{s}n}{a}_{n}^{2}}{2\pi {C}_{Pn}r}\frac{[\mathrm{sin}{\widehat{q}}_{n}-{\widehat{q}}_{n}\phantom{\rule{0.16em}{0ex}}\mathrm{cos}{\widehat{q}}_{n}]{e}^{i{k}_{\mathrm{f}}(r-{a}_{n}-{v}_{\mathrm{f}}t)}}{{\widehat{q}}_{n}^{2}[(1-{\widehat{\rho}}_{n})(\mathrm{sin}{\widehat{q}}_{n}/{\widehat{q}}_{n})-\mathrm{cos}{\widehat{q}}_{n}+i{\widehat{\rho}}_{n}{\widehat{v}}_{n}\phantom{\rule{0.16em}{0ex}}\mathrm{sin}{\widehat{q}}_{n}]}{e}^{-i{\mathbf{k}}_{\mathrm{f}}\xb7{\mathbf{r}}_{\mathbf{n}}}\right].\end{array}$$Equations 5, 6 can also be used to model the pressure field generated by a mixture of different populations of fluid particles. Thus, the equations can be employed for a mixture of RBCOs and RBCDs, if one approximates erythrocytes as fluid spheres. Moreover, it can be assumed that the size and biophysical properties of RBCOs and RBCDs are identical. With these approximations, Eq. 5 reduces to

## 7

[TeX:] \documentclass[12pt]{minimal}\begin{document} \begin{eqnarray} p_{\rm f}^{\rm {ensemble}}(r, k_{\rm f}) \approx \frac{i\beta I_0 v_{\rm s}a^2}{C_{P}r} \frac{ [\sin \hat{q}-\hat{q}\,\cos \hat{q}] e^{ik_{\rm f}(r-a)}}{\hat{q}^2[(1-\hat{\rho })(\sin \hat{q}/\hat{q})-\cos \hat{q}+i\hat{\rho }\hat{v}\,\sin \hat{q}]} \left[\mu _{\rm O} \sum _{n=1}^{N_{\rm O}} e^{-i {\bf k}_{\rm f} \cdot {\bf {r}_n}}+\mu _{\rm D} \sum _{n=N_{\rm O}+1}^{N} e^{-i {\bf k}_{\rm f} \cdot {\bf {r}_n}} \right], \end{eqnarray} \end{document} $$\begin{array}{c}\hfill {p}_{\mathrm{f}}^{\mathrm{ensemble}}(r,{k}_{\mathrm{f}})\approx \frac{i\beta {I}_{0}{v}_{\mathrm{s}}{a}^{2}}{{C}_{P}r}\frac{[\mathrm{sin}\widehat{q}-\widehat{q}\phantom{\rule{0.16em}{0ex}}\mathrm{cos}\widehat{q}]{e}^{i{k}_{\mathrm{f}}(r-a)}}{{\widehat{q}}^{2}[(1-\widehat{\rho})(\mathrm{sin}\widehat{q}/\widehat{q})-\mathrm{cos}\widehat{q}+i\widehat{\rho}\widehat{v}\phantom{\rule{0.16em}{0ex}}\mathrm{sin}\widehat{q}]}\left[{\mu}_{\mathrm{O}}\sum _{n=1}^{{N}_{\mathrm{O}}}{e}^{-i{\mathbf{k}}_{\mathrm{f}}\xb7{\mathbf{r}}_{\mathbf{n}}}+{\mu}_{\mathrm{D}}\sum _{n={N}_{\mathrm{O}}+1}^{N}{e}^{-i{\mathbf{k}}_{\mathrm{f}}\xb7{\mathbf{r}}_{\mathbf{n}}}\right],\end{array}$$_{O}and μ

_{D}are the absorption coefficients for the RBCOs and RBCDs, respectively. Here,

*N*

_{O}and

*N*

_{D}=

*N*−

*N*

_{O}are the numbers of RBCOs and RBCDs, respectively. Furthermore, for a mixture of such cells, Eq. 6 simplifies into

## 8

[TeX:] \documentclass[12pt]{minimal}\begin{document} \begin{eqnarray} &&\hspace*{-15pt}p_{\rm f}^{\rm {ensemble}}(r, t) \approx \frac{i\beta F v_{\rm s}a^2}{2\pi C_Pr}\nonumber\\ &&\times \int _{-\infty }^{\infty } d\omega \frac{[\sin \hat{q}-\hat{q}\,\cos \hat{q}]e^{ik_{\rm f}(r-a-v_{\rm f}t)}}{\hat{q}^2[(1-\hat{\rho })(\sin \hat{q}/\hat{q})-\cos \hat{q}+i\hat{\rho }\hat{v}\,\sin \hat{q}]}\nonumber\\ &&\times \left[\mu _{\rm O} \sum _{n=1}^{N_{\rm O}} e^{-i {\bf k}_{\rm f} \cdot {\bf {r}_n}}+\mu _{\rm D} \sum _{n=N_O+1}^{N} e^{-i {\bf k}_{\rm f} \cdot {\bf {r}_{\rm n}}} \right]. \end{eqnarray} \end{document} $$\begin{array}{ccc}& & {p}_{\mathrm{f}}^{\mathrm{ensemble}}(r,t)\approx \frac{i\beta F{v}_{\mathrm{s}}{a}^{2}}{2\pi {C}_{P}r}\hfill \\ & & \times {\int}_{-\infty}^{\infty}d\omega \frac{[\mathrm{sin}\widehat{q}-\widehat{q}\phantom{\rule{0.16em}{0ex}}\mathrm{cos}\widehat{q}]{e}^{i{k}_{\mathrm{f}}(r-a-{v}_{\mathrm{f}}t)}}{{\widehat{q}}^{2}[(1-\widehat{\rho})(\mathrm{sin}\widehat{q}/\widehat{q})-\mathrm{cos}\widehat{q}+i\widehat{\rho}\widehat{v}\phantom{\rule{0.16em}{0ex}}\mathrm{sin}\widehat{q}]}\hfill \\ & & \times \left[{\mu}_{\mathrm{O}}\sum _{n=1}^{{N}_{\mathrm{O}}}{e}^{-i{\mathbf{k}}_{\mathrm{f}}\xb7{\mathbf{r}}_{\mathbf{n}}}+{\mu}_{\mathrm{D}}\sum _{n={N}_{O}+1}^{N}{e}^{-i{\mathbf{k}}_{\mathrm{f}}\xb7{\mathbf{r}}_{\mathrm{n}}}\right].\hfill \end{array}$$## 3.

## Simulation Method

## 3.1.

### Simulation of Spatial Organizations of Cells

In this study, OA signals were simulated from a collection of RBCs distributed in a 2-D space and illuminated by optical radiation. Simulations using 2-D distributions are computationally less expensive compared to that of 3-D as well as capable of providing physically meaningful results. For example, our group used 2-D simulations to study the OA signal properties of nonaggregated and aggregated RBCs.^{13} In the context of ultrasound scattering, the effects of size variance of various cancer cells on backscattering were also studied employing 2-D simulations.^{17} In both cases, we observed that the simulated trends were in good agreement with those of experiments. The size of the 2-D region was taken as 500 × 500 μm^{2} and was illuminated by a laser source. This region could be viewed as a 2-D crosssection of a 3-D sample consisting of different blood corpuscles. The white blood cells and platelets do not play a significant role in the OA signals. This is because they are much less in numbers than the RBCs and have negligible optical absorption at the wavelengths of interest compared to that of RBCs.^{10} Therefore, the RBCs can be hypothesized to behave as the main OA sources in blood. The 2-D sample was assumed to be occupied by a collection of spatially random distribution of RBCs at a 45% packing fraction.^{16} The packing fraction is generally termed in the literature as the hematocrit and defined by the fractional area occupied by the cells for a 2-D sample. A 45% packing fraction was chosen because this corresponds to the normal level of hematocrit in human blood.^{16}

Spatial organizations of RBCs were generated by employing a Monte Carlo algorithm referred to as the random sequential adsorption (RSA) technique.^{18} In this method, a spatial location of a cell was randomly proposed within the 2-D region and accepted if it did not overlap with the existing cells under periodic boundary conditions. If it did not maintain the nonoverlapping condition, then this move was canceled and a new position was proposed. In this way, position coordinates of the cells were generated. Once the spatial positioning was completed, some of the cells (randomly chosen) were considered as RBCOs (*N*
_{O}) and other cells were treated as RBCDs (*N*
_{D}). Accordingly, the oxygen saturation (SO_{2}) level of that tissue realization was determined by defining, SO_{2} = *N*
_{O}/(*N*
_{O} + *N*
_{D}). As examples, possible random arrangements of RBCs are displayed in Figs. 2a, 2b. A smaller [compared to the actual size (i.e., 500 × 500 μm^{2})] 2-D region is shown for figure clarity. Each filled circle represents a RBC in these figures. In Fig. 2a, 30% cells are RBCOs (red circles) and 70% cells are RBCDs (cyan circles). However, in Fig. 2b these numbers are 80 and 20%, respectively.

## 3.2.

### Simulation Parameters for Cells

In this study, each erythrocyte was approximated as a homogeneous fluid sphere. The volume of the equivalent sphere of a RBC was fixed to 87 μm^{3} with radius *a* = 2.75 μm.^{13} The density and speed of sound of the fluid medium inside a RBC were taken as ρ_{s} = 1092 kg/m^{3} and *v*
_{s} = 1639 m/s, respectively.^{19} The same quantities for the surrounding fluid medium (blood plasma) were chosen as ρ_{f} = 1005 kg/m^{3} and *v*
_{f} = 1498 m/s, respectively.^{19} Moreover, the attenuation (due to scattering by the cells) of the light beam during its propagation through the medium was neglected. Consequently, it was assumed that the incident light intensity and the fluence (taken as, *F* = 1 J m^{−2}) were same for all cells. Furthermore, the thermal properties of the cells were considered as constants (β = 1 K^{−1} and *C*
_{p} = 1 J kg^{−1} K^{−1}) in this study. It may be noted from Eqs.
7, 8 that these physical parameters would control the amplitude of the OA signals without affecting the frequency content of the signals.

Because the oxygen saturation of hemoglobin molecules were different for RBCOs and RBCDs, the light absorption coefficients were not generally similar for these two populations (except at the isosbestic points). As a result, the emitted OA signal strengths would differ. In this study, it was assumed that oxygen saturation level was 100% for RBCOs. That means, four oxygen molecules were attached to each hemoglobin molecule in a RBCO forming an oxyhemoglobin molecule. Approximately 280×10^{6} hemoglobin molecules are packaged within a volume of 87 μm^{3} in an erythrocyte,^{20} and accordingly, the concentration of hemoglobin molecules is nearly 5.34 × 10^{−3} moles/l. Therefore, the absorption coefficient of the medium inside a RBCO could be computed by multiplying the molar concentration of oxyhemoglobin molecules and the molar extinction coefficient at a light wavelength of interest.^{21, 22} Similarly, the absorption coefficient of a RBCD could be estimated by using the appropriate molar extinction coefficient for deoxyhemoglobin molecules at that wavelength. Note that the oxygen saturation level associated with RBCDs was assumed to be 0%, meaning that hemoglobin molecules did not bind to any oxygen molecules. The estimated absorption coefficients of RBCOs and RBCDs are presented in Table 1 at three wavelengths (700, 1000, and 1064 nm). The first two wavelengths were chosen because, at these wavelengths, optical absorption coefficients of RBCOs and RBCDs would differ significantly providing distinguishable simulated signal strengths. The third wavelength was considered to compare the simulation results to that of experiments, which employed a laser source operating at 1064 nm.^{11}

## Table 1

Absorption coefficients of oxygenated and deoxygenated cells at various incident laser wavelengths.

Wavelength (nm) | μO (m−1) | μD (m−1) |
---|---|---|

700 | 356.64 | 2206.60 |

1000 | 1259.31 | 254.30 |

1064 | 407.97 | 49.13 |

## 3.3.

### Simulation and Analysis of Optoacoustic Signals

A computer code was written in C to simulate OA signals and ran on a remote computer cluster (Operating System: Linux CentOS 5, RAM: 16 GB, Processor: Intel Xenon E5462 Quad Core 2.8 GHz for each node). Equation 8 was evaluated by implementing the trapezoidal rule and using the above-mentioned numerical values for different parameters for a mixture of RBCOs and RBCDs distributed randomly in 2-D space. The OA pressure was computed at a large distance (≈6000 μm) from the center of a simulated tissue configuration and in the backward (i.e., opposite to the direction of light propagation) direction at different time points. The OA pressure computed at each time point using Eq. 8 was a complex quantity because Eq. 8 provided an analytic signal and for which the imaginary part is the Hilbert transform of the real part.^{23} The real parts of the OA time series data provided the radio frequency (RF) signal. The signal amplitudes, as a function of time, were obtained from the square root of the sum of squares of the real and imaginary parts. The OA signal properties were studied for a series of SO_{2} levels (0–100%). For each oxygen saturation level, 250 RF lines were computed from 250 possible tissue realizations. The average signal amplitude [indicated by (OA signal amplitude) in Figs. 3d, 4d], signal envelope histogram, and average power spectrum were obtained from those RF lines to study OA signal properties at each SO_{2} level. The generated envelope histogram was fitted with a Rayleigh probability distribution curve by using the fminsearch function of MATLAB R2009b. This function uses the Nelder–Mead simplex algorithm and optimizes the simulated histogram and the Rayleigh probability distribution function in order to obtain the best fit curve.^{24}

## 4.

## Simulation Results

A representative OA RF line is presented in Fig. 3a for a 700-nm illuminating laser source. In this case, 50% of cells in the sample were RBCOs and 50% were RBCDs. The hematocrit of the sample was fixed to 45%. The corresponding envelope histogram is shown in Fig. 3b. The best fit Rayleigh distribution curve is also presented in Fig. 3b. The mean power spectra of OA signals at several SO_{2} levels are displayed in Fig. 3c over a large frequency range (megahertz to gigahertz). The spectral power generally decreases as the SO_{2} level increases. The mean OA signal amplitude is plotted in Fig. 3d as a function of SO_{2}. The OA amplitude decreased monotonically with SO_{2}. At this wavelength, the absorption coefficient of a deoxygenated cell is greater than that of a oxygenated cell. Therefore, the OA amplitude is maximum at SO_{2} = 0%, where the sample contains only RBCDs. It is minimum at SO_{2} = 100% because constituent cells are only RBCOs. Also, the OA amplitude at SO_{2} = 0% is almost six times higher than that of 100%. This can be attributed to the fact that for this incident optical radiation μ_{O}:μ_{D} ≈ 1: 6.

A simulated RF line is shown in Fig. 4a for a sample with equal numbers of RBCOs and RBCDs (*N*
_{O}:*N*
_{D} = 1:1) and at 45% hematocrit for the incident optical radiation from a 1000-nm laser source. The signal amplitude is reduced compared to the previous case [see Fig. 3a]. At 700 nm, the RBCDs are the dominating OA sources and their OA signal strengths are greater than that of the RBCOs, which are the dominating OA sources at 1000 nm. The signal envelope histogram and the Rayleigh fitted curve are presented in Fig. 4b. The width of the Rayleigh distribution curve has decreased, and the peak has also shifted toward the left in comparison to Fig. 3b, because the signal amplitude is lower in this case. Figure 4c illustrates the power spectrum curves for various SO_{2} levels. Note that as the SO_{2} level increases spectral power also increases in the low- (<10 MHz) as well as in the high-frequency (>100 MHz) ranges. Between 10 and 100 MHz, spectral curves overlap for some cases. The variation of the OA mean signal amplitude with SO_{2} is plotted in Fig. 4d. The OA amplitude monotonically increases with SO_{2}. The absorption coefficient of RBCOs are more than that of RBCDs and that results in greater amplitude at higher SO_{2}. The OA amplitude at SO_{2} = 100% is nearly five times stronger than that of 0%, because μ_{O}:μ_{D} ≈ 5: 1.

We compared simulated data to experimental data reproduced from Fig. 8(c) of Ref. 11 (*in vitro* measurement at 15% hematocrit and 1064-nm optical source). The measured amplitudes (±SD) have been normalized by the amplitude of the fitted line at 100% oxygenation level and are plotted in Fig. 5. The simulated results (obtained at the same hematocrit and laser wavelength with that of experiment) have also been normalized by the mean OA amplitude computed at SO_{2} = 100% and plotted in Fig. 5. Both experimental and simulated data demonstrate a monotonic rise with SO_{2}. The fitted line to the experimental data approximately passes through the two simulated data points corresponding to SO_{2} = 0 and 100%. Note that the relative change in the OA signal amplitude between 0 and 100% SO_{2} is similar in both the theoretical and experimental data. The simulated trend is slightly nonlinear and the deviations from the linear fit are maximal between SO_{2} = 20–40%. Nevertheless, reasonably good agreement can be observed between theory and experiment.

The accuracy of the OA technique to estimate blood SO_{2} has also been assessed using the predictions of the simulations. The SO_{2} of a blood sample can be estimated using the following

## 9

[TeX:] \documentclass[12pt]{minimal}\begin{document} \begin{eqnarray} \rm {SO}_2=\frac{P_{\rm OA}^{\lambda _2}\epsilon _{\rm Hb}^{\lambda _1}-P_{\rm OA}^{\lambda _1}\epsilon _{\rm Hb}^{\lambda _2}}{P_{\rm OA}^{\lambda _1}\Delta \epsilon _{\rm Hb}^{\lambda _2}-P_{\rm OA}^{\lambda _2}\Delta \epsilon _{\rm Hb}^{\lambda _1}}, \end{eqnarray} \end{document} $$\begin{array}{c}\hfill {\mathrm{SO}}_{2}=\frac{{\mathrm{P}}_{\mathrm{OA}}^{{\lambda}_{2}}{\epsilon}_{\mathrm{Hb}}^{{\lambda}_{1}}-{\mathrm{P}}_{\mathrm{OA}}^{{\lambda}_{1}}{\epsilon}_{\mathrm{Hb}}^{{\lambda}_{2}}}{{\mathrm{P}}_{\mathrm{OA}}^{{\lambda}_{1}}\Delta {\epsilon}_{\mathrm{Hb}}^{{\lambda}_{2}}-{\mathrm{P}}_{\mathrm{OA}}^{{\lambda}_{2}}\Delta {\epsilon}_{\mathrm{Hb}}^{{\lambda}_{1}}},\end{array}$$^{10}Here, P

_{OA}is the mean OA signal amplitude. Furthermore, [TeX:] \documentclass[12pt]{minimal}\begin{document}$\epsilon _{{\rm HbO}_2}$\end{document} ${\epsilon}_{{\mathrm{HbO}}_{2}}$ and ε

_{Hb}are the known molar extinction coefficients for oxyhemoglobin and deoxyhemoglobin molecules, and [TeX:] \documentclass[12pt]{minimal}\begin{document}$\Delta \epsilon _{\rm Hb}\break=\epsilon _{\rm HbO_2}-\epsilon _{\rm Hb}$\end{document} $\Delta {\epsilon}_{\mathrm{Hb}}={\epsilon}_{{\mathrm{HbO}}_{2}}-{\epsilon}_{\mathrm{Hb}}$ . The superscripts denote the laser wavelengths at which these quantities are measured. In this study, we used 700 and 1000 nm wavelengths. The estimated SO

_{2}is plotted against the actual SO

_{2}and presented in Fig. 6. An identity line is also shown in Fig. 6. The OA method using simulated amplitudes overestimates SO

_{2}between 0 and 50% and underestimates between 50 and 100%. The nonlinear variations of the OA amplitude with SO

_{2}at both laser wavelengths might have contributed to this pattern.

## 5.

## Discussion

This paper discusses the generalized version of a theoretical model describing the OA field generated by a collection of light-absorbing objects with different physical properties and heterogeneous sizes. The model has been used here to examine how OA field would vary with SO_{2} level of a blood sample. In this work, we only considered mixtures of RBCOs and RBCDs. In such a mixture, oxygen saturation level of a single RBC was treated either as 100% or 0% and no intermediate oxygen saturation states were considered. However, in reality, a hemoglobin molecule can bind to any number (ranging from 4 to 1) of oxygen molecules leading to a fully or partially oxygen saturated state and, accordingly, this would define the absorption coefficient of that RBC. The model presented here can also accommodate such situations because it takes into account the optical properties of individual cells.

Each RBC was considered as an OA source and as it was excited by a δ function heating pulse; therefore, the emitted pressure pulse would contain all frequency pressure waves. The resultant pressure field for a sample composed of many cells randomly distributed in space was simulated by summing the pressure waves that originated from those cells. Hence, overall signal patterns were governed by the interference of waves. Also, note that both constructive and destructive interference took place for frequencies for which wavelengths of the pressure waves were smaller than the sample length. For such frequencies, phase variations between the OA signals from individual RBCs would be random (ranging from 0 to 2π) and that might have contributed to the origin of variations in the signal envelope that give rise to speckles in the signals [see Figs. 3a and 4a]. That is why the Rayleigh distribution also provided good fits to the histograms. However, experimental OA signals (e.g., Fig. 1 of Ref. 12) are generally free from speckle and consisted of a strong peak corresponding to the sample front surface (where most of the light energy was absorbed) when measured by a transducer with ≈2 MHz as the center frequency.^{12} This may be due to fact that the finite receiving bandwidth associated with the transducer filtered out high-frequency components, which may result in speckle-free signals.

It can be noted from Fig. 3c that the power spectral curves for SO_{2} = 0 and 100% looked very similar, although their magnitudes were different. For each SO_{2} level, the interfering pressure waves had identical strength because they originated from cells with the same oxygen saturation state. Therefore, contributions of the interference term (i.e.,
[TeX:]
\documentclass[12pt]{minimal}\begin{document}$\sum _{n=1}^{N} e^{-i {\bf k}_{\rm f} \cdot {\bf {r}_n}}$\end{document}
${\sum}_{n=1}^{N}{e}^{-i{\mathbf{k}}_{\mathrm{f}}\xb7{\mathbf{r}}_{\mathbf{n}}}$
) arising from the random spatial locations of cells became comparable in both cases. That resulted in similar spectral slopes between two power spectral curves over various frequency bandwidths [see Fig. 3c]. However, the magnitudes of those power spectral curves at each frequency would differ because the absorption coefficients of RBCOs (for SO_{2} = 100%) and RBCDs (for SO_{2} = 0%) were different for the 700-nm laser source. The power at each frequency for SO_{2} = 0% was at least 36 times greater than that of SO_{2} = 100% because, in this case, μ_{O}:μ_{D} ≈ 1: 6. However, for other saturation levels, the cells were not only randomly positioned in space but also the strengths of emitted pressure waves from RBCOs and RBCDs were different. Thus, power spectral curves exhibited different spectral slopes at some bandwidths compared to the curves corresponding to SO_{2} = 0 and 100%. Similar spectral characteristics could also be observed from Fig. 4c for the 1000-nm laser source.

The monotonic rise of the OA signal mean amplitude with SO_{2} is in agreement with that of an *in vitro* experimental work carried out by Esenaliev
^{11} for a laser source operating at 1064 nm; as shown in Fig. 5. In another study,^{12} the same group presented *in vivo* experimental results. Blood oxygenation was estimated using OA signals measured from external jugular vein in the neck area of a sheep at three laser wavelengths. The OA signal amplitude measured at 700 nm and normalized by that of 800 nm exhibited a linear decrease when the blood oxygenation increased from ≈20 to 95% [Figs. 3(a) and 3(b) in Ref. 12]. In contrast, the OA signal amplitude obtained at 1064 nm and normalized by that of 800 nm showed a linear increase for the same variation of blood oxygenation [see Figs. 3(a) and 3(b) in Ref. 12). Similar trends can also be observed in this study [Figs. 3d, 4d]. However, simulated patterns demonstrated slight nonlinear variations as a function of SO_{2}, whereas experimental trends were obtained by providing linear fits to the measured values. The simulated mean amplitude were obtained from signals containing a wide range of frequencies (megahertz to gigahertz). Therefore, it would not be straightforward to speculate how simulated signal amplitude for each bandwidth of interest would vary with SO_{2}. Nevertheless, at very low frequencies a coherent addition of signals might be anticipated leading to an expected linear response (also clear to some extent from the spectral lines) and that was not examined in detail in this study.

It would be of interest in future to incorporate attenuation of the optical radiation and finite nature of the ultrasound receiver bandwidth within this theoretical model and study the effects of these factors on OA signals. Further investigation is also required to clarify the nonlinear variations of theoretical curves with SO_{2}. Another important study that might be carried out in the future is to simulate signals from a collection of RBCs distributed in 3-D space and then to investigate the sensitivity and specificity of the OA technique to estimate total hemoglobin concentration in blood.

## 6.

## Summary and Conclusions

The effects of erythrocyte oxygenation on OA signals have been studied using a theoretical model. The theoretical model presented here has been developed by extending a recently proposed theoretical framework.^{13} In this model, each erythrocyte was approximated as a fluid sphere and its optical absorption property was defined by the oxygenation states of the constituent hemoglobin molecules. It was assumed that oxygen saturation level was 100% for oxygenated cells and 0% for deoxygenated cells. The oxygen saturation level of a cell controlled the strength of its OA emission. The field generated by each cell was obtained by solving the wave equation for OA field in the frequency domain and by imposing the appropriate boundary conditions.^{14} The resultant pressure field and, thus, the OA RF signal for a sample composed of a collection of cells was simulated by using the linear superposition principle for the pressure waves emitted by the individual cells. A Monte Carlo algorithm was implemented to generate spatially random distributions of RBCOs and RBCDs. Both signal envelope statistics and spectral features for such samples were examined as a function of SO_{2} level.

It was found that that the simulated mean OA signal amplitude decreased monotonically as the SO_{2} level increased for the 700-nm optical source. The same quantity for the 1000-nm incident laser radiation exhibited a monotonic rise as the SO_{2} level increased. The computed OA amplitude demonstrated nearly sixfold decrease and fivefold increase, respectively, at those wavelengths when SO_{2} level varied from 0 to 100%. These results are consistent with the experimental results. Spectral power in the low-frequency range (<10 MHz) also decreased for the 700-nm laser source and increased for the 1000-nm optical radiation with increasing SO_{2}. However, these trends were not distinctively observed between 10 and 100 MHz at SO_{2} levels >30% for the 1000-nm laser. The theoretical model presented here can be used to study OA signal properties of blood at different oxygen saturation levels.

## Acknowledgments

The financial support was provided by the Canadian Institutes of Health Research (CIHR) (Grant No. MOP-97959), CIHR grant from the Terry Fox Foundation, and the Canada Research Chairs Program awarded to M. C. Kolios. We thank “Réseau Québécois de Calcul de Haute Performance” (RQCHP), for allowing us to use the computing platform for simulations. The authors also acknowledge Dr. Min Rui, Dr. Sankar Narshimhan, and Eno Hysi for many helpful discussions and questions.

## References

*in vivo*imaging,” Nat. Biotechnol., 24 (7), 848 –851 (2006). 10.1038/nbt1220 Google Scholar

*In vivo*, noninvasive, label-free detection and eradication of circulating metastatic melanoma cells using two-color photoacoustic flow cytometry with a diode laser,” Cancer Res., 69 (20), 7926 –7934 (2009). 10.1158/0008-5472.CAN-08-4900 Google Scholar

*in vivo*,” IEEE Trans. Med. Imaging, 24 (4), 436 –440 (2005). 10.1109/TMI.2004.843199 Google Scholar

*In vivo*monitoring of blood oxygenation in large veins with a triple-wavelength optoacoustic system,” Opt. Express, 15 (24), 16261 –16269 (2007). 10.1364/OE.15.016261 Google Scholar

*in vivo*by picosecond time-resolved spectrophotometry,” J. Biochem., 109 455 –461 (1991). Google Scholar