Photoacoustic and acousto-optic tomography for quantitative and functional imaging

Optical absorption contrast, large imaging depth at ultrasonic resolution, and the potential of functional/quantitative imaging are the features driving the rapid development of photoacoustic (PA) imaging. For quantitative and functional PA imaging, the fluence distribution is required to transform a PA image to a map of absorption coefficients. A suitable method to estimate the fluence for in vivo applications must not require a priori knowledge of the medium optical properties, should work for various illumination conditions, and must be applicable to different PA imaging geometries. Existing methods of estimating fluence in tissue do not meet these requirements simultaneously. Here we present a method to measure the fluence distribution in tissue using acousto-optics (AO) that meets all the above requirements. We developed a PA and AO tomography system for small-animal imaging to investigate the potential and the feasibility of fluence-corrected PA imaging using our method in a single instrument. We performed experiments on phantoms, an ex vivo tissue sample, and freshly sacrificed mice. We demonstrate that the correction for spatial and spectral fluence variations in PA images establishes the direct relation between image value and optical absorption, which in turns improves the quantitative estimation of blood oxygen saturation. © 2018 Optical Society of


INTRODUCTION
Photoacoustic (PA) imaging bridges the resolution gap between optical microscopy and diffuse optical tomography (DOT).It benefits from the rich optical absorption contrast of tissue and offers ultrasonic resolution at an imaging depth comparable to DOT.PA imaging has several variants, such as photoacoustic tomography (PAT), with applications in small-animal imaging for preclinical research [1][2][3], mammography [4][5][6], imaging of extremities of the human body [7], endoscopy [8], and hand-held systems in combination with linear array ultrasound (US) systems [9].In PA imaging, tissue is excited by a short (∼ns) laser pulse; as a result, the local absorption of light by tissue chromophores leads to a local pressure rise, which subsequently relaxes by the emission of broadband acoustic waves.These acoustic waves then propagate to the surface of the medium and are measured using US transducers.Assuming a perfect acoustic reconstruction method, the initial pressure distribution σ o can be formed based on the measured acoustic signals, which reads σ o r, λ ΓE a λ, r Γφλ, rμ a r, λ, where r is the position, λ is the wavelength of the excitation light, Γ is the Gruneisen parameter, φ is the optical light fluence depending on the μ a and μ 0 s of the surrounding tissue, with μ 0 s the reduced scattering coefficient and μ a the absorption coefficient.The absorption coefficient μ a r, λ P n i C i rε i λ is the linear combination of the concentrations C of the n chromophores present with extinction coefficients ε.
As seen from Eq. ( 1), the PA image is proportional to the absorbed energy density E a , which is a product of the absorption coefficient μ a (r, λ) and the optical fluence distribution φr, λ in the sample.The optical fluence distribution is an unknown, which presents a challenge in quantifying the image.The quantification refers to transformation of the PA image into concentration maps of the constituent chromophores, either in absolute terms or percentiles of constituent chromophores.A lower degree of quantification is to establish maps of the relative or absolute absorption coefficients.Given a certain illumination configuration, the optical fluence distribution depends on the distribution of the absorption and scattering coefficients in the sample, which has two major implications.The first is the so-called "spectral coloring;" as a result of wavelength-dependent absorption and scattering coefficients, the spectrum of the fluence at a point inside the sample is different from the incident spectrum.The second is "absorption and scattering nonuniqueness," meaning it is not possible to recover absorption and scattering distributions uniquely and simultaneously from a single measurement.As a result, it is not possible to distinguish whether the difference in calculated (i.e., based on a model) and measured absorbed energy density (E a ) is due to an error in the local absorption coefficient or in the background scattering coefficient.
If the fluence φr, λ distribution is known, multispectral PA imaging can provide a direct measurement of the spectrally and spatially resolved concentrations of tissue chromophores.Quantification efforts in PA imaging, generally speaking, can be divided into two categories: the model-based inversion schemes and the experimental methods [10].The model-based inversion schemes aim to recover the absorption coefficient μ a r, λ using measured PA data in combination with light transport models for estimation of the light fluence distribution [11][12][13][14][15][16][17][18][19][20][21], with the exception [22,23] where authors proposed to directly recover the absorption coefficient using a stochastic approach and [24] where authors modeled the light fluence as an affine function of reference spectra.Widely used light transport models include the radiative transfer equation (RTE), the diffusion equation (i.e., simplified form of RTE in the diffusive regime), and numerical simulations based on the Monte Carlo method.The inversion schemes that make assumptions regarding optical properties of the medium may have limited scope when applied to tissue [25].Common assumptions listed in the order of crude to more sophisticated include: background scattering and bulk absorption coefficients are known [11]; fluence distribution is independent of wavelength [12]; bulk scattering and absorption coefficients are known and the absorption and/or scattering perturbations are small [13,16]; the influence of these perturbations is insignificant on the fluence distribution [17,19]; and the wavelength dependence of the scattering is known [18].Experimental approaches aim at fluence correction in PA imaging without the need of a priori knowledge about the medium optical properties.Bauer et al. [26] use DOT to first measure the absorption and scattering of the sample and then use it to calculate the fluence distribution.Dean Ben et al. use fluence-dependent photo-switchable contrast agents [25].Previously, we have shown two experimental approaches to solve the problem of fluence compensation in PA.The first approach combines a transmission mode acousto-optic (AO) measurement with two PA measurements and provides fluence compensation of PA signals [27].The two PA measurements are done by subsequently illuminating the medium from excitation and detection points during AO measurements.The second approach enables a direct measurement of fluence variation using reflection mode AO only, for a single-point illumination and detection [28].Since experimental methods do not require a priori knowledge of optical properties of the medium, it makes them potentially advantageous for model-based approaches.However, the methods we have presented so far have limited application in real-world imaging systems, since the validity of these methods is limited to either transmission or reflection mode experimental configurations.Furthermore, since the allowed illumination and detection is limited to one or two single points, the optical etendue of the conceived AO detection system will be inherently limited.
Here we present a method to measure the fluence distribution using ultrasonic modulation of light that is general in terms of illumination conditions, and therefore, it is independent of the experimental geometry.Further, we present a combined PA and AO tomography system for small-animal imaging to demonstrate the feasibility of our method to enable fluence-corrected PA imaging.

A. Theoretical Aspect, Measurement of Optical Fluence Inside Highly Scattering Media
Acousto-optics, a hybrid of light and US, is proposed as a tissue imaging modality [29], where localization is achieved using focused US.US experiences negligible scattering in tissue compared to light, allowing the formation of a well-defined focus at a desired location.Upon illuminating the tissue with coherent light, fractions of light that escape from the medium will form speckle patterns.Light that interacts with the US focus gets phase-modulated due to refractive index and scattering displacement modulation caused by the pressure gradient of US [30].The nonmodulated and frequency-modulated light diffuses farther to the boundary of the medium.Light modulation causes a fluctuation of the speckle patterns, from which the modulated part can be isolated using different detection methods [31][32][33][34][35]. Thus, US enables noninvasive tagging of light at a region of interest (ROI) inside the tissue.
Here we describe the mathematical formulation for estimating the fluence inside a scattering medium with unknown optical properties based on ultrasonic tagging of light.Consider a scattering medium that is simultaneously illuminated at x n surface points, such that the input power at each point is P in ; then the fluence rate φ xT at internal point x T can be written as where Prx i , x T is the probability per unit area per solid angle of a photon starting at point i to cross an aperture placed at x T .Here the condition of equal input power P in at all the illumination points is to balance the number of unknowns and the number of equations at a later stage of mathematical derivation.However, physically it does not result in any limitation, since by using a large number of discrete input beams of the same power, we can construct any desired light intensity distribution at the surface of the medium, as demonstrated in our experiments.
Considering that the light at internal point x T is tagged with efficiency η using focused US and reinjected into the medium, the power of the ultrasonically tagged light P L,T reinjected into the medium can be written, with the help of Eq. ( 2), as where A T is the frontal area of the tagging volume (US focus) at point x T .The concept of the frontal area is borrowed from aerodynamics and is convenient to describe the area projected by irregular shapes.The frontal area is the area of the shadow cast by irregular shapes when illuminated by parallel light rays.In our case, we use the frontal area averaged for all possible directions of illumination.The internally injected ultrasonically tagged light at power P L,T propagates further through the medium, leading to detection of ultrasonically tagged light at the boundary of the medium.The power of detected ultrasonically tagged light (P L,D ), when n detection points coincide with n excitation points, and when each detection point has an area A and solid opening angles Ω, can be written as The first term on the right-hand side of Eq. ( 4) accounts for the diffusely reflected tagged light (i.e., light injected at i, tagged at T , and detected at i); the second term when i ≠ j is the tagged transmitted light (light injected at i, tagged at T , and detected at j).
According to the photon path reversibility [27], we can write Considering the photon path reversibility, and using it in Eq. ( 4) results in Taking the square on both sides of Eq. (2) gives Comparison of Eqs. ( 6) and ( 7) reveals that the terms in the brackets on the right-hand side are identical; hence, solving both equations together yields where P L,D is the normalized total power of the tagged light.The final expression for local fluence when illumination and detection are done at multiple surface points of the medium is identical to the one derived in [28] for point illumination and detection in reflection mode, but now is derived without restriction in terms of the number of points for light illumination and detection.The parameters η and A T in the context of ultrasonic modulation of light are difficult to quantify, which presents a challenge in measuring the light fluence in absolute terms.Assuming η and A T are constant throughout the tissue, Eq. ( 8) may provide a means to map the distribution of fluence in relative terms: However, the validity of the assumption regarding η and A T being constant throughout the tissue is critical in order to reliably transform the measured AO signal into fluence variation.The justification of these assumptions is given in Section 2.C on "Ultrasonic Tagging of Light and Speckle Contrast-Based Detection of Tagged Light."

B. Experimental Setup
Figure 1 shows a schematic of the combined PA and AO tomographic imaging system.The sample is held stationary in the imaging tank containing water for acoustic coupling.The 64-element US transducer array is custom-fabricated by IMASONIC Inc. (France); in azimuthal direction, it has an angular aperture of 19.2°, a pitch of 0.3 mm, and an element size of 0.2 × 10 mm, with a radius of curvature of 40 mm.The individual elements have a convex geometry with a focus at 40 mm in the elevation plane and a center frequency of 5 MHz with a −6 dB bandwidth of 100%.From a PA imaging perspective, the choice of the center frequency of the US transducer used defines the trade-off between imaging and depth.We choose 5 MHz, with the goal of using our system for small-animal imaging, where 5 MHz is a suitable frequency since it provides ∼150 mm resolution and 2-3 cm imaging depth.From an AO perspective, in our previous work [28], we showed that at 5 MHz it is possible to accurately measure light fluence using AO in a medium with a reduced scattering coefficient variation from 6-15 cm −1 .
A commercially available clinical US scanner (MyLab One, ESAOTE Europe) is used for transmission and acquisition of US and PA signals at sampling frequency of 50 MHz.For PA, a frequency-doubled 6 ns pulse width Nd:YAG laser (Quanta-Ray Lab, Spectra-Physics) is used to pump the optical parametric oscillator (VersaScan L-532, GWU, wavelength tuning range: 680-910 nm).For AO, we used a single-mode continuous wave (CW) laser (Coherent MBR 110 EL, Ti:Sapphire) pumped by a Coherent Verdi 6 laser.An AO modulator (Gooch & Housego, R23080-3-LTD) is used in the beam of the CW laser for converting it into a quasi-pulsed laser with pulse length of 1.5 μs light.A pulsed illumination in combination with short US bursts enables axially resolved imaging in AO [36].The light from the lasers is coupled into a custom-made (LightGuideOptics Germany GmbH) multimode optical fiber bundle (diameter = 1 cm) that splits into six smaller bundles (diameter = 0.4 cm) to illuminate the imaging plane of the sample.Each fiber bundle is 6 mm in diameter, and the bundles are placed 60°apart from each other.A polarizing beam splitter (PBS) is used to reflect the diffuse light emanating from the sample, collected by the optical fibers, on to a CCD camera for recording the speckle pattern to measure ultrasonically modulated light.For all the phantoms and ex vivo sample measurements, we used the Stingray F046 CCD camera (15 frames per second), whereas for the mouse measurements, we used the Mako G-30 CCD camera (300 frames per second).The higher frame rate offered by Mako G-30 reduces the total measurement time by a factor of 20.The US array and illuminating optical fiber bundles are attached to a single linear translation stage allowing cross-sectional imaging across the length of the sample (max = 20 cm).A LabVIEW program controls the scanning and synchronization of the US scanner, laser pulses, and the CCD camera using a function generator (Tektronix, AFG 3102).

C. Ultrasonic Tagging of Light and Speckle Contrast-Based Detection of Tagged Light
The measurement of ultrasonically modulated light in our experiments is done using a method based on measuring speckle contrast change [35].The speckle contrast C σ∕hI i is the ratio of the standard deviation σ of the light intensity in the speckle pattern to the average light intensity.The modulation of optical path lengths due to US pressure causes speckles to fluctuate in time, which when integrated over time (within the integration time of CCD) results in blurring of the speckle pattern, leading to a reduction in contrast.The signal in this method is the speckle contrast change (ΔC) between ON and OFF states of US.In typical AO experiments for a given US pressure and integration time of the CCD camera, the measured ΔC is proportional to the power of ultrasonically modulated light (P L,D , P L,D κΔC) emanating from the sample [35]; κ is proportionality constant involving experimental parameters.Hence, when using a speckle contrast-based detection method, Eq. ( 9) for estimating the fluence can be written as In our experiments, we used a US transducer array to transmit a focused US for tagging of the light.A maximum pressure of 2.7 MPa is possible to be transmitted using this probe in combination with a MyLab One US scanner.Figure 2(a) shows the measured pressure distribution using a calibrated needle hydrophone in the imaging plane of the US transducer array.In our theory, we assumed that the active frontal area of the US focus remains unchanged; this parameter is related to the shape and size of the US focus.To justify this assumption, we imaged the US focus using Schlieren imaging after it has propagated through various media.Figure 2(b) shows the Schlieren images of bursts of five US cycles after propagation through 2.6 cm of water (top), agar and Intralipid phantom (middle), and chicken breast tissue (bottom).These results show that the US focus shape and size remain unaffected after propagation through different tissue-like media.The second assumption is regarding the tagging efficiency (η) being constant throughout the tissue.Physical analysis of AO modulation reveals that η is dependent on optical scattering of the tissue as well as the local parameters of the US used for modulation of light, such as pressure and frequency [37].While different strategies can be applied to minimize the variation in US parameters, e.g., [38] inside the tissue, the scattering properties of the tissue cannot be assumed constant.However, the codependence of the η on US parameters and scattering properties of the tissue [30] might be exploited to minimize the dependence of η on scattering properties of the tissue.In our previous work, we experimentally studied the effect of scattering coefficient variations on the ability of AO tagging-based method to measure the fluence variations, where we used a 5 MHz US for the modulation of light.The results of the study show that within the range of reduced scattering coefficient of 6 cm −1 to 15 cm −1 (which covers common tissue types for PA imaging within the optical window [39]) our method can reliably measure the fluence variations.This in essence proves that dependence of η on reduced scattering coefficient within the above-mentioned range is negligible.Our hypothesis is that by choosing short wavelength US, the speckle contrast detection process can be saturated at relatively low scattering coefficient using US pressures relevant for medical US imaging.However, the validity of this assumption dictates the accuracy of the results; therefore, care must observed in choosing the US parameters and detection method.Here we study the speckle contrast change as a function of the US pressure in the focus in our setup.
We performed this experiment in a scattering phantom made of agar and Intralipid, with a reduced scattering coefficient of 6 cm −1 .We measured the change in contrast as a function of US pressure, and all the other experimental parameters, such as input light power, position of the US focus inside the medium, and the integration time of the CCD, were kept constant.Figure 2(c) shows that the measured ΔC has a quadratic dependence on the US pressure in the focus up till 2.5 MPa, while above this pressure, ΔC begins to plateau.The quadratic dependence of ΔC on US pressure is in agreement with theory [37].

PA Imaging
For PA measurements, the multiple projection data for a single cross section is acquired by rotating the probe around the sample using a stepper motor with a step size of 30°.The total number of projections in presented experiments varies, and it is stated in the corresponding section.Multiple cross sections in animal experiments are acquired by linearly translating the US array using a stepper motor.The PA images are reconstructed using a filtered back-projection algorithm [40,41], after de-convolving the data with the measured impulse response of the individual elements of each element of the US probe.The impulse response was measured by illuminating a 50 μm diameter surgical suture (black color) with the above-mentioned pulsed laser operating at 750 nm wavelength and recording the generated PA signal with the US probe connected to the MyLabOne scanner.

US Imaging
US imaging is performed using pulse echo tomography.We acquire nine projections for a given cross section with a step size of 40°.Each projection is acquired by transmitting multiple plane waves (−20°to 20°in steps of 0.5°) and then coherently adding the reconstructed images.The US images are reconstructed a using Fourier transform-based reconstruction algorithm [42], where the curved array is approximated to be linear, under the assumption that the ratio of the phased array aperture to the radius of curvature is small.Multiple projections are then combined to form an image of the object cross section.

Fluence Mapping
Mapping of fluence variations using AO is a point-by-point scanning process.The entire imaging plane is scanned by the US burst, and the resulting speckle contrast is measured using a CCD camera, which is then transformed into a relative fluence map using Eq.(10).The scanning along the propagation direction of US is done by controlling the delay between US bursts (5 cycles = 1.5 mm) and the laser pulses.The transmission delay applied to the phased array US probe allows formation of a focus at the desired depth.In our experiments, the US focus was formed at a fixed depth of 34 mm, which is 6 mm before the natural focus (40 mm) of the phased array and its center of rotation.The center of rotation is also approximately the center of the sample in our experimental setup.The scanning along the transverse direction is done by shifting the active aperture (32 elements) in steps of 0.6 mm (see Fig. 1).The active aperture means the set of transducer elements that transmits US to form the focus.With the shift in active aperture as described above, the transmission angle of the US beam is also changed such that the scan area is a rectangle of 11.5 × 18 mm 2 , instead of a segment of a circle, which would result from emitting a symmetric beam from the active aperture.This scanning approach reduces the redundant measurements close to the center of circle, thereby significantly reducing the measurement time in comparison with mechanical scanning of the US focus.
The measurement speed in this scanning approach is currently limited by the frame rate of the CCD camera.With the current camera (Mako G-030), the total measurement time in our system is 29 min to scan an area of 530 mm 2 with 30 averages and spatial sampling of 0.6 mm.

E. Estimation of Blood Oxygen Saturation
As a result of spectral coloring, the estimation of blood oxygen saturation (sO 2 ) using multiwavelength PA imaging is challenging.In our previous work [43], we demonstrated that the combination of PA and AO enables measurement of absolute oxygen saturation, with improved accuracy compared to a method based on PA alone.The method used in our previous work was based on the transmission mode AO.The method required two PA measurements per wavelength such that the PA excitation points coincided with the optodes of AO measurements.This was a very specific and not generally applicable procedure that did not yield optical fluence as a separately measured quantity.In contrast, in our current work we introduce a more generalized illumination scheme.In our current study, the sample is illuminated by six equally spaced fiber bundles as an example of the many possible schemes.Furthermore, we now directly measuring fluence variations at each wavelength and use them to correct for the effect of spectral coloring in PA images.
Assuming that the two dominant chromophores in blood are oxyhemoglobin and de-oxyhemoglobin, the absolute sO 2 of blood can be estimated using the following expression: where M σ∕Γκ 2 ffiffiffiffiffiffiffiffiffiffiffiffiffi ffi λ 2 κΔC p , is the fluence-normalized PA image, ε Hb is the molar extinction coefficient of de-oxyhemoglobin and Δε ε HbO − ε Hb , where ε HbO is the molar extinction coefficient of oxyhemoglobin.Indices λ 1 and λ 2 indicate the two wavelengths used.In the above expression of M , ΔC is multiplied with λ 2 to account for the dependence of tagging efficiency of US on the wavelength of light [43].The value of sO 2 ranges from 0 to 1, where 0 means the blood consists of only de-oxyhemoglobin and 1 means the blood consists of only oxyhemoglobin.

F. Blood Samples for Oxygen Saturation Measurements
Blood samples were obtained from anonymous healthy volunteers organized by the donor service at the Experimental Centre of Technical Medicine at the University of Twente.Blood samples with different sO 2 values were used for measurements.Different values of sO 2 were obtained by adding sodium dithionite (Na 2 S 2 O 4 ) to the whole blood.Oxygen saturation values of the blood samples were measured using a blood oximeter (Avoximeter 1000E, ITC point of care) before and after each experiment for comparison with photoacoustically estimated sO 2 .The measured hemoglobin concentration of the blood sample used was 15.51 0.37 g∕dL.

RESULTS AND DISCUSSION
A. Two-Dimensional Fluence-Compensated PA Imaging in Phantom Figure 3(a) shows an illustration of the cross section of a phantom with three identical absorbing inclusions embedded at different depths.The phantom is an optically scattering 26 mm diameter cylinder with a reduced scattering coefficient of 6 cm −1 and is made by mixing 3% agar (W/V) and 5.1% Intralipid (V/V) into water.Figure 3(b) shows the PA image of the corresponding cross section, acquired using 12 projections, where the PA image values from identical absorbers are different from each other due to variations in fluence within the sample.Figure 3(c) shows the measured fluence distribution in this sample; each measurement point is averaged over 30 measurements.The fluence map is acquired by scanning the US focus, with total acquisition time of 29 min.The scanning step size along the propagation direction of US as well as in transverse direction is 0.6 mm.The region where the sample is present is marked with the dotted circle.Six bright spots around the sample correspond to the illumination points of the six optical fiber bundles around the sample.The measured fluence appears to increase with depth for the first 2 mm from the surface.This is because of the finite size of the US burst; as a result, the maximum signal is measured when the complete US burst is inside the medium, and the results are only valid from this depth onward.
Figure 3(d) shows the PA image normalized with the measured fluence distribution [Fig.3(c)], where differences in image values due to fluence variations are compensated, and the PA image value linearly relates to the optical absorption coefficient.Figure 3(e) shows the PA image value along a vertical line through the center of three absorbers, before and after the fluence correction, for quantitative comparison.These results demonstrate the feasibility of a noninvasive experimental approach to measure a two-dimensional fluence distribution inside scattering media.

B. Estimation of Blood Oxygen Saturation sO 2
Figure 4(a) shows the illustration of the cross section of the phantom used.The phantom is a 26 mm diameter cylinder, where approximately half the cylinder is only scattering, and the other half is scattering and absorbing, to enhance the spectral coloring effect.The phantom is made of agar (3% W/V) and Intralipid (8% V/V), where Intralipid provides optical scattering and the optical absorption is introduced by adding 15 μL black ink (Royal Talens, Ecoline 700) in 100 mL of water.The phantom has a background-reduced scattering coefficient of 9.3 cm −1 and 8.6 cm −1 and absorption coefficient of 0.017 cm −1 and 0.013 cm −1 at 750 and 800 nm, respectively.Two nylon tubes (inner diameter 0.6 mm and outer diameter 0.9 mm) are placed in the phantom, one in nonabsorbing background (NonAbs.BG) and one in absorbing background (Abs.BG), to mimic blood vessels.Figures 4(b) and 4(c) show PA images normalized with input pulse energy at 750 and 800 nm, respectively, without the fluence correction.The PA images are acquired using three projections with step size of 30°.Both the tubes contain blood with sO 2 43%.However, the PA image value for the tube placed in NonAbs.BG is higher than for the tube in the Abs.BG due to differences in local fluence.Further, the ratio of the PA image values at two wavelengths is different for both tubes and does not agree with the expected ratio for 43% oxygenated blood.
These differences are a result of spectral coloring.Figure 4(d) shows the measured fluence variation in the Abs.BG and Non Abs.BG along the lines indicated with the solid line [Fig.4(a)].The fluence variations shown in Fig. 4(d) are corrected for wavelength-dependent tagging efficiency as described in the methods section.The measured fluence distribution clearly shows that the decay rate is dependent on the background properties of the phantom.The decay in fluence with depth is faster for 750 nm wavelength than 800 nm wavelength.Figures 4(e) and 4(f ) are the PA images corrected for fluence variations at 750 nm and 800 nm, respectively.Here the PA image values from both the tubes within an image at given wavelength are almost identical, and the effect of spectral coloring is removed.In Fig. 4(g), sO 2 of different blood samples measured using PA only and fluence-corrected PA is plotted against sO 2 measured using the oximeter.The estimation of sO 2 using PA only and  fluence-corrected PA images plotted in Fig. 4(g) is the average value within an ROI around the blood-filled tubes.The ROI around the blood-filled tubes is defined by considering the pixels with image value higher than 0.25 times the maximum value.The threshold criteria are chosen such that all the pixels within the blood-filled tube are included in the analysis.The spectral coloring causes erroneous measurement of sO 2 using PA only, on average 90% and 129% over estimation in Abs.BG and Non Abs.BG, respectively.The error is dependent on the optical properties of the background sample.The fluence correction of PA images based on our method, enables more accurate measurement of sO 2 , on average 8% over and 4% under estimation in Abs.BG and Non Abs.BG, respectively.The precision of the sO 2 estimation after fluence correction is represented by the error bars, which is 12% in Abs.BG and 7% in Non Abs.BG.The influence of the optical properties of the background sample is reduced considerably.These results demonstrate the ability of the method to correct for spectral coloring in an optically heterogeneous medium and to provide accurate measure of blood sO 2 .
Biological tissue is not only optically but also acoustically heterogeneous.The acoustic heterogeneities may lead to inaccurate reconstruction of PA images and may distort the US focus and thereby affect the accuracy of fluence correction using our method.Although the acoustic heterogeneity in the phantom experiment in the form of nylon tubes is present at the target location, the deformation in the shape and size of US focus caused by the tube might not be significant because the US focus deformation is an accumulative process with the propagation distance.To assess the feasibility and accuracy of quantitative sO 2 measurements using our method in acoustically heterogeneous tissue, we repeated the above experiments to measure blood sO 2 in porcine tissue [Fig.5(a)].The porcine tissue sample consists of layers of fat and muscles.We placed two nylon tubes to mimic blood vessels at different depths.The tissue sample is fixed in a 26 mm diameter cylinder made of agar gel. Figure 5(b) shows the US image of this sample, where both blood-filled nylon tubes are visible alongside the contrast of muscle and fat. Figure 5(c) shows an example of a fluence-compensated PA image at 750 nm of this sample.The image is acquired using four projections.In this example, the image value from the two tubes is slightly different, even after the fluence compensation.The difference may be due to a reconstruction artifact, as can be seen from the zoomed-in part of the figure; the shape of tube 2 is well reconstructed, leading to a higher image value, whereas the shape of tube 1 is more dispersed, leading to a lower image value resulting from acoustic heterogeneities.Figure 5(d) plots the measured fluence variation as a function of depth along the dotted line through the tubes [US image, Fig. 5(b)], for the two wavelengths used: 750 and 800 nm. Figure 5(e) plots the measured sO 2 based on uncompensated and fluence compensated PA against the sO 2 measured with the oximeter.The precision of the sO 2 estimation using our method is represented by the error bars.The results show that the wavelength-dependent fluence compensation of PA using our method improves the accuracy of sO 2 estimation compared to PA alone, even in tissue.The improvement in the estimation of oxygen saturation after wavelength-dependent fluence compensation of PA compared to PA alone in tissue is not as significant as in the phantom experiment.While US focus distortions might present a challenge for correcting PA images for spatial variations in fluence in tissue, its effect on the compensation of spectral coloring should be minimal.This is because the correction for spectral coloring is local, and the shape of the US focus at a given location inside the sample remains unchanged for two excitation wavelengths of light.

C. Two-Dimensional Fluence-Compensated PA Imaging in Sacrificed Mice
Further, we investigated the feasibility of acousto-optically mapping the fluence to facilitate fluence-compensated PA tomography in small animals.We used two nude mice in our experiments; mice were taken from another study at University Medical Center Groningen and properly housed and taken care of according to the guidelines of the Central Committee Animal Experimentation of the Netherlands (CCD).
Quantitative validation of the fluence compensation in small animals is challenging because of the unknown absorption coefficients of endogenous contrast sources, such as blood vessels.Therefore, to quantitatively assess the validity of our fluence compensation method in small animals, we surgically inserted two nylon tubes inside the tumor of the first mouse.The tubes contained an identically absorbing solution of black ink (Royal Talens, Ecoline 700) and water.The spectroscopically measured absorption coefficient of the solution was 0.75 mm −1 .For this experiment, we used a freshly sacrificed xenografted nude mouse with human breast cancer cell line MDA-MB-231, depicted in Fig. 6(a).We performed PA imaging of this mouse and measured the fluence distribution in a limited region of the mouse cross section using AO.Both PA and AO measurements were performed using 785 nm wavelength light.Figure 6(b) shows the reconstructed PA image of the mouse cross section in the proximity of the dotted line in Fig. 6(a).In the PA image, the two tubes (T 1 and T 2 ) have different image values.The image value from tube 1 is approximately 73% lower than that from tube 2, mainly due to differences in local fluence at the locations of the  We measured the fluence in a small ROI in this case because the azimuthal scanning is done by mechanically rotating the US array, which is very time-consuming.In the fluence map, the bright spot corresponding to the light illumination point appears to be approximately 2 mm under the skin.This is because the US focus has a finite size, and the maximum signal is observed when the US focus is completely inside the sample.After fluence correction, the peak image value of tube 2 is 10% higher than that of tube 1, which can be explained by a reconstruction artifact.Tube 1 is reconstructed nicely into a round shape, whereas tube 2 has an elongated shape, causing the peak image value to be lower.
We demonstrated the feasibility of mapping fluence variations in the entire cross section of the mouse, thereby enabling fluencecompensated PA imaging of a small animal.In this experiment we used a sacrificed xenografted nude mouse with human breast cancer cell line BT474.Figure 7(a) shows the two cross-sectional tomographic US images, one centered at the tumor and the second 9 mm towards the head, where kidneys and spleen are visible.The PA images are acquired using 12 projections, with steps of 30°.The PA images show vasculature inside the tumor and organs as well as the large blood vessels (aorta and inferior vena cava).In the cross section centered at the tumor, the aorta (A) and inferior vena cava (VC) are visible as two absorbing features, whereas in the cross section close to the kidneys, they are not very well resolved, which could be due to a reconstruction artifact.The measured fluence variation maps for both locations show six regions of high fluence under the skin corresponding to the illumination spots of the six optical fiber bundles.The light fluence decreases with the radial distance from the surface towards the center as well as with azimuthal distance away from the illumination spots.The rate at which fluence decreases is different in the two cross sections.There is more fluence in the center of the image when the cross section is through the tumor in comparison with the cross section through the spleen and kidneys.This difference could be due to the fact that the kidneys and the spleen are strongly absorbing organs, causing fluence to decay rapidly.Another difference is that in the cross section through the tumor, there is a large urine-filled bladder that has very weak scattering and absorption compared to tissue; as a result, more light can reach the center.However, the measured AO signal inside the bladder cannot be related to the local fluence, since due to very low scattering, the tagging efficiency η and the isotropy of the light must be affected.The effect of this can be observed in the results, as the signal from the illumination spot on the bladder  is low compared to other illumination spots.Another significant difference between the fluence maps of the two cross sections is that, in the cross section with the tumor, the fluence is more diffused from 60°to 150°.This can also be associated with the fact that there is more blood in the organs resulting in a higher absorption of light, which leads to a faster decay of fluence in the azimuthal direction.The comparison of Figs.For quantitative comparison, Fig. 7(e) plots the image values along the horizontal axis in between two solid lines in Figs.7(b) and 7(d).The image value of A/VC in the cross section through the tumor is 4 times higher than in the cross section through the kidneys.Assuming the amount of blood in A/VC at both locations is the same, the difference in image value should have been compensated by normalizing the PA images with fluence maps.However, other than fluence, the possible cause for differences in image value of A/VC could be the acoustic attenuation.In the cross section through the tumor, the PA signal from A/VC partly propagates through the large urine-filled bladder and partly through muscle and soft tissue; whereas in the cross section through the kidneys, the PA signal propagates partly through organs and intestines and partly through muscle and soft tissue.The acoustic attenuation by the urine-filled bladder is roughly 100 times less compared to soft tissue and muscle at 5 MHz.The kidneys have approximately 2 times higher attenuation than muscles at 5 MHz [44].As a result, the acoustic attenuation of the cross section through the kidneys is significantly higher, which results in relatively more absorption of the acoustic signal during outward propagation, leading to a lower image value of A/VC in this slice.We confirmed this hypothesis by reconstructing the PA images of both slices using only five projections, such that the first projection was at 0°and the fifth projection was at 150°.This way we excluded the projections in which the US propagated through the tissue with high attenuation difference in both slices.The results of this analysis showed that the PA image value of A/VC in the slice with the tumor decreased, whereas in the slice with the kidney, the value remained the same.This trend confirms our hypothesis that the large difference in image value of A/VC in both slices is due to the acoustic attenuation difference of tissue and not the error in fluence correction.Furthermore, by considering only five projections, the PA image value of the A/VC in the slice with the tumor is 14% higher than the value in the slice with the kidney.
In Fig. 7, the effect of fluence correction on the appearance of PA images is not prominent due to the large variation in PA image values.To highlight the effect of fluence correction on the appearance of PA image, Fig. 8 shows a comparison of the PA image before and after the fluence correction in a region marked by the dotted rectangle in Fig. 7(d).We considered this region because within this region there are no features with very large PA image value.In Fig. 8, features are marked with arrows to aid the visual comparison.
The implication of these results is that measuring fluence variations experimentally using AO in small animals is feasible.The acoustic attenuation may potentially affect the accuracy of fluence measurements, since ultrasonic tagging efficiency is pressure-dependent.However, this effect of acoustic attenuation of US tagging efficiency can be circumvented by performing AO measurements using US pressure high enough such that modulation process is saturated, as shown in Fig. 2(c).Further, while the fluence compensation in small-animal imaging improves the quantitative interpretation of the PA images and enhances the visibility of deeper lying contrast sources, for a complete quantitative picture, it is necessary to consider acoustic properties of the tissue in the image reconstruction process.

CONCLUSION
We have developed a method to measure fluence using ultrasonic tagging of light that, for the first time, is integrated directly into a PA tomography imaging system.Our implementation has general applicability for both tomographic and reflection mode PA setups.This new technology enables acquisition of two-dimensional PA images and fluence maps in a single instrument in a synergetic manner.The main advantage of our method is that it does not require a priori knowledge about optical properties of the medium, and it is applicable to optically and acoustically heterogeneous media.The method can be easily implemented in different PA imaging geometries.
The results presented here show the feasibility of this technology to provide correction of PA images for spatial as well as spectral variations in fluence.Moreover, the fluence correction of PA images using our method in phantom experiments establishes a linear relation between the PA image value and the absorption coefficients of chromophores.Based on the results of the phantom experiments, the accuracy of the fluence correction is within 8%.The results on estimation of blood oxygen saturation show that fluence correction of PA images using our method provides more accurate estimation compared to uncorrected PA images.Further, the results of fluence compensation in a mouse show the potential of the method to correct for differences in image value caused by the fluence variations.
Future research will focus on implementing our methodology in vivo.There are two challenges for the in vivo application of our method: one is the speckle decorrelation caused by the tissue dynamics and the second is the instrumental complexity because of using different types of lasers for AO and PA measurements.Speckle decorrelation is caused by the long measurement time needed to collect sufficient light to form a speckle image when using CW lasers.A previous study has demonstrated that the use of a pair of coherent pulses with pulse length of 5 ns and a temporal separation within one US cycle can be used to perform AO measurements in a fast decorrelating turbid medium [34].This is because a pulsed laser can deliver sufficient light in a single pulse to form a speckle image, thereby shortening the measurement time relative to the tissue decorrelation time.Such a laser then can be used to perform both PA and AO measurements, thereby reducing the instrumental complexity of our system.Hereafter, we will take into account the acoustic heterogeneities, either by experimentally measuring or using literature values assigned to identifiable organs in the US images during the PA image reconstruction and investigate the full advantage of fluence compensation in PA imaging.

Fig. 1 .
Fig. 1.Illustration of the experimental setup.AOM, acousto-optic modulator; PBS, polarization beam splitter; CCD, charged-coupled device (camera); rotation (360°); translation (200 mm).The inset marked with dotted rectangle shows the illustration of the method for AO scanning in transverse direction of US focus propagation.

Fig. 2 .
Fig. 2. Measured pressure distribution in the imaging plane of the US array; (b) Schlieren images of the US focus after 2.6 cm propagation through water (top), agar plus Intralipid phantom (middle) and chicken breast tissue (bottom); (c) measured change in contrast as a function of acoustic pressure; solid line is the quadratic fit, and squares are measured data.The error bars on measured data are standard deviations after 10 measurements.Trans.axis,transverse axis; P, pressure in MPa; a, constant; b, constant.

Fig. 4 .
Fig. 4. Comparison of absolute blood oxygen saturation (sO 2 ) estimation using PA without and with fluence correction in an optically heterogeneous phantom; (a) illustration of the phantom cross section with two blood-filled nylon tubes; (b) and (c) uncorrected PA images at 750 and 800 nm wavelengths; (d) measured fluence variations for 750 and 800 nm wavelengths along the lines marked by solid lines in (a), where the error bars represent the standard deviation of the 10 measurements; (e) and (f) corresponding fluence-corrected PA images; (g) estimated blood sO 2 using PA measurements before and after fluence correction plotted against the ground truth sO 2 measured using an oximeter (Avoximeter 1000E), where vertical error bars are based on the error propagated from fluence measurement (d) and horizontal errors bars are based on the specification of (Avoximeter 1000E).

Fig. 3 .
Fig. 3. Demonstration of 2D fluence compensation in PA images.(a) Photograph of the phantom with three identical absorbing inclusions at different locations; (b) PA image acquired using six projections; (c) measured fluence distribution using ultrasonically modulated light; (d) fluence-normalized PA image; (e) line profile through the absorbers (1, 2 and 3) before (broken line) and after (solid line) fluence correction.

Fig. 5 .
Fig. 5. Comparison of absolute blood oxygen saturation (sO 2 ) estimation using PA without and with fluence correction in porcine tissue sample.(a) Photograph of the sample; (b) US image of the sample with two blood-filled tubes embedded; (c) example of fluence-corrected PA image at 750 nm; (d) measured fluence variation along the radial line through the tubes, where the error bars represent the standard deviation of the five measurements; (e) estimated blood sO 2 using PA measurements before and after fluence correction plotted against the ground truth sO 2 measured using an oximeter (Avoximeter 1000E), where vertical error bars are based on the error propagated from fluence measurement (d) and horizontal errors bars are based on the specification of (Avoximeter 1000E).
tubes. Figure 6(c) shows the acousto-optically measured fluence distribution in the region marked by the dotted lines in Fig. 6(b).

Figure 6 (
d) shows the PA image normalized with the acousto-optically measured fluence map.It is apparent that after the normalization of the PA image with the fluence, the difference in image value of the two tubes containing identically absorbing black Ecoline ink (Royal Talens, 700) is accounted for.For quantitative comparison, Fig. 6(b) shows the line profiles of corrected and uncorrected PA images along the lines in Figs.6(b) and 6(d).

Fig. 6 .
Fig. 6.(a) Photograph of the xenografted nude mouse with two nylon tubes (T 1 and T 2 ) embedded in the tumor; (b) PA image of the mouse cross section in the proximity of the dotted line in (a); (c) measured fluence variation map in the region marked by dotted polygon in (b); (d) fluence-normalized PA image; (e) line profile of the uncorrected (b) and corrected (d) PA images along the line marked though the tubes.

Fig. 7 .
Fig. 7. Two-dimensional fluence compensation in a sacrificed mouse.(a) Cross-sectional US tomographic images at two locations 9 mm apart from each other; (b) uncorrected PA images; (c) acousto-optically measured fluence variation maps; (d) fluence-normalized PA images; (e) line profiles of the uncorrected and fluence-corrected PA images for quantitative comparison.Description of markers in (a) and (b): T, tumor; Sp, spine; B, bladder; A/VC, aorta and inferior vena cava (the two are not distinguishable in the images); K, kidney; and S, spleen.
7(b) and 7(d) qualitatively shows that the visibility of absorbing structures at larger depths in fluence-corrected PA images is improved.The image value of the feature A in both cross sections is increased in the fluence-corrected image compared to the PA image alone.After fluence correction, kidney and spleen appear to have more uniform image values throughout the entire organ.

Fig. 8 .
Fig. 8. Zoomed in on a region marked by the dotted rectangle in Fig. 7(d) for the visual comparison of the fluence compensation effect on PA image.(a) Before fluence compensation and (b) after fluence compensation.