Quantitative photoacoustic tomography

In this paper, several algorithms that allow for quantitative photoacoustic reconstruction of tissue optical, acoustic and physiological properties are described in a finite-element method based framework. These quantitative reconstruction algorithms are compared, and the merits and limitations associated with these methods are discussed. In addition, a multispectral approach is presented for concurrent reconstructions of multiple parameters including deoxyhaemoglobin, oxyhaemoglobin and water concentrations as well as acoustic speed. Simulation and in vivo experiments are used to demonstrate the effectiveness of the reconstruction algorithms presented.


Introduction
Biomedical photoacoustic tomography (PAT) is a potentially powerful modality that can offer high-resolution structural and functional imaging of tissue (Kruger & Liu 1994;Kruger et al. 1999;Oraevsky et al. 1999;Viator et al. 2002;Norton & Vo-Dinh 2003). In PAT, a short-pulse laser source is used to irradiate the tissue of interest. The laser-produced temperature rise and subsequent thermoelastic expansion of tissues generate an acoustic wave that is detected by ultrasound transducers along multiple boundary positions. A reconstruction algorithm is used to recover the photoacoustic (PA) images quantitatively or qualitatively. PAT has shown the potential to detect breast cancer, to probe brain functioning in small animals, and to assess vascular and skin diseases (Hoelen et al. 1998;Wang et al. 2002;Yin et al. 2004;Niederhauser et al. 2005;Manohar et al. 2007).
However, conventional PAT reconstruction methods are mostly qualitatively based, and can image only the distribution of absorbed light energy density, which is the product of both the local optical absorption coefficient and the optical fluence distribution within the irradiated medium. It is well known that it is the tissue absorption coefficient that directly correlates with tissue structural and functional information such as haemoglobin and blood oxygenation. Several recent studies have suggested that it is possible to recover the quantitative optical absorption coefficient map when conventional PAT is combined with a light transport model (Ripoll & Ntziachristos 2005;Cox et al. 2006;Yin et al. 2007;). These studies have opened a new avenue to realize truly quantitative PAT by exploiting the spectral characteristics of specific chromophores in tissue, thus providing spatially resolved quantitative physiological and molecular information for valuable diagnostic decision making.
In addition to the optical properties and physiologically relevant tissue parameters, including deoxyhaemoglobin (Hb), oxyhaemoglobin (HbO 2 ) and water (H 2 O), PAT can also capture tissue mechanical parameters such as acoustic velocity. Our recent studies in quantitative PAT have shown that the images of both acoustic velocity and absorbed light energy density can be recovered simultaneously . This ability to recover both optical and acoustic properties not only provides more accurate reconstruction of optical properties than for qualitative PAT because of the elimination of the assumption of homogeneous acoustic velocity built into the qualitative PAT methods, but also adds the potential to better differentiate benign from malignant lesions as it is known that there exist significant differences in acoustic properties between normal and tumour tissues (Greenleaf & Bahn 1981).
In this paper, we will demonstrate and compare our quantitative PAT approaches that are based on the finite-element method. We will discuss the main characteristics of the reconstruction algorithms that allow for quantitative recovery of the optical absorption coefficient as well as the simultaneous reconstruction of the acoustic velocity and absorbed energy density.
In addition, advances in diffuse optical tomography (DOT) have shown that absorbers can be recovered with high accuracy when a priori spectral information provided by near-infrared spectroscopy is incorporated into the DOT reconstruction (Corlu et al. 2003(Corlu et al. , 2005Li et al. 2004). Similarly, we are able to develop a new spectral approach that allows for direct simultaneous reconstruction of tissue chromophores and acoustic velocity using multiplewavelength PA measurements. Owing to the use of multiple laser wavelengths, this approach provides more accurate parameter reconstruction, especially for acoustic velocity, as well as direct recovery of physiological parameters with enhanced accuracy over the Beer's law-based fitting methods (Laufer et al. 2007). We will demonstrate this multispectral PAT approach using simulations and in vivo experiments.

Overview of finite-element-based quantitative PAT reconstruction methods
(a ) Simultaneous recovery of acoustic velocity and absorbed energy density Since our reconstruction algorithm for the simultaneous reconstruction of acoustic property and absorbed light energy density has been described in detail elsewhere , we give only a brief outline here for context. The frequency-domain Helmholtz wave equation in an acoustically heterogeneous medium can be stated as V 2 pðr; uÞ C k 2 0 ð1 C OÞpðr; uÞ Z ik 0 v 0 bJðrÞ C p : ð2:1Þ Here p is the acoustic pressure; k 0 Zu/v 0 is the wavenumber, described by the angular frequency, u, and the speed of the acoustic wave in a reference or coupling medium, v 0 ; b is the thermal expansion coefficient; C p is the specific heat; J is the absorbed light energy density, which is the product of the optical absorption coefficient and incident optical fluence; and O is a coefficient that depends on both acoustic speed and attenuation as in which v is the speed of the acoustic wave in the medium/tissue; and a is the acoustic attenuation coefficient. The finite-element discretization of equation (2.1) is stated as Ap Z bJ; ð2:3Þ where the elements of the matrix A and b are expressed as Here J, O and p are expressed by their real and imaginary components; j i is the basis function; h i indicates integration over the problem domain; and the second-order absorption boundary conditions are employed, where hZ ðKik K3=2rC i3=8kr 2 Þ=ð1K i=krÞ; gZ ðKi=2kr 2 Þ=ð1K i=krÞ; and q is the angular coordinate at radial position r. To obtain a matrix equation capable of inverse solution, a combined Marquardt and Tikhonov regularization scheme is used for the inverse calculation, .; p c M Þ T ; p o i and p c i are observed and computed complex acoustic field data for iZ1, 2, ., M boundary measurement locations; Dc is the update vector for the optical and acoustic properties; J is the Jacobian matrix formed at the boundary measurement sites; l 0 is a scalar; and I is the identity matrix. Thus, here the image formation task is to update the absorbed light energy density and acoustic property distributions via iterative solution of equations (2.3) and (2.6), so that a weighted sum of the squared difference between the computed and measured acoustic pressure can be minimized.
For most cases where the heterogeneity/object has contrast in both the absorbed energy density and acoustic velocity, the energy density and acoustic velocity can be recovered at a single wavelength. This is very similar to what is seen in DOT, where both absorption and scattering coefficients can be reconstructed at a single wavelength when the objects have contrast in both the absorption and scattering. For such cases, crosstalk certainly exists, but it is not strong. This crosstalk is also dependent on the level of contrast for the two parameters. For the numerical and experimental cases we have previously investigated in PAT , the level of contrast used was relatively low, which explains the insignificant crosstalk that occurred in our reconstruction.
However, for the cases where acoustic heterogeneity is quite different from the distribution of absorbed light energy density, our simulation tests indicate that the crosstalk cannot be minimized using one wavelength measurement. For our simulation tests, the circular background region (15 mm in radius) contains three circular targets (2 mm in radius each) positioned at 2 (top right), 9 (top left) and 6 o'clock (bottom). The top left inclusion has contrast for velocity, while the bottom and top right inclusions have contrast for absorbed energy density. The velocities for the targets and background are 1.54!10 6 and 1.44!10 6 mm s K1 , respectively, while the absorbed energy densities for the background and targets are 1 and 10 mJ mm K3 , respectively. We can see from the recovered images shown in figure 1 that the absorbed energy density distribution has some crosstalk effect on acoustic velocity, while the acoustic velocity distribution has no significant effect on the recovered absorbed energy density due to its low contrast.

(b ) Recovery of absolute absorption coefficient
We have developed two powerful algorithms (Yin et al. 2007;) that are able to recover the optical absorption coefficient from absorbed light energy density obtained by the conventional PAT method. These two algorithms are both based on the finite-element solution to the PA wave equation coupled with the photon diffusion equation.
(i) Algorithm 1 This algorithm includes two steps (Yin et al. 2007). The first is to obtain the map of absorbed light energy density through a model-based reconstruction algorithm that is based on the finite-element solution to the PA wave equation in the frequency domain. The second step is to obtain the distribution of optical fluence using a photon diffusion equation-based optimization procedure and to recover the distribution of the optical absorption coefficient from the optical fluence and the absorbed energy density obtained from the first step. In this method, the optical fluence is obtained using diffuse optical measurements along the boundary of the tissue/medium and through the iterative calculation of the following equations (Iftimia & Jiang 2000): V$DðrÞVFðrÞK m a ðrÞFðrÞ ZKSðrÞ; KDVF$n Z aF; ð2:7Þ are the measured and calculated optical fluence for iZ1, 2, ., M boundary locations; m a is the optical absorption coefficient; and D is the optical diffusion coefficient. The goal of the second step is to obtain the distribution of the optical fluence within the entire imaging domain through an optimization procedure based on equations (2.7) and (2.8), which then allows the separation of the optical absorption coefficient from the optical fluence using the absorbed energy density obtained from the first step.
While algorithm 1 is very effective in recovering the absorption coefficient from the PAT data, there are some limitations associated with this approach. First, it requires additional diffuse light measurements, making the implementation of the hardware more complex. Second, the recovered results may depend on the accuracy of the distribution of the absolute absorbed energy density from conventional PAT.
(ii) Algorithm 2 To overcome these limitations mentioned above, a novel reconstruction approach that combines conventional PAT with the diffusion equation-based regularized Newton method for accurate recovery of optical properties has been proposed ). This approach is based on the following photon diffusion equation as well as the Robin boundary conditions, in terms of the absorbed energy density, JZm a F:

V$DðrÞVðEðrÞJðrÞÞKJðrÞ ZKSðrÞ;
ð2:9Þ KDVðEðrÞJÞ$n Z EðrÞaJ: ð2:10Þ Here E(r)Z1/m a (r); S(r) is the incident point or distributed source term; and D is the optical diffusion coefficient. For the inverse computation, the so-called Tikhonov regularization sets up a weighted term as well as a penalty term in order to minimize the squared differences between computed and measured absorbed energy density values, where L is the regularization matrix or filter matrix (Yuan et al. 2008 i is the absorbed energy density obtained from PAT; and J c i is the absorbed energy density computed from equations (2.9) and (2.10) for iZ1, 2, ., N locations within the entire PAT reconstruction domain. The initial estimate of the absorption coefficient can be updated based on the iterative Newton method as follows with bZ1: The basic idea of this approach is to incorporate the high-resolution PAT images into the DOT-like reconstruction, so that both the resolution and quantitative accuracy of optical image reconstruction are enhanced. However, the accuracy of the recovered absorption coefficient depends on the structural information obtained from PAT. Fortunately, our regularization-based diffuse optical reconstruction algorithm is immune to imperfect a priori spatial information (Yuan et al. 2008).
Further, it is noted that the scattering coefficient is assumed as constant/ homogeneous for the use of both algorithms 1 and 2. While reconstructing only the absorption coefficient minimizes the non-uniqueness issue, it certainly imposes a limitation on our algorithms. This suggests that our methods may be applicable to cases where the scattering contrast is low or the objects having scattering contrast are small in size; otherwise, we may have to use DOT as we have indicated in our published work (Yin et al. 2007).

Development of spectrally resolved PAT
Early work in PAT focused on reconstructing tissue optical properties at several selected wavelengths (Laufer et al. 2007). A least-squares fitting algorithm was then used to estimate chromophore concentrations based on the recovered optical properties and Beer's law. However, our simulation tests indicate that the water content image cannot be effectively recovered using a fitting method. In our study, multispectral PAT is proposed to directly reconstruct chromophore concentrations in tissue without the use of a fitting method.
In multispectral PAT, the frequency-domain Helmholtz wave equation in an acoustically heterogeneous medium is written as V 2 pðr; u; lÞ C k 2 0 ð1 C OÞpðr; u; lÞ Z ik 0 v 0 bm a ðr; lÞFðr; lÞ C p ; ð3:1Þ where l is the wavelength of the incident light. In addition, according to Beer's law, the wavelength-dependent tissue absorption can be expressed as where c i is the concentration and 3 i (l) is the extinction coefficient of the ith chromophore (HbO 2 , Hb or H 2 O) at wavelength l. In light of equation (3.2), equation (3.1) can be rewritten as Thus, similar to equation (2.6), the inverse solution can be obtained by solving the equation ð3:4Þ in which DcZð6O 1 ;6O 2 ;.;6O n ;6c 1;1 ;6c 1;2 ;.;6c 1;n ;6c 2;1 ;6c 2;2 ;.;6c 2;n ; 6c 3;1 ;6c 3;2 ;.;6c 3;n Þ T is the update vector for chromophores and acoustic velocity; l 0 is the regularization parameter determined by combined Marquardt and Tikhonov regularization schemes; and p o i and p c i are measured and calculated data for iZ1, 2, ., M boundary locations. For each acoustic frequency u within each incident optical wavelength l, these are written as The Jacobian matrix, J, is denoted as J Z ½J O;l;u ;J c i ;l;u , whereJ O;l;u andJ c i ;l;u represent the Jacobian submatrix for acoustic velocity and different chromophores, respectively. It should be noted that, in this reconstruction algorithm, we assume that the optical fluence F can be estimated through a photon diffusion equation during the image reconstruction procedure .

Results and discussion
In this section, we show reconstruction results that demonstrate the feasibility of the multispectral PAT described above. The multispectral PAT approach is evaluated using several simulations and small animal experiments.
For the simulations, we use a circular background (15 mm in radius) containing three circular targets (2 mm in radius each), where each inclusion had different contrast for each given parameter, namely Hb, HbO 2 , H 2 O and acoustic velocity. The chromophore concentrations and velocity used for all the test cases are listed in table 1. Six optical wavelengths (633, 682, 723, 850, 854 and 930 nm) were chosen for the spectrally resolved PAT. The 'measured' data were generated using a forward computation with 2 per cent random Gaussian noise. A total of 120 receivers were equally distributed along the boundary of the circular background. Fifty frequencies ranging from 50 to 540 kHz were used at each optical wavelength for each receiver. The extinction coefficient for each chromophore used was taken from S. Pahl show the true locations of the targets including Hb, HbO 2 , H 2 O and acoustic velocity. The reconstructed images shown in figure 2b(i)-(iv) correspond to Hb, HbO 2 , H 2 O and acoustic velocity. We see from figure 2 that the crosstalk errors between the chromophore concentrations and acoustic velocity have been effectively reduced using our quantitative PAT approach.
Animal experiments were performed using a PAT system that has been described in detail elsewhere ). In this system (figure 3), pulsed light from a Ti:sapphire laser (Symphotic Tii Corp., Camarillo, CA) is  used to irradiate the tissue and generates an acoustic pressure wave. A 1 MHz transducer (GE Panametrics, Waltham, MA) is used to receive the acoustic signals. The transducer and the animal (oxygen provided) are immersed in a water tank. A rotary stage rotates the receiver relative to the centre of the tank. One set of data is taken at 120/360 positions when the receiver is scanned circularly over 3608. The complex wavefield signal is first amplified by a preamplifier, and then amplified further by a pulser/receiver. A data acquisition board converts it into a digital signal, which is fed to a computer. The radius of the receiver motion path was 100 mm in the animal experiments conducted. The entire data acquisition is realized through C programming. In the current version of the system, data collection for a total of 120 measurements requires approximately 2 min. Figure 4 shows a photograph of a typical mouse with an implanted subcutaneous tumour. Five optical wavelengths (755, 800, 860, 900 and 930 nm) were used for the measurements. Figure 5a-d presents the reconstructed in vivo Hb, HbO 2 , H 2 O and acoustic velocity images. We see that the tumour is remarkably imaged with the highest contrast in Hb and HbO 2 . In addition, the H 2 O image is also effectively recovered with better quality compared with that from diffuse optical imaging (Corlu et al. 2005). The contrast in velocity between the tumour and surroundings (figure 5d ) is apparently lower than that for the physiological parameter images (figure 5a-c). The in vivo results shown here further validate the merits of the multispectral PAT method. For all the tests, multispectral PAT can effectively minimize the crosstalk and improve the image quality. The direct reconstruction method has demonstrated that improved quantification of chromophore concentrations is achieved with multi-wavelength measurements when spectral priors are incorporated into the inverse reconstruction.

Conclusions
In this paper, we have described and compared several algorithms for realizing quantitative PAT. In addition, we have demonstrated experimental evidence that both physiological and acoustic property images can be obtained using our multispectral PAT approaches. The multispectral PAT approach presented provides an efficient means of concurrent reconstruction of multiple parameters including different chromophore concentrations and acoustic velocity. Images generated from both noisy simulated data and measured in vivo data have been provided. It is clear from the reconstructed results that the crosstalk issue for recovery of multiple parameters can be reduced or minimized using multispectral PAT approaches. We also found that the selection of the number of wavelengths (Corlu et al. 2003(Corlu et al. , 2005 as well as the wavelength range is critical in obtaining reliable and accurate reconstruction of multiple parameters, and the findings will be discussed in detail elsewhere (Yuan & Jiang 2009). Nonetheless, the ability to reconstruct multiple-parameter images may provide us with a convenient tool to quantify physiological function, disease progression or response to intervention.