1.

Introduction

Photoacoustic imaging (PAT) is a recent biomedical imaging modality that can provide high-resolution images of optical contrast of heterogeneous media such as biological tissues.^1–13 In a typical PAT experiment, a short pulse of near-infrared (NIR) photons is sent into the medium to be probed. Photons then propagate inside the medium and a portion of them gets absorbed during the propagation process. The energy absorbed by the medium leads to the heating of the medium, and the heating then forces the medium to expand. The medium cools down when the remaining photons exit, and the cooling process leads to the contraction of the medium. The expansion and contraction of the medium initializes a pressure change inside the medium, which then propagates in the form of ultrasound waves. The ultrasound signals on the surface of the medium are then measured. These measurements are used to infer information on the optical properties of the medium.^14–31

We study here the problem of using PAT to image two-photon absorption properties of tissue-like heterogeneous media.^32–42 Here by “two-photon absorption,” we mean the absorption event where an electron transfers to an excited state after simultaneously absorbing two photons whose total energy exceed the electronic energy band gap.^43–45 Even though it occurs less frequently in normal biological tissues than its single-photon counterpart (i.e., the absorption event where an electron transfers to an excited state after absorbing the energy of a single photon), two-photon absorption is extremely useful in practice. In recent years, various types of materials with strong two-photon absorptions have been proposed and engineered as exogenous contrast agents for different optical imaging modalities.^46–49 Many such materials can be tuned to be associated with specific molecular signatures. Therefore, they can be used to visualize particular cellular functions and molecular processes inside biological tissues.

There have been extensive experimental investigations on measuring two-photon absorption properties of various materials using PAT.^{33,34,39–42,50,51} Even though many of these studies use special experimental techniques, they do show collectively the feasibility of experimentally detecting contributions from two-photon absorptions to measured photoacoustic signals. That is, these studies showed that it is indeed possible to have strong enough photoacoustic effect from two-photon absorption that can be experimentally detected. However, it has not been satisfactorily demonstrated so far how to, if it is possible at all, separate the photoacoustic effect due to single-photon absorption from that due to two-photon absorption. In the rest of this paper, we demonstrate, computationally, through a model-based reconstruction algorithm, that it is possible to get quantitative reconstructions of both single-photon and twophoton absorptions and therefore separate them, if indeed the two-photon absorption contribution to the photoacoustic effect can be detected in the measured acoustic data.

2.

Mathematical Models

The main physical processes involved in a two-photon photoacoustic tomography (TP-PAT) experiment are the propagation of NIR photons and the propagation of ultrasound signals in the underlying medium. In optically heterogeneous media such as the biological tissues, it is well established now that the propagation of NIR photons can be modeled with a diffusion equation for the local density of photons.^17,23,30,52 The main difference between TP-PAT and the regular PAT is that two-photon absorption, in addition to single-photon absorption, needs to be considered in the model for light propagation. Let us denote by $\mathrm{\Omega }\subseteq {\mathbb{R}}^d$ ($d\ge 2$) the medium to be probed and denote by $u(\mathbf{x})$ the density of photons at position $\mathbf{x}\in \mathrm{\Omega }$. We then have that $u(\mathbf{x})$ solves the following semilinear diffusion equation:⁵³

The main difference between the nonlinear diffusion equation (1) and the classical linear diffusion equation for modeling light propagation in PAT is the extraterm $\mu (\mathbf{x})|u|u(\mathbf{x})$ that models the two-photon absorption mechanism. It is not the objective of this work to carefully derive the current model from first principles for light propagation in diffusive media. However, we do want to offer the following insights. First, in probabilistic representation, the solution to the diffusion equation, the photon density $u(\mathbf{x})$ (after being normalized by the source strength), is the probability of finding a photon at the location $\mathbf{x}$. Therefore $uu$ is the probability of finding two photons simultaneously at $\mathbf{x}$. This probability, multiplied by the absorption rate $\mu$, is the total two-photon absorption at $\mathbf{x}$. The reason we use $|u|u$, not $uu$, is to ensure that the solution to the resulting nonlinear diffusion equation, that is $u$, is nonnegative, as is required by physics; see more discussions in Ref. 53. Second, one can indeed derive our nonlinear diffusion model from the nonlinear radiative transfer model used in Ref. 32 if the underlying medium is highly scattering media, following classical results in kinetic theory.⁵⁴ Therefore, Eq. (1) is not a completely artificial model.

The nonlinear term $\mu |u|u$ makes the model Eq. (1) harder to solve computationally than the classical linear diffusion model. It is important to emphasize that this nonlinear diffusion model is indeed a well-posed mathematical model that admits a unique solution for a given illumination source $g$ under classical assumptions on regularities of the coefficients and the domain. Classical numerical discretization schemes, such as finite element and finite difference methods, can be used to discretize the equation. Iterative schemes such as the Newton’s method can be used to solve the resulting nonlinear algebraic system; see more detailed discussions in Ref. 53.

The initial pressure field generated by the photoacoustic effect in TP-PAT is the product of the Grüneisen coefficient of the medium, $\mathrm{\Gamma }$, and the total energy absorbed locally by the medium, $\sigma u+\mu |u|u$. Note that here the total absorbed energy consists of two components, the contribution from single-photon absorption, $\sigma u$, and the contribution from two-photon absorption, $\mu |u|u$. Therefore, we write the initial pressure field as⁵³

The change of pressure field generates ultrasound waves that propagate following the standard acoustic wave equation, the same model equation for ultrasound propagation in the regular PAT¹¹

The acoustic datum measured in TP-PAT is the ultrasound signal on the surface of the medium, $p_{|(0,T]\times \partial \mathrm{\Omega }}$, for time $T$ sufficiently long, and very often, we need to measure data that are generated from multiple illumination sources. From the measured data, we are interested in reconstructing the physical coefficients $(\mathrm{\Gamma },\gamma ,\sigma ,\mu )$ of the underlying medium. Note that among all the coefficients, the two-photon absorption coefficient $\mu$ is the only new coefficient that appears in TP-PAT. The coefficients $(\mathrm{\Gamma },\gamma ,\sigma )$ are also quantities to be reconstructed in the regular PAT.^{3,5,13,17,30,52,55–61} Mathematical analysis in Ref. 53 shows that one can not simultaneously reconstruct all four coefficients $(\mathrm{\Gamma },\gamma ,\sigma ,\mu )$ even from data collected from multiple illuminations, if all illumination sources have the same optical wavelength. We will therefore only focus on the two absorption coefficients $(\sigma ,\mu )$ in the rest of the paper. The reconstruction of all four coefficients using multispectral data following the ideas in Refs. 3, 21, 55, 6263.64.65.–66 will be the subject of a future work.

3.

A Two-Step Reconstruction Method

We now present a numerical method for the reconstruction of the absorption coefficients $(\sigma ,\mu )$. We follow the standard two-step procedure in quantitative PAT image reconstructions. In the first step, the qualitative step, we reconstruct the initial pressure field, $H$ in Eq. (2), from measured acoustic data using the wave equation (3). In the second step, the quantitative step, we reconstruct the absorption coefficients from the initial pressure field $H$ using the nonlinear diffusion equation (1). We emphasize that the reason for choosing this two-step reconstruction strategy, instead of the more recent one-step algorithms such as those in Refs. 6768.–69, is that the two-step method allows us to avoid solving the nonlinear diffusion equation (1) in the quantitative reconstruction step in this specific setup; see more discussions in Sec. 3.2.

4.

Numerical Implementations

We now provide some details on the numerical implementation of the two-step algorithm for quantitative TP-PAT image reconstructions. We limit ourselves to two-dimensional simulations. Nothing changes in three-dimensional case besides the increasing of the computational cost. We use the notational convention $\mathbf{x}=(x,y)$ for the spatial variable.

Since the units of the physical quantities in the nonlinear diffusion model (1) and the acoustic wave Eq. (3) are very different, we first normalize the problems by taking the following convention. We take the spatial domain $\mathrm{\Omega }$ to be the unit square $\mathrm{\Omega }={[\mathrm{0,1}]}^2$ and set the ultrasound speed $c=1$. We set the time interval where the measurements are taken as $(0,T]$ with $T=3$. This convention means that if we take the size of $\mathrm{\Omega }$ in units of ${\mathrm{cm}}^2$, the ultrasound speed at $1.5\times {10}^5\mathrm{cm}/\mathrm{s}$, then $T=3$ equals $20\mu \mathrm{s}$. We observed in our numerical simulations (see discussions in the next section) that these choices of $\mathrm{\Omega }$, $c$, and $T$ for the wave equation are sufficient to capture all of the physical wave signals generated by the initial conditions $H(\mathbf{x})$ we have tested, at least up to the numerical discretization errors. For the diffusion problem, we set $U_0={10}^{11}$ as the characteristic photon density involved in the system, and normalize the solution $u$ and the boundary illumination against $U_0$.

The wave equation (3) is posed on ${\mathbb{R}}^2$, not inside $\mathrm{\Omega }$. We therefore have to make a truncation to have a finite domain for the wave simulation. We do this by using the technique of perfectly matched layers (PML).^87,88 We surround our physical domain $\mathrm{\Omega }$ with a PML region of thickness 0.2 to have the computational domain $\tilde{\mathrm{\Omega }}={[-0.2,1.2]}^2$. We use a split-field PML scheme (see, e.g., Refs. 87 and 88). This scheme reduces to the undamped wave equation in the physical domain $\mathrm{\Omega }$, coupled with a damped wave split-field scheme in the PML region $\tilde{\mathrm{\Omega }}\setminus \mathrm{\Omega }$. Ultimately, we end up solving the system of equations, assuming again that $c=1$

We discretize these equations using standard second-order finite differences in space, and first-order finite differences in time on uniform spatial-temporal grids. The spatial grid covering $\tilde{\mathrm{\Omega }}$ consists of $141\times 141$ spatial points

To solve the adjoint wave equation (7), we first perform the change of variable $t^{\prime }=T-t$ to transform the equation into an initial value (instead of final value) problem. We then apply the same type of spatial-temporal discretizations to the new equation.

The nonlinear diffusion equation (1) and the linear diffusion equations involved in the reconstruction process, mainly in Eq. (9), are all discretized using a standard first-order finite element method with about 12,000 elements on a triangular mesh of $\mathrm{\Omega }$. The nonlinear algebra system resulting from the discretization of Eq. (1) is solved with the Newton’s method. In the forward simulation, the initial pressure field $H$ that is needed in the acoustic wave equation (3) is linearly interpolated from the quadrature points of the triangular elements. In the reconstruction process, the initial pressure field $H$ reconstructed in the qualitative step, i.e., the first step, is interpolated back to the quadrature points of the triangular elements as datum for the quantitative reconstruction step. These interpolation processes induce additional noise in the reconstruction process besides the artificial white noise we add to the acoustic data that we discuss below.

To generate synthetic data, we solve the nonlinear diffusion equation (1) with the true physical coefficients to generate $H$ and then solve the acoustic wave equation (3) to produce $p{(t,\mathbf{x})}_{|(0,T]\times \partial \mathrm{\Omega }}$. To add noise to the synthetic data, we use the following strategy. At each point $\mathbf{x}\in \partial \mathrm{\Omega }$ where the ultrasound signal is measured, we generate an independent Gaussian “white-noise” process $w_{\mathbf{x}}(t)$, $t\in (0,T]$, that satisfies

In our numerical simulations in the upcoming section, we set the tolerance level $ϵ={10}^{-6}$ in Algorithm 1 and run this algorithm for a maximum number of $K=1000$ iterations. The quantity ${\Vert \mathcal{A}{H}_k-p^{*}\Vert }_{L^2((0,T]\times \partial \mathrm{\Omega })}^2$ is usually guaranteed to decrease monotonically [see, e.g., the discussion in (Ref. 83, Sec. 7.4)]. However, it is clear from our description above that our implementation of the operators $\mathcal{A}$ and ${\mathcal{A}}^{\top }$ constitute only approximate adjoints of one another due to errors associated with the finite difference approximations and PML region in the wave solver. Therefore, we also force Algorithm 1 to exit if nonmonotonic behavior of ${\Vert \mathcal{A}{H}_k-p^{*}\Vert }_{L^2((0,T]\times \partial \mathrm{\Omega })}^2$ is encountered.

5.

Numerical Simulations

We now present some numerical simulations to demonstrate the performance of our quantitative reconstruction strategy. Our main objective is to show that the photoacoustic data measured in TP-PAT allow quantitative separation of the single-photon and the two-photon absorption coefficients. In all the simulation results below, we set the Grüneisen coefficient $\mathrm{\Gamma }=1$ and the diffusion coefficient $\gamma (\mathbf{x})=0.02+0.01\sin (2\pi y)$. Moreover, we use data collected from four different illumination sources. The first source function takes a constant value the top and right sides of the boundary and is zero everywhere else. That is,

To measure the quality of the reconstructions, we use relative $L^2$ errors. Let $f$ be a quantity to be reconstructed, $f_t$ its true value, and $f_r$ the reconstructed value. Then, the relative $L^2$ error of the reconstruction, denoted by ${\mathcal{E}}_{L^2}(f)$, is the ratio between the size of the error in the reconstruction and the size of the true quantity. That is,

5.1.

Experiment I

In the first group of numerical simulations, we attempt to reconstruct the absorption coefficients $(\sigma ,\mu )$ shown in Fig. 1. We first perform reconstructions, using Algorithm 1, on the initial pressure field $H$ generated by the four illumination sources ${\{g_i\}}_{i=1}^4$ from acoustic data of different noise levels. The quality of the reconstructions, in terms of the relative $L^2$ errors, is summarized in Table 1, third column. We observed, as has been confirmed by many works in the PAT community, the qualitative reconstruction is of high quality. We show in Figs. 2 and 3 some reconstructions with illuminations $g_1$ and $g_2$, respectively. Shown are the true initial pressure field $H$, the reconstructed $H$ using clean data (NSR $\eta =0.0$) and noisy data with NSR $\eta =0.1$. The relative $L^2$ errors of the reconstruction in Fig. 2 are 0.04% for the case of $\eta =0.0$ and 2.7% for the case of $\eta =0.1$, whereas those for the reconstructions in Fig. 3 are 0.06% for the case of $\eta =0.0$ and 2.6% for the case of $\eta =0.1$.

Fig. 1

The true absorption coefficients, (a) $\sigma$ ${(\mathrm{cm}}^{-1})$ and (b) $\mu \times U_0$ $({\mathrm{cm}}^{-1})$, used to generate synthetic data in experiment I.

Table 1

Quality of reconstructions in experiment I. Shown are relative L2 errors in the reconstructions of various initial pressure fields (third column) and the corresponding absorption coefficients in Fig. 1 (fourth column) from ultrasound data with different noise levels (controlled with the NSR η).

NSR	Illumination	EL2(H)	[EL2(σ),EL2(μ)]
$\eta =0.00$	$g_1$	$4.05\times {10}^{-4}$	$(0.46,3.33)\times {10}^{-2}$
	$g_2$	$6.66\times {10}^{-4}$
	$g_3$	$6.22\times {10}^{-4}$
	$g_4$	$4.96\times {10}^{-4}$
$\eta =0.01$	$g_1$	$7.13\times {10}^{-3}$	$(1.71,4.08)\times {10}^{-2}$
	$g_2$	$7.30\times {10}^{-3}$
	$g_3$	$8.25\times {10}^{-3}$
	$g_4$	$8.00\times {10}^{-3}$
$\eta =0.05$	$g_1$	$1.78\times {10}^{-2}$	$(3.80,6.33)\times {10}^{-2}$
	$g_2$	$1.80\times {10}^{-2}$
	$g_3$	$1.77\times {10}^{-2}$
	$g_4$	$1.93\times {10}^{-2}$
$\eta =0.10$	$g_1$	$2.68\times {10}^{-2}$	$(5.15,8.01)\times {10}^{-2}$
	$g_2$	$2.63\times {10}^{-2}$
	$g_3$	$2.48\times {10}^{-2}$
	$g_4$	$2.55\times {10}^{-2}$

Fig. 2

(a) The initial pressure field $H(\mathbf{x})$ generated from illumination $g_1$ using the true absorption coefficients in Fig. 1, (b) as well as the reconstructions of $H$ using ultrasound data with NSR $\eta =0.0$ (clean data) and (c) NSR $\eta =0.1$.

Fig. 3

The same as Fig. 2 except that the illumination source used is $g_2$.

Let us mention that even though we can see clearly the two-photon absorbing inclusions in the true initial pressure field $H$ and the reconstructed $H$, the true absorption coefficients in Fig. 1 are very different from the $H$ in Figs. 2 and 3. In other words, knowing $H$ does not provide us enough information about the true absorption coefficients unless we perform the next step, the quantitative step, of the reconstruction.

In Fig. 4, we show the reconstructions of the coefficient pair $(\sigma ,\mu )$ in Fig. 1 from the initial pressure fields we obtained in the qualitative step, using Algorithm 2. Shown, from left to right, are reconstructions using noisy data with NSR $\eta =0.00$, $\eta =0.05$, and $\eta =0.10$, respectively. The quality of the reconstructions is very high with relative $L^2$ errors $(1.71,4.08)\times {10}^{-2}$, $(3.80,6.33)\times {10}^{-2}$, and $(5.15,8.01)\times {10}^{-2}$, respectively; see the fourth column of Table 1 for the reconstruction result using data with $\eta =0.01$, which we did not show here since it is too similar to the case with $\eta =0.00$.

Fig. 4

The absorption coefficients (a) $\sigma$ and (b) $\mu \times U_0$ reconstructed using noisy data with NSR $\eta =0.00$, $\eta =0.05$, and $\eta =0.10$ (from left to right).

The reconstructions in Fig. 4 show that by performing quantitative reconstructions, we can separate the two-photon absorption coefficient from the single-photon absorption coefficient from the initial pressure field. This is clearly important for practical applications of TP-PAT where two-photon absorption is the main quantity of interest.

5.2.

Experiment II

In the second group of numerical simulations, we study the reconstruction of the absorption coefficients $(\sigma ,\mu )$ shown in Fig. 5. In Fig. 6, we present the true initial pressure filed $H$ computed with illumination source $g_1$, and the reconstructions of this $H$ with clean ultrasound data (NSR $\eta =0.00$) and noisy data (NSR $\eta =0.10$). By visual inspection, we can see the presence of both the single-photon absorption and the two-photon absorption inclusions in $H$. The reconstructions are impressively good, with relative $L^2$ errors $9.30\times {10}^{-4}$ and $2.51\times {10}^{-2}$ for $\eta =0.00$ and $\eta =0.10$, respectively, even when $H$ is this complicated. We have also performed similar reconstructions for $H$ generated from the other three illumination sources $g_2$, $g_3$, and $g_4$. The relative errors in the reconstructions are summarized in Table 2, third column. Despite its slight degeneration as noise level increases, the quality of the reconstructions of $H$ remains high at moderate noise levels.

Fig. 5

The true absorption coefficients, (a) $\sigma$ ${(\mathrm{cm}}^{-1})$ and (b) $\mu \times U_0$ $({\mathrm{cm}}^{-1})$, used to generate synthetic data in experiment II.

Fig. 6

(a) The initial pressure field $H(\mathbf{x})$ generated from illumination $g_1$ using the true absorption coefficients in Fig. 5, (b) as well as the reconstructions of $H$ using ultrasound data with NSR $\eta =0.00$ (clean data) and (c) NSR $\eta =0.10$.

Table 2

Quality of reconstructions in experiment II. Shown are relative L2 errors in the reconstructions of various initial pressure fields (third column) and the corresponding absorption coefficients in Fig. 5 (fourth column) from ultrasound data with different noise levels.

NSR	Illumination	EL2(H)	[EL2(σ),EL2(μ)]
$\eta =0.00$	$g_1$	$9.30\times {10}^{-4}$	$(2.19,5.26)\times {10}^{-2}$
	$g_2$	$1.16\times {10}^{-3}$
	$g_3$	$1.17\times {10}^{-3}$
	$g_4$	$9.38\times {10}^{-4}$
$\eta =0.01$	$g_1$	$8.09\times {10}^{-3}$	$(2.77,5.78)\times {10}^{-2}$
	$g_2$	$8.25\times {10}^{-3}$
	$g_3$	$7.79\times {10}^{-3}$
	$g_4$	$7.68\times {10}^{-3}$
$\eta =0.05$	$g_1$	$1.72\times {10}^{-2}$	$(4.31,7.73)\times {10}^{-2}$
	$g_2$	$1.82\times {10}^{-2}$
	$g_3$	$1.82\times {10}^{-2}$
	$g_4$	$1.78\times {10}^{-2}$
$\eta =0.10$	$g_1$	$2.51\times {10}^{-2}$	$(5.60,9.34)\times {10}^{-2}$
	$g_2$	$2.48\times {10}^{-2}$
	$g_3$	$2.40\times {10}^{-2}$
	$g_4$	$2.79\times {10}^{-2}$

To separate $\mu$ from $\sigma$ in the initial pressure fields, we perform quantitative reconstructions using Algorithm 2. In Figs. 7(a)–7(c), we show the reconstructions from data with NSR $\eta =0.00$, $\eta =0.05$, and $\eta =0.10$, respectively. The relative $L^2$ errors in the reconstructions are $(2.19,5.26)\times {10}^{-2}$, $(4.31,7.73)\times {10}^{-2}$, and $(5.60,9.34)\times {10}^{-2}$, respectively. Once again, we see good separation of the two different absorption coefficients, which were mixed together in the initial pressure fields $H$ in Fig. 6. The quantitative reconstruction results are summarized in Table 2, fourth column.

Fig. 7

Reconstruction of the absorption coefficient pair (a) $\sigma$ and (b) $\mu \times U_0$ in Fig. 5 using data at different noise levels (NSR $\eta =0.00$, $\eta =0.05$, and $\eta =0.10$ from left to right).

6.

Concluding Remarks

We studied in this paper quantitative image reconstructions in TP-PAT, aiming at reconstructing the single-photon absorption and the two-photon absorption coefficients of biological tissues from measured ultrasound signals generated by the photoacoustic effect of light absorption. We introduced a nonlinear diffusion equation as the model for light propagation in TP-PAT and presented a two-step image reconstruction strategy, including a noniterative quantitative reconstruction step, based on this model. We showed, with computational simulations, that while single-photon absorption and two-photon absorption are mixed in the images of the initial pressure fields, they can be stably separated from each other through the quantitative reconstruction step, using Algorithm 2, even when the ultrasound data contain relatively high level of random noises. Our numerical simulations confirm the results of mathematical analysis of the problem in a previous publication.⁵³

Compared to the case in the regular PAT, quantitative image reconstruction in TP-PAT is far less investigated, theoretically or computationally, to date. Our numerical simulations show great promise in the quantitative imaging of the two-photon absorption. However, there are still a lot of issues that need to be addressed. For instance, it would be very interesting to test the two-step reconstruction method we proposed against experimentally measured data to see what types of resolution and contrast we can get for the two-photon absorption coefficient. It would also be important to develop efficient algorithms to reconstruct the diffusion coefficient $\gamma$ in addition to the absorption coefficients. Last but not the least, reconstructing the Grüneisen coefficient with multispectral data, following for instance the ideas in Refs. 3, 21, 55, 6263.64.65.–66, could also be extremely useful as well.