1 November 2011 Effects of erythrocyte oxygenation on optoacoustic signals
Author Affiliations +
J. of Biomedical Optics, 16(11), 115003 (2011). doi:10.1117/1.3655355
A theoretical model examining the effects of erythrocyte oxygenation on optoacoustic (OA) signals is presented. Each erythrocyte is considered as a fluid sphere and its optical absorption is defined by its oxygen saturation state. The OA field generated by a cell is computed by solving the wave equation in the frequency domain with appropriate boundary conditions. The resultant field from many cells is simulated by summing the pressure waves emitted by individual cells. A Monte Carlo algorithm generates 2-D spatially random distributions of oxygenated and deoxygenated erythrocytes. Oxygen saturation levels of oxygenated cells a assumed to be 100% and 0% for deoxygenated cells. The OA signal amplitude decreases monotonically for the 700-nm laser source and increases monotonically for 1000 nm optical radiation when blood oxygen saturation varies from 0 to 100%. An approximately sixfold decrease and fivefold increase of the OA signal amplitude were computed at those wavelengths, respectively. The OA spectral power in the low-frequency range (<10 MHz) and in the very high-frequency range (<100 MHz) decreases for 700 nm and increases for 1000 nm with increasing blood oxygen saturation. This model provides a theoretical framework to study the erythrocyte oxygenation-dependent OA signals.
Saha and Kolios: Effects of erythrocyte oxygenation on optoacoustic signals



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.




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


[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} 2p1vs22pt2=βCPHt,
where β and 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


[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} 2p+k2p=iωμβI0CPinsidetheabsorber,0outsidetheabsorber,
where 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


[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} pfsingle(r,kf)=iμβI0vsa2CPr[sinq̂q̂cosq̂]eikf(ra)q̂2[(1ρ̂)(sinq̂/q̂)cosq̂+iρ̂v̂sinq̂],
where [TeX:] \documentclass[12pt]{minimal}\begin{document}$\hat{q}=\omega a/v_{\rm s}=k_{\rm s}a$\end{document} q̂=ωa/vs=ksa , [TeX:] \documentclass[12pt]{minimal}\begin{document}$\hat{\rho }=\rho _{\rm s}/\rho _{\rm f}$\end{document} ρ̂=ρs/ρf , [TeX:] \documentclass[12pt]{minimal}\begin{document}$\hat{v}=v_{\rm s}/v_{\rm f}$\end{document} v̂=vs/vf , and 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


[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} pfsingle(r,t)=iμβFvsa22πCPr×dω[sinq̂q̂cosq̂]q̂2[(1ρ̂)(sinq̂/q̂)cosq̂+iρ̂v̂sinq̂]eikf(ravft),
where F is the optical radiation fluence.


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


[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} pfensemble(r,kf)=n=1NiμnβnI0vsnan2CPn[sinq̂nq̂ncosq̂n]q̂n2[(1ρ̂n)(sinq̂n/q̂n)cosq̂n+iρ̂nv̂nsinq̂n]eikf(|rrn|an)|rrn|n=1NiμnβnI0vsnan2CPnr[sinq̂nq̂ncosq̂n]eikf(ran)q̂n2[(1ρ̂n)(sinq̂n/q̂n)cosq̂n+iρ̂nv̂nsinq̂n]eikf·rn.
The superscript ensemble refers to a collection of absorbers. The subscript 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


[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} pfensemble(r,t)dωn=1NiμnβnFvsnan22πCPnr[sinq̂nq̂ncosq̂n]eikf(ranvft)q̂n2[(1ρ̂n)(sinq̂n/q̂n)cosq̂n+iρ̂nv̂nsinq̂n]eikf·rn.
This is the general expression of the time-dependent pressure field generated by a collection of particles distributed in space with different sizes and physical properties.

Fig. 1

Schematic diagram illustrating the OA setup geometry.


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


[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} pfensemble(r,kf)iβI0vsa2CPr[sinq̂q̂cosq̂]eikf(ra)q̂2[(1ρ̂)(sinq̂/q̂)cosq̂+iρ̂v̂sinq̂]μOn=1NOeikf·rn+μDn=NO+1Neikf·rn,
where μO and μD are the absorption coefficients for the RBCOs and RBCDs, respectively. Here, N O and N D = NN O are the numbers of RBCOs and RBCDs, respectively. Furthermore, for a mixture of such cells, Eq. 6 simplifies into


[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} pfensemble(r,t)iβFvsa22πCPr×dω[sinq̂q̂cosq̂]eikf(ravft)q̂2[(1ρ̂)(sinq̂/q̂)cosq̂+iρ̂v̂sinq̂]×μOn=1NOeikf·rn+μDn=NO+1Neikf·rn.
The integral in Eq. 8 has been evaluated numerically in this work to compute the OA fields generated by blood samples having different proportions of RBCOs and RBCDs.


Simulation Method


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 μm2 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 (SO2) level of that tissue realization was determined by defining, SO2 = 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 μm2)] 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.

Fig. 2

(a) Simulated spatial organization of red blood cells at 45% packing fraction. Here, 30% cells are oxygenated (red circles) and 70% cells are deoxygenated (cyan circles). (b) Same as (a), but the numbers of oxygenated and deoxygenated corpuscles are 80 and 20%, respectively.



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 μm3 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/m3 and v s = 1639 m/s, respectively.19 The same quantities for the surrounding fluid medium (blood plasma) were chosen as ρf = 1005 kg/m3 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×106 hemoglobin molecules are packaged within a volume of 87 μm3 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)


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 SO2 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 SO2 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

Fig. 3

(a) Representative OA RF line computed from a sample with N O: N D = 1: 1 at 45% hematocrit and illuminated by a 700-nm laser source, (b) plots of the envelope histogram using 250 RF lines and the best fit Rayleigh distribution curve, (c) power spectral lines at different SO2 levels, and (d) variation of mean OA signal amplitude (±SD) with SO2.


Fig. 4

(a) Representative OA RF line computed from a sample with N O: N D = 1: 1 at 45% hematocrit and illuminated by a 1000-nm laser source, (b) plots of the envelope histogram using 250 RF lines and the best fit Rayleigh distribution curve, (c) power spectral lines at different SO2 levels and (d) variation of mean OA signal amplitude (±SD) with SO2.



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 SO2 levels are displayed in Fig. 3c over a large frequency range (megahertz to gigahertz). The spectral power generally decreases as the SO2 level increases. The mean OA signal amplitude is plotted in Fig. 3d as a function of SO2. The OA amplitude decreased monotonically with SO2. 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 SO2 = 0%, where the sample contains only RBCDs. It is minimum at SO2 = 100% because constituent cells are only RBCOs. Also, the OA amplitude at SO2 = 0% is almost six times higher than that of 100%. This can be attributed to the fact that for this incident optical radiation μOD ≈ 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 SO2 levels. Note that as the SO2 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 SO2 is plotted in Fig. 4d. The OA amplitude monotonically increases with SO2. The absorption coefficient of RBCOs are more than that of RBCDs and that results in greater amplitude at higher SO2. The OA amplitude at SO2 = 100% is nearly five times stronger than that of 0%, because μOD ≈ 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 SO2 = 100% and plotted in Fig. 5. Both experimental and simulated data demonstrate a monotonic rise with SO2. The fitted line to the experimental data approximately passes through the two simulated data points corresponding to SO2 = 0 and 100%. Note that the relative change in the OA signal amplitude between 0 and 100% SO2 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 SO2 = 20–40%. Nevertheless, reasonably good agreement can be observed between theory and experiment.

Fig. 5

Comparison of theoretical and experimental results. Experimental data (±SD) were reproduced from Fig. 8(c) of Ref. 11. Theoretical data were generated at 15% hematocrit and for a 1064-nm laser source to match the experiment.


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


[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} SO2=POAλ2εHbλ1POAλ1εHbλ2POAλ1ΔεHbλ2POAλ2ΔεHbλ1,
if the OA signals at two laser wavelengths are known.10 Here, POA is the mean OA signal amplitude. Furthermore, [TeX:] \documentclass[12pt]{minimal}\begin{document}$\epsilon _{{\rm HbO}_2}$\end{document} εHbO2 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} ΔεHb=εHbO2ε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 SO2 is plotted against the actual SO2 and presented in Fig. 6. An identity line is also shown in Fig. 6. The OA method using simulated amplitudes overestimates SO2 between 0 and 50% and underestimates between 50 and 100%. The nonlinear variations of the OA amplitude with SO2 at both laser wavelengths might have contributed to this pattern.

Fig. 6

Plot of the OA estimation of blood SO2 as a function actual SO2.




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 SO2 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 SO2 = 0 and 100% looked very similar, although their magnitudes were different. For each SO2 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} n=1Neikf·rn ) 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 SO2 = 100%) and RBCDs (for SO2 = 0%) were different for the 700-nm laser source. The power at each frequency for SO2 = 0% was at least 36 times greater than that of SO2 = 100% because, in this case, μOD ≈ 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 SO2 = 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 SO2 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 SO2, 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 SO2. 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 SO2. 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.


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 SO2 level.

It was found that that the simulated mean OA signal amplitude decreased monotonically as the SO2 level increased for the 700-nm optical source. The same quantity for the 1000-nm incident laser radiation exhibited a monotonic rise as the SO2 level increased. The computed OA amplitude demonstrated nearly sixfold decrease and fivefold increase, respectively, at those wavelengths when SO2 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 SO2. However, these trends were not distinctively observed between 10 and 100 MHz at SO2 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.


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.


1.  J. K. J. Li, T. Wang, and H. Zhang, “Rapid noninvasive continuous monitoring of oxygenation in cerebral ishemia and hypoxia,” Cardiovasc Eng. 10, 213–217 (2010). 10.1007/s10558-010-9116-x Google Scholar

2.  L. V. Wang, “Prospects of photoacoustic tomography,” Med. Phys. 35, 5758–5767 (2008). 10.1118/1.3013698 Google Scholar

3.  H. F. Zhang, K. Maslov, G. Stoica, and L. V. Wang, “Funtional photoacoustic microscopy for high-resolution and noninvasive in vivo imaging,” Nat. Biotechnol. 24(7), 848–851 (2006). 10.1038/nbt1220 Google Scholar

4.  M. Eghtedari, M. Motamedi, V. L. Popov, N. A. Kotov, and A. A. Oraevsky, “Optoacoustic imaging of gold nanoparticles targeted to breast cancer cells,” Proc. SPIE 5320, 21–28 (2004). 10.1117/12.534464 Google Scholar

5.  R. M. Weight, J. A. Viator, P. S. Dale, C. W. Caldwell, and A. E. Lisle, “Photoacoustic detection of metastatic melanoma cells in the human circulatory system,” Opt. Lett. 31(20), 2998–3000 (2006). 10.1364/OL.31.002998 Google Scholar

6.  E. I. Galanzha, E. V. Shashkov, P. M. Spring, J. Y. Suen, and V. P. Zharov, “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

7.  J. A. Viator, G. Au, G. Paltauf, S. L. Jacques, S. A. Prahl, H. Ren, Z. Chen, and J. S. Nelson, “Clinical testing of a photoacoustic probe for port wine stain depth determination,” Lasers Surg. Med. 30(2), 141–148 (2002). 10.1002/lsm.10015 Google Scholar

8.  M. Rui, W. Bost, E. Weiss, R. Lemor, and M. C. Kolios, “Photoacoustic microscopy and spectroscopy of individual red blood cells,” in Biomedical Optics, OSA Technical Diegest (CD), Optical Society of America, Paper No. BSuD93 (2010). Google Scholar

9.  J. J. Niederhauser, M. Jaeger, R. Lemor, P. Weber, and M. Frenz, “Combined ultrasound and optoacoustic system for real-time high-contrast vascualr imaging in vivo,” IEEE Trans. Med. Imaging 24(4), 436–440 (2005). 10.1109/TMI.2004.843199 Google Scholar

10.  X. Wang, X. Xie, G. Ku, L. V. Wang, and G. Stoica, “Noninvasive imaging of hemoglobin concentration and oxygenation in the rat brain using high-resolution photoacoustic tomography,” J. Biomed. Opt. 11(2), 024015 (2006). 10.1117/1.2192804 Google Scholar

11.  R. O. Esenaliev, I. V. Larina, K. V. Larin, D. J. Deyo, M. Motamedi, and D. S. Prough, “Optoacoustic technique for noninvasive monitoring of blood oxygenation: a feasibility study,” Appl. Opt. 41(22), 4722–4731 (2002). 10.1364/AO.41.004722 Google Scholar

12.  H. P. Brecht, D. S. Prough, Y. Y. Petrov, I. Patrikeev, I. Y. Petrova, D. J. Deyo, I. Cicenaite, and R. O. Esenaliev, “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

13.  R. K. Saha and M. C. Kolios, “A simulation study on photoacoustic signals from red blood cells,” J. Acoust. Soc. Am. 129(5), 2935–2943 (2011). 10.1121/1.3570946 Google Scholar

14.  G. J. Diebold, “Photoacoustic monopole radiation: waves from objects with symmetry in one, two and three dimensions,” Chapter 1 in Photo-acoustic Imaging and Spectroscopy, L. V. Wang, Ed., pp. 3–17, Taylor and Francis (2009). Google Scholar

15.  A. Ishimaru, Chapter 4 in Wave Propagation and Scattering in Random Media, pp. 69–92, Academic Press, New York (1978). Google Scholar

16.  K. K. Shung and G. A. Thieme, Ultrasound Scattering in Biological Tissues, pp. 1–499, CRC Press, Boca Raton (1993). Google Scholar

17.  R. M. Vlad, R. K. Saha, N. M. Alajez, S. Ranieri, G. J. Czarnota, and M. C. Kolios, “An increase in cellular size variance contributes to the increase in ultrasound backscatter during cell death,” Ultrasound Med. Biol. 36(9), 1546–1558 (2010). 10.1016/j.ultrasmedbio.2010.05.025 Google Scholar

18.  E. L. Hinrichsen, J. Feder, and T. Jøssang, “Random packing of disks in two dimensions,” Phys. Rev. A 41, 4199–4209 (1990). 10.1103/PhysRevA.41.4199 Google Scholar

19.  K. K. Shung, Y. W. Yuan, D. Y. Fei, and J. M. Tarbell, “Effect of flow disturbance on ultrasonic backscatter from blood,” J. Acoust. Soc. Am. 75, 1265–1272 (1984). 10.1121/1.390733 Google Scholar

20.  S. I. Fox, Chapter 13 in Human Physiology, pp. 388–430, 10th ed., McGraw-Hill, New York (2008). Google Scholar

21.  S. Prahl, Tabulated data from various sources,  http://omlc.ogi.edu/spectra, accessed March 10, 2011. Google Scholar

22.  Y. Nomura and M. Tamura, “Quantitative analysis of the hemoglobin oxygenation state of rat brain in vivo by picosecond time-resolved spectrophotometry,” J. Biochem. 109, 455–461 (1991). Google Scholar

23.  R. N. Bracewell, “Relatives of the Fourier transform,” Chapter 12 in The Fourier Transform and Its Applications, 2nd ed., pp. 241–274, McGraw-Hill Kogakusha Ltd., Tokyo (1978). Google Scholar

24.  J. C. Lagarias, J. A. Reeds, M. H. Wright, and P. E. Wright, “Convergence properties of the Nelder-Mead simplex method in low dimensions,” SIAM J. Optimi. 9 (1), 112–147 (1998). 10.1137/S1052623496303470 Google Scholar

© 2011 Society of Photo-Optical Instrumentation Engineers (SPIE)
Ratan K. Saha, Michael C. Kolios, "Effects of erythrocyte oxygenation on optoacoustic signals," Journal of Biomedical Optics 16(11), 115003 (1 November 2011). https://doi.org/10.1117/1.3655355



Monte Carlo methods


Laser sources


Computer simulations

Back to Top