## Abstract

In biomedical photoacoustic imaging the images are proportional to the absorbed optical energy density, and not the optical absorption, which makes it difficult to obtain a quantitatively accurate image showing the concentration of a particular absorbing chromophore from photoacoustic measurements alone. Here it is shown that the spatially varying concentration of a chromophore whose absorption becomes zero above a threshold light fluence can be estimated from photoacoustic images obtained at increasing illumination strengths. This technique provides an alternative to model-based multiwavelength approaches to quantitative photoacoustic imaging, and a new approach to photoacoustic molecular and functional imaging.

©2010 Optical Society of America

## 1. Introduction

Photoacoustic tomography (PAT) can provide images related to the
distribution of optical absorption within highly scattering media such as
biological tissue. Its many applications to date include, among other
things, preclinical studies of tumor growth, breast imaging, imaging of
the microvasculature and small animal imaging [1, 2, 3, 4, 5]. Photoacoustic
images, while useful, and while closely related to the distribution of
optical absorption, do not give quantitative measures of absorption as
they are weighted by the unknown, non-uniform, light fluence in the
tissue. Rather, the image amplitude for any given point is proportional to
the absorbed optical energy density there, *h* =
*μ _{a}ϕ*, where

*ϕ*is the local light fluence and

*μ*the optical absorption coefficient. A serious implication is that multiwavelength spectroscopic methods cannot be applied directly to PAT images, which significantly restricts its use to only a qualitative imaging modality for many applications. If this dependence on the light fluence could be removed then PAT would quickly become a very powerful tool for molecular and functional imaging, particularly for preclinical studies.

_{a}Attempts to estimate the absorption coefficient distribution, or the concentration distribution of a chromophore, have so far mostly been model-based, ie. the optical coefficients in a model of diffuse light transport were iteratively adjusted until the modelled image matched that measured. [6, 7, 8, 9, 10] One complication of the model-based approach is the non-uniqueness introduced if the optical scattering is also unknown and spatially-varying, although this may be overcome either by using multiple illumination patterns, [11], or prior knowledge of the spatial geometry of the absorbers [10, 12], by assuming knowledge of the wavelength dependence of the scattering and the specific absorption spectra of the chromophores. [8, 9]

Another potential drawback to model-based methods is the size of the
inverse problem due to the large image size/spatial resolution ratio
available in PAT, particularly a problem when using Hessian-based
inversion methods. [13] A different
approach, using complex geometrical optics, has also been proposed, [11] but this shares a requirement of
the model-based methods that *μ _{a}* > 0 everywhere (and that in practice the image amplitude is
above the noise everywhere). This is necessary for these methods because
otherwise the images (the data) is not carrying sufficient information
about the fluence.

It would be much simpler, and would avoid all these problems, if it were
possible to measure the fluence within the tissue directly. Unfortunately
there are currently no modalities that can do this with the same depth and
resolution as PAT. However, it may be possible *indirectly*
to measure the fluence at the points of interest (ie. where the PAT image
is non-zero). Here we suggest the use of a contrast agent whose absorption
falls quickly to zero above some known fluence. In this case it may be
possible to determine the fluence inside the tissue with the same
resolution as a PAT image. We show how this may allow the spatially
varying concentration distribution of the chromophore to be recovered.

## 2. Fluence Dependent Chromophores

The specific absorption coefficient of many absorbers varies, often
irreversibly, with the light fluence, *α* =
*α*(*ϕ*). For example, many organic dyes
show fluence-dependent changes in their absorption due to optical
saturation [14], transient
formation of isomers [15], or even
permanent photobleaching [16]. Gold
nanorods undergo shape changes above a certain fluence level, which result
in significant changes to the absorption spectrum in the near-infrared.
[17, 18] Indeed, at a given wavelength, the absorption of
a gold nanorod can decrease abruptly as soon as the fluence reaches a
threshold related to its melting point.

The question of interest here is: are there any types of fluence dependence
(any functions *α*(*ϕ*)) which would allow
the concentration of the chromophore with this dependence to be recovered
from photoacoustic images? The idea would be to obtain images at multiple
illumination strengths and from these, and knowledge of
*α*(*ϕ*), to recover the chromophore’s
concentration at every point. This could be considered analogous to
estimating concentration using images at multiple
*wavelengths* due to the wavelength dependence of
absorption *α* = *α*(*λ*),
except that - as was mentioned above - such spectroscopy is not possible
with PAT images without knowledge of the fluence. So what is different
here? The difference is that the absorption of a fluence dependent
chromophore could, indirectly, give information about the values of the
fluence inside the tissue. This extra information cannot be accessed by
varying the wavelength and would be sufficient to allow the concentration
to be estimated.

## 3. Estimating Spatially Varying Chromophore Concentrations

Here, one type of fluence dependence is considered: a chromophore whose
specific absorption coefficient falls rapidly to zero above some threshold
fluence. The specific absorption coefficient for this ‘switching-off type
of absorber can be idealized as *α*(*x*) =
*α*
_{0}
*U*(*ϕ _{th}* -

*ϕ*(

*x*)) where

*U*(

*y*) is the unit step function (0 for

*y*< 0, 1 for

*y*≥ 0) and

*α*

_{0}is the specific absorption coefficient at fluences below the known threshold fluence

*ϕ*. (Gold nanorods that undergo a rapid shape change above a certain fluence, as mentioned above, are potential candidates.) Consider an example consisting of a distribution of a chromophore with this switching property in a turbid medium with background absorption

_{th}*μ*

_{a0}. (It is not necessary to assume that this is spatially uniform.) Here it is shown that the chromophore concentration distribution,

*C*(

*x*), can be estimated using a series of PAT images generated with gradually increasing illumination strengths. (The incident light intensity is increased, in small steps, for each subsequent image.) In general, the total absorption at a point

*x*in the tissue may be written as the combination

*μ*(

_{a}*x*) =

*μ*

_{a0}(

*x*)+

*C*(

*x*)

*α*(

*ϕ*(

*x*)) where the second term is a contribution from the fluence dependent chromophore, and

*C*(

*x*) is its concentration. Substituting the expression for

*α*from above gives the absorption coefficient as

When the fluence everywhere is below the threshold *ϕ _{th}*, the magnitude of the PAT image

*h*(

*x*) will increase linearly with the illumination intensity. Just before the threshold fluence is reached somewhere, the image amplitude at a point

*x*can be written as

(Here we will assume that the absorbed energy images can be reconstructed
exactly. In practice, artefacts and uncertainties will, of course, affect
the accuracy with which the chromophore concentrations can be estimated.)
If the illumination intensity is increased by a small factor
*k* so that the fluence reaches the threshold in some
places, then in those regions the chromophore will become non-absorbing
(will ‘switch-off’).

It is helpful at this stage to divide the domain into three non-overlapping
regions, *A _{n}*,

*B*, and

_{n}*C*, which change with each increase in fluence, such that on step

_{n}*n*:

*A*contains the points where the fluence is below the threshold and always has been,

_{n}*A*= {

_{n}*x*:

*ϕ*(

_{m}*x*) <

*ϕ*for all

_{th}*m*≤

*n*} (chromophore still absorbing);

*B*contains the points where the fluence exceeded the threshold for the first time on this step,

_{n}*B*= {

_{n}*x*:

*ϕ*(

_{m}*x*) <

*ϕ*for all

_{th}*m*<

*n*- 1 and

*ϕ*(

_{n}*x*) ≥

*ϕ*} (chromophore just switched off);

_{th}*C*contains the points where the fluence exceeded the threshold on some previous step, irrespective of its current value,

_{n}*C*= {

_{n}*x*:

*ϕ*(

_{m}*x*) ≥

*ϕ*for at least one

_{th}*m*≤

*n*} (chromophore switched off on a previous step). In practice there may be a delay of one step between reaching the threshold and the chromophore not absorbing, but this does not affect the argument significantly.

The image at an arbitrary step, *n* (the fluence has been
increased *n* times) can now be written as

When *Cα*
_{0}, the absorption change in *B _{n}*, is small enough that the fluence at regions outside

*B*will not be changed significantly by the chromophore switching off, the fluence can be approximated by

_{n}In other words, within the regions *A _{n}* and

*C*the fluence distribution,

_{n}*ϕ*, can be approximated by scaling the previous fluence

_{n}*ϕ*

_{n-1}, but in region

*B*, where the absorption has changed, there will also be a perturbation to the fluence. The difference between an image

_{n}*h*and a scaled version of the previous image

_{n}*kh*

_{n-1}can now be written as

The term *μ*
_{a0}
*δϕ*
_{n} was disregarded because *δϕ*
_{n} is small, by design, as the illumination strength is only increased
slightly. (*μ*
_{a0} may also be small.) By making the
approximation *kϕ*
_{n-1}(*x*) ≈
*ϕ _{th}*,

*x*∈

*B*, an expression for the concentration of the chromophore in the region

_{n}*B*can be obtained:

_{n}Initially, any point *x* is contained in region
*A*
_{0}, so *B*
_{0} = *C*
_{0} = {}. As the incident illumination is increased (and assuming
the optical properties take physiologically realistic values),
*x* will go from being in *A _{n}* to

*B*and eventually to

_{n}*C*. In other words, region

_{n}*A*will shrink,

_{n}*A*=

_{n}*A*

_{n-1}-

*B*, and region

_{n}*C*will grow,

_{n}*C*=

_{n}*C*

_{n-1}+

*B*

_{n-1}until all

*x*lie in

*C*for some sufficiently large

_{n}*n*. Because of this, the set

*B*

_{1}∪

*B*

_{2}∪… will eventually contain every point in the domain once. This, and the fact that the concentration estimate above is zero outside

*B*allows the chromophore concentration to be obtained everywhere by summing the estimates obtained at each step:

_{n}## 4. Numerical Simulation

In order to investigate this method of recovering chromophore
concentration, a finite element model of the diffusion approximation to
the light transport equation was used to simulate the above proposal.
[19] A point source was located
centrally at 0.1 mm depth in a 15mm × 15mm domain (represented by a
100×100 pixel mesh), and the background scattering was set by defining the
optical diffusion constant to be *D* = 1.5 mm. A boundary
condition was imposed to ensure there was no incoming photon flux from
outside the domain. The anisotropy factor was set to *g* =
0.9. Gaussian random noise was added to the image of absorbed energy
density at 5% of the mean value. At each step the source strength was
increased by 10% (*k* = 1.1). The amplitude of the fluence
was set initially so that it was everywhere below the value of the
threshold fluence *ϕ _{th}*. The background absorption coefficient was

*μ*

_{a0}= 0.02 mm

^{-1}. The nonlinear absorbers were restricted to three regions: two squares, in which the absorption coefficient (before switching off) was

*μ*

_{a,sq}=

*C*

*α*

_{0}= 0.04 mm

^{-1}, and a circle where

*μ*

_{a,circ}=

*Cα*

_{0}= 0.02 mm

^{-1}. In fact, for convenience and without loss of generality,

*α*

_{0}was set to 0.04 (squares) and 0.02 (circle) so that the concentration distribution

*C*(

*x*) took values from 0 to 1 (see Fig. 2A). (Deliberately, no specific units of concentration have been indicated, to avoid the erroneous idea that these numbers somehow indicate a minimum detectable concentration. The minimum concentration that could be recovered in practice will depend on both the molar absorption coefficient of the actual contrast agent used, and the SNR, which itself will be a factor of the instrumentation.)

Figure 1 shows the PAT image
amplitude at points (3.75mm,3.75mm), dashed line, and (7.5 mm,7.5mm),
solid line. In the former case the absorption coefficient remains the same
throughout the experiment, so the image amplitude just increases by 10% at
each step as *k* = 1.1. At the latter point - located at
the centre of the domain - the nonlinear chromophore is present, and at
step 22 the fluence there has reached the threshold fluence,
*ϕ _{th}*, thus rendering the chromophore non-absorbing from then on. The
abrupt jump in the photoacoustic image amplitude that results can be seen
clearly in the solid curve of Fig.
1. As the threshold fluence is known (it is a property of the
chromophore) it is features of this sort that give an indication of the
fluence within the tissue.

Figure 2A shows the true
concentration of the nonlinear chromophore, and Fig. 2B is the estimate obtained using Eq. (4). Figure 3 shows a profile through both Figs. 2A (dashed line) and 2B (solid
line). There is clearly good agreement between the estimated concentration
and the true values. This will improve further with increases in the
signal-to-noise ratio, but decrease near the edges of objects with a large
value of *Cα*
_{0}, as Eq. (2)
becomes increasingly approximate.

## 5. Discussion

#### 5.1. Residual absorption

The model of the contrast agent used here - one in which the absorption
falls abruptly to zero at a threshold fluence - is clearly an
idealisation. In practice there is likely to be a remaining residual
absorption, due in part to the probabilistic nature of molecular
changes. The details concerning this residual absorption, and its
dependence on parameters such as the illumination duration, will
depend on the mechanism by which absorption changes occur; if the
absorption change is due to some kind of photochemical effect then the
situation is likely to be different from when the absorption change is
due to a thermal mechanism. If the residual molar absorption
coefficient, *α*
_{res} (which could be determined experimentally beforehand
for a given pulse energy and duration) is not affected by subsequent
illuminations, then it can be incorporated straightforwardly into the
scheme of Section 3 by redefining the background absorption to include
it: *μ*
_{a0}
^{′}(*x*) = *μ*
_{a0}(*x*) +
*C*(*x*)*α*
_{res} and *α*
_{0} in Eq. (3)
is replaced by *α*
_{0} - *α*
_{res}. It is more likely, though that the residual absorption
will decrease due to the light from subsequent steps, which would
complicate the concentration estimation, and it would be pragmatic to
look for measures to minimise the residual absorption should be taken.
One way to do this might be to optimise the combination of pulse
energy and illumination pulse duration.

Consider the case of a heat-induced shape change in a gold nanorod. (a)
*Fixed pulse energy*: for any illuminating pulse
whose duration is short enough that thermal diffusion is negligible
then, for a fixed pulse energy, the same degree of heating will occur
in the nanorods leading to the same absorption change. However, the
thermal diffusion time will depend on a number of factors (the tissue
type, the density of the gold nanorods within the tissue [20], whether they are clumped,
any coating they are carrying [21], etc) and for a single nanorod may well be less than
the laser pulse duration, in which case thermal diffusion may become
important and shorter pulse durations could lead to larger absorption
changes for a given pulse energy. (b) *Fixed peak
power*: if the fluence rate is high enough to deliver
enough energy to overcome thermal diffusion and cause sufficient
temperature rise to cause the require shape changes, then increasing
the pulse duration while keeping the peak power constant could reduce
the residual absorption simply by virtue of supplying more energy. One
potential drawback of a longer pulse duration is that it might result
in a longer scan, and therefore an increased chance that the
background absorption will change during the measurements.

#### 5.2. Fluence dependence in real contrast agents

In general, the molar absorption coefficient *α* of a
contrast agent is likely to depend on the fluence in a more
complicated way than has been assumed here. As a simple example: if
the nanorod population had a range of sizes then there may be a
transitional range of fluences from *ϕ _{L}* to

*ϕ*, over which

_{H}*α*decreases as the fluence increases, rather than the abrupt change assumed above. If the molar absorption as a function of fluence does not contain clear features (abrupt changes or changes of gradient, for instance) associated with known fluences - in the way that the cut-off was associated with the threshold fluence

*ϕ*- then recovering the chromophore concentration using a variant of the approach described here becomes more difficult. Where such features are present though, it may be possible to devise a means by which to estimate the concentrations.

_{th}#### 5.3. Fluence increment and SNR

There is a trade-off at the heart of the method, related to the size of
the fluence increment. A larger increment could result in a greater
signal-to-noise ratio (SNR) in the difference image but may violate
the assumption that *δ _{ϕ}* is small. Also, when the fluence level is low, in the
initial steps, the SNR is likely to be worse than at later steps where
the fluence is greater. The optimal increment may therefore change
with the fluence level, and it may be advantageous to make larger
increases in fluence at the low levels than at the higher levels later
on. The multiplicative increase factor

*k*could be adjusted at each step, and these

*k*taken into account when estimating the fluence.

_{n}## 6. Summary

This paper has introduced the idea of using the fluence dependence of the
specific absorption coefficient of a chromophore to estimate its
concentration. It has been shown that the concentration of a chromophore
may be recovered from photoacoustic images obtained at increasing
illumination strengths if the chromophore stops absorbing above a critical
fluence threshold. It is also known that an equivalent result holds for a
chromophore that switches *on* at a fluence threshold.
[22] The following general
questions remain open: what characteristics must the fluence dependent
specific absorption coefficient *α*(*ϕ*)
have if it is to enable chromophore concentration estimation? Is it
possible to recover the concentrations of two or more chromophores with
different fluence dependences? These results could lead to the design of
new contrast agents that would facilitate single wavelength quantitative
(and thus molecular) photoacoustic imaging.

## Acknowledgements

The authors would like to thank the anonymous reviewers for their insightful comments. This work was supported by the Engineering and Physical Sciences Research Council, UK.

## References and links

**1. **R.A. Kruger, K.D. Miller, H.E. Reynolds, W.L. Kiser, D.R. Reinecke, and G.A. Kruger, “Contrast enhancement of breast cancer
in vivo using thermoacoustic CT at 434 MHz - feasibility
study,” Radiology **216**, 279–283
(2000)

**2. **X. Wang, Y. Pang, G. Ku, X. Xie, G. Stocia, and L. V. Wang, “Noninvasive laser-induced
photoacoustic tomography for structural and functional in vivo
imaging of the brain”, Nature
Biotech. **21**(7), 803–806
(2003)

**3. **E.Z. Zhang, J.G. Laufer, R.B. Pedley, and P.C. Beard, “In vivo high-resolution 3D
photoacoustic imaging of superficial vascular
anatomy,” Phys. Med. Biol. **54**, 1035–1046
(2009)

**4. **L.V. Wang, ed., *Photoacoustic Imaging and
Spectroscopy*, CRC Press,
2009.

**5. **A.A. Oraevsky and L.V. Wang, eds., *Photons Plus Ultrasound: Imaging and
Sensing*, Proc. SPIE **7564** (2010)

**6. **B.T. Cox, S. R. Arridge, K. Köstli, and P.C. Beard, “Two-dimensional quantitative
photoacoustic image reconstruction of absorption distributions in
scattering media by use of a simple iterative
method,” Appl. Opt. **45**, 1866–1875
(2006)

**7. **H. Jiang, Z. Yuan, and X. Gu, “Spatially varying optical and
acoustic property reconstruction using finite-element-based
photoacoustic tomography,” J. Opt. Soc.
Am. A **23**(4), 878–888
(2006)

**8. **B.T. Cox, S. R. Arridge, and P.C. Beard, “Estimating chromophore distributions
from multiwavelength photoacoustic images,”
J. Opt. Soc. Am. A , **26**,
443–455
(2009)

**9. **J.G. Laufer, B. T. Cox, E.Z. Zhang, and P.C. Beard, “Quantitative determination of
chromophore concentrations from 2D photoacoustic images using a
nonlinear model-based inversion scheme,”
Appl. Opt. **49**, 1219–1233
(2010)

**10. **L. Yao, Y. Sun, and H. Jiang, “Transport-based quantitative
photoacoustic tomography: simulations and
experiments,” Phys. Med. Biol. **55**, 1917–1934
(2010)

**11. **G. Bal and G. Uhlmann, “Inverse diffusion theory of
photoacoustics,” arXiv: 0910.2503v0911 [math.AP]
(2009)

**12. **A. Rosenthal, D. Razansky, and V. Ntziachristos, “Quantitative Optoacoustic Signal
Extraction Using Sparse Signal Representation,”
IEEE Trans. Med. Imag. **28**(12), 1997–2006
(2009)

**13. **B.T. Cox, J.G. Laufer, and P.C. Beard, “The challenges for quantitative
photoacoustic imaging,” Proc.
SPIE **7177**, 717713
(2009)

**14. **A. Marcano, N. Melikechi, and G. Verde, “Shift of the absorption spectrum of
organic dyes due to saturation,” J. Chem.
Phys. **113**(14), 5830–5835
(2000)

**15. **A. Mishra, R.K. Behera, P.K. Behera, B.K. Mishra, and G.B. Behera, “Cyanines during the 1990s: A
review,” Chem. Rev. **100**, 1973–2011
(2000)

**16. **C. Eggeling, J. Widengren, R. Rigler, and C.A.M. Seidel, “Photobleaching of Fluorescent Dyes
under Conditions Used for Single-Molecule Detection: Evidence of
Two-Step Photolysis,” Anal. Chem. **70**(13), 2651–2659
(1998)

**17. **S.-S. Chang, C.-W. Shih, C.-D. Chen, W.-C. Lai, and C.R.C. Wang, “The Shape Transition of Gold
Nanorods,” Langmuir **15**(3), 701–709
(1998)

**18. **S. Link, C. Burda, B. Nikoobakht, and M.A. El-Sayed, “Laser-induced shape changes of
colloidal gold nanorods using femtosecond and nanosecond laser
pulses,” J. Phys. Chem. B **104**, 6152–6163
(2000)

**19. **S.R. Arridge, M. Schweiger, M. Hiraoka, and D.T. Delpy, “A finite element approach for
modelling photon transport in tissue,”
Med. Phys. **20**, 299–309
(1993)

**20. **J.L. Jiménez Pérez, R. Gutierrez Fuentes, J.F. Sanchez Ramirez, and A. Cruz-Orea, “Study of gold nanoparticles effect on
thermal diffusivity of nanofluids based on various solvents by
using thermal lens spectroscopy,” Eur.
Phys. J. Special Topics **153**, 159161
(2008)

**21. **J. Alper and K. Hamad-Schifferli, ”Effect of Ligands on Thermal
Dissipation from Gold Nanorods,”
Langmuir **26**(6), 37863789
(2010)

**22. **B.T. Cox, “Quantitative Photoacoustic Tomography
with Fluence-Dependent Absorbers,” in
*Biomedical Optics*, OSA Technical Digest (CD)
(Optical Society of America,
2010), paper BWG3.