Translator Disclaimer
2 January 2018 Quantitative photoacoustic imaging of two-photon absorption
Author Affiliations +
Photoacoustic tomography (PAT) is a hybrid imaging modality where we intend to reconstruct optical properties of heterogeneous media from measured ultrasound signals generated by the photoacoustic effect. In recent years, there have been considerable interests in using PAT to image two-photon absorption, in addition to the usual single-photon absorption, inside diffusive media. We present a mathematical model for quantitative image reconstruction in two-photon photoacoustic tomography (TP-PAT). We propose a computational strategy for the reconstruction of the optical absorption coefficients and provide some numerical evidences based on synthetic photoacoustic acoustic data to demonstrate the feasibility of quantitative reconstructions in TP-PAT.



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.113 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.1431

We study here the problem of using PAT to image two-photon absorption properties of tissue-like heterogeneous media.3242 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.4345 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.4649 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,3942,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.


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 ΩRd (d2) the medium to be probed and denote by u(x) the density of photons at position xΩ. We then have that u(x) solves the following semilinear diffusion equation:53

Eq. (1)

·γ(x)u(x)+σ(x)u(x)+μ(x)|u|u(x)=0,in  Ωu+κγuν=g(x),on  Ω,
where is the usual gradient operator with respect to the spatial variable x, and the function g(x) models the incoming NIR illumination source on the boundary Ω. The function γ(x) is the diffusion coefficient of the medium, and the functions σ(x) and μ(x) are the single-photon absorption and the (intrinsic) two-photon absorption coefficients, respectively. The unit outer normal vector at point x on the boundary Ω is denoted by ν(x), and the notation uν=ν·u is used in the Robin boundary condition. The coupling parameter κ in the boundary condition is the rescaled extrapolation length. Its value depends on many parameters and can be explicitly calculated in specific settings.54

The main difference between the nonlinear diffusion equation (1) and the classical linear diffusion equation for modeling light propagation in PAT is the extraterm μ(x)|u|u(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(x) (after being normalized by the source strength), is the probability of finding a photon at the location x. Therefore uu is the probability of finding two photons simultaneously at x. This probability, multiplied by the absorption rate μ, is the total two-photon absorption at 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.54 Therefore, Eq. (1) is not a completely artificial model.

The nonlinear term μ|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, Γ, and the total energy absorbed locally by the medium, σu+μ|u|u. Note that here the total absorbed energy consists of two components, the contribution from single-photon absorption, σu, and the contribution from two-photon absorption, μ|u|u. Therefore, we write the initial pressure field as53

Eq. (2)

where the Grüneisen coefficient is nondimensionalized, and it describes the efficiency of the photoacoustic effect of the underlying medium.

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 PAT11

Eq. (3)

1c2(x)2pt2Δp=0,in  (0,+)×Rdp(t,x)=χΩH(x),pt(t,x)=0,in  {t=0}×Rd,
where p(t,x) is the pressure field, and c(x) is the speed of the ultrasound waves. In most biological applications of PAT, the ultrasound speed c is assumed known and is often taken as the speed of ultrasound in water since most biological tissues behave like water to ultrasound waves. The function χΩ is the characteristic function of the domain Ω. It should be understood as the extension operator that extends the initial pressure field inside the medium Ω to the whole space Rd, that is,

The acoustic datum measured in TP-PAT is the ultrasound signal on the surface of the medium, p|(0,T]×Ω, 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 (Γ,γ,σ,μ) of the underlying medium. Note that among all the coefficients, the two-photon absorption coefficient μ is the only new coefficient that appears in TP-PAT. The coefficients (Γ,γ,σ) are also quantities to be reconstructed in the regular PAT.3,5,13,17,30,52,5561 Mathematical analysis in Ref. 53 shows that one can not simultaneously reconstruct all four coefficients (Γ,γ,σ,μ) 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 (σ,μ) 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.


A Two-Step Reconstruction Method

We now present a numerical method for the reconstruction of the absorption coefficients (σ,μ). 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.


Qualitative Step: Reconstructing Initial Pressure Fields

In the qualitative step of TP-PAT, we aim at reconstructing the initial pressure field H from measured datum p|(0,T]×Ω. This step is the same as that in the regular PAT, which has been extensively studied in the past. Many efficient algorithms have been proposed; see for instance Refs. 4, 27, 31, 7071. for an incomplete list of works in this direction. We implement here a simple least-square-based algorithm for the reconstruction.

To simplify the presentation, let us denote by A the linear operator that takes the initial pressure field H(x) to the acoustic field on the boundary Ω, i.e.,

Eq. (4)

Our objective is to invert the operator A to find H for a given measurement p*(t,x)|(0,T]×Ω. We solve this problem in the least-square sense, that is, we search for H as the minimizer of the misfit functional

Eq. (5)

Standard least-square theory then implies that the minimizer H solves the normal equation

Eq. (6)

where A denotes the L2-adjoint of the operator A. We therefore need to invert the self-adjoint operator AA to find H.

We solve the normal equation (6) using the conjugate gradient method.83 In a nutshell, the method seeks a solution of the normal equation by iteratively choosing “conjugate” or “AA-orthogonal” directions, and minimizing the magnitude of the residual, AHp*L2((0,T]×Ω)2, in each of these conjugate directions.

More precisely, let Hk be the value of H at iteration k and let {hi}i=1k be the set of “AA-orthogonal” directions constructed in the first k iterations. The directions {hi}i=1k satisfy the AA orthogonality relation,   2ik

hi,AAhjL2(Ω)=0,  1ji1.

We now search for an update of Hk in the direction hk such that the residual is minimized after the update. That is, we minimize the residual Ψ(α) over α with

where it was recalled p*=AH. The optimality condition immediately gives that the step length at iteration k is
αk=hk,AA(HkH)L2(Ω)hk,AAhkL2(Ω)=hk,skL2(Ω)hk,AAhkL2(Ω),with  sk=A[A(HkH)].
Note that if we define rk=A(HkH)=AHkp* as the residual of the original problem at step k, then sk=Ark is simply the so-called normal residual corresponding to rk. The updated value Hk+1 is then obtained as
while the normal equation residual sk is updated as

The conjugate gradient method updates the search direction following:

We summarize the conjugate gradient method in Algorithm 1 following the routine in Ref. 83, with an accuracy tolerance parameter ϵ>0 and the maximal number of iteration K.

Algorithm 1

CG algorithm for qualitative reconstruction.

1: Set parameters ϵ and K; set k=0
2: Set initial guess H=0
3: Evaluate the residual r=p*AH and the normal residual s=Ar
4: Set initial search directions h=s
5: Evaluate the size of normal residual γ=sL2(Ω)2
6: whilekK and γ/Ap*L2(Ω)2>ϵdo
7: g=Ah
8: α=γ/gL2((0,T]×Ω)2
9: H=H+αh
10: r=rαg, s=Ar
11: β=sL2(Ω)2/γ
12: γ=sL2(Ω)2
13: h=s+βh
14: k=k+1
15: end while

It is often the case that a regularization term is added to the misfit functional ϕ(H) in Eq. (5). In our implementation, we did not include a regularization term in ϕ(H). When it is needed, the two algorithmic parameters ϵ and K can both serve as mechanisms to regularize the reconstruction. We did not pursue in this direction in this study, but we understand that tuning regularization can refine some of the reconstruction results that we show in the next section.

Let us mention that even though the operator A and its adjoint operator A are called in each iteration of the algorithm, these operators are never explicitly formed in the numerical implementation. We only need to know the actions of these operators on given vectors. For instance, to evaluate Ah for a given h(x), we solve the acoustic wave equation (3) with initial condition p(0,x)=h(x) and record the solution on the boundary of Ω: p(t,x)|(0,T]×Ω. To evaluate Ar for a given r(t,x), we first solve the following adjoint wave equation:

Eq. (7)

1c2(x)2vt2Δv=0,in  (0,T)×Ωv(t,x)=0,vν(t,x)=r,in  (0,T)×Ωv(t,x)=0,vt(t,x)=0,in  {t=T}×Ω.
We then take Ar=vt(0,x). The derivation of Eq. (7) is straightforward and has been documented previously,67,84 so we omit the details here.


Quantitative Step: Reconstructing Absorption Coefficients

The second step, the quantitative step, is to reconstruct the optical coefficients from the initial pressure field H recovered in the first step. In recent years, this step was the subject of many computational studies in the case of the regular PAT; see, for instance, Refs. 3, 14, 16, 17, 23, 30, 52, 5556., 8485.86 for a partial list of references.

Our main objective here is to develop an algorithm to reconstruct the absorption coefficients (σ,μ) in TP-PAT to show that we can separate two-photon absorption from single-photon absorption. We assume that both the Grüneisen coefficient Γ and the diffusion coefficient γ are known already, for instance from a regular PAT reconstruction.

We assume that we have reconstructed initial pressure fields generated from J2 illuminations sources. We denote by {Hj=Γ[σuj+μ|uj|uj]}j=1J those initial pressure fields, where uj denotes the solution of the nonlinear diffusion equation with illumination sources gj (1jJ).

We first reconstruct from the initial pressure fields {Hj}j=1J, using the fact that Γ and γ are known, the quantities

Eq. (8)

This allows us to replace the term σuj+μ|uj|uj in the nonlinear diffusion equation (1) for source j to obtain the following linear diffusion equation for uj (1jJ)

Eq. (9)

·(γuj)=HjΓ  in  Ω,uj+κujν=gj  on  Ω.
We can solve this linear elliptic equation to reconstruct uj, again since Γ and γ are known. Therefore, we can reconstruct the quantities

Eq. (10)

Therefore, at each point xΩ, we have the following system to determine σ and μ
We then reconstruct (σ,μ) by solving this small linear system, in least-square sense, at each point xΩ, to get

Eq. (11)

We have assumed here that the small 2×2 matrix
in Eq. (11) is invertible at each point xΩ. Theoretical analysis in Ref. 53 shows that one can indeed invert this matrix if the illuminations are selected carefully, that is, the illuminations are sufficiently different from each other. In our numerical experiments, we observe that this matrix is invertible for almost all illuminations that we have tried.

We now summarize the quantitative reconstruction step in Algorithm 2.

Algorithm 2

Noniterative algorithm for quantitative reconstruction.

1: forj1,Jdo
2: Reconstruct the quantities σuj+μ|uj|uj following Eq. (8)
3: end for
4: forj1,Jdo
5: Solve the diffusion equation (9) to reconstruct uj
6: end for
7: forj1,Jdo
8: Reconstruct the quantities σ+μ|uj| following Eq. (10)
9: end for
10: for each point xΩdo
11: Evaluate ω=j=1J|uj(x)|, θ=j=1J|uj(x)|2, ξ=j=1JHjΓuj and ζ=j=1JHj|uj|Γuj
12: Evaluate (σ,μ) using the formula: [σ(x)μ(x)]=(Jωωθ)1(ξζ)
13: end for

Let us emphasize two important features of the quantitative reconstruction algorithm. First, even though the reconstruction of the coefficients (σ,μ) from initial pressure field H is a nonlinear inverse problem, our reconstruction method is noniterative. Therefore, there is no convergence issues at all. The method is guaranteed to give the correct reconstruction result. Second, the reconstruction algorithm is computationally cheap. The major computational cost of the reconstruction algorithm is the solution of the J linear diffusion equations in Eq. (9). The cost in dealing with the algebraic calculations in the rest of the algorithm, in Eqs. (8), (10), and (11), is almost negligible.


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 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 Ω to be the unit square Ω=[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 Ω in units of cm2, the ultrasound speed at 1.5×105  cm/s, then T=3 equals 20  μs. We observed in our numerical simulations (see discussions in the next section) that these choices of Ω, c, and T for the wave equation are sufficient to capture all of the physical wave signals generated by the initial conditions H(x) we have tested, at least up to the numerical discretization errors. For the diffusion problem, we set U0=1011 as the characteristic photon density involved in the system, and normalize the solution u and the boundary illumination against U0.

The wave equation (3) is posed on R2, not inside Ω. 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 Ω with a PML region of thickness 0.2 to have the computational domain Ω˜=[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 Ω, coupled with a damped wave split-field scheme in the PML region Ω˜Ω. Ultimately, we end up solving the system of equations, assuming again that c=1

where p=px+py, τx(x) and τy(x) are absorptive terms supported only in the PML region Ω˜Ω. In our simulations, we use τx(x)=τx(x,y)=χx>1(x)α(x1)2+χx<0(x)αx2 with α a given constant. Similarly, τy(x,y)=τx(y,x). Initial and boundary conditions can be transformed into this first-order formulation in a straightforward way.

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 Ω˜ consists of 141×141 spatial points

and the temporal grid covering [0,T]=[0,3] consists of 3001 grid points
The velocity fields vx and vy are solved for at staggered half-time steps, i.e., at times {tk+1/2}, using values of the split pressure fields px and py at the usual time steps {tk}. The pressure field p=px+py is then updated using the velocity fields at these staggered times. Hence, the scheme is known as a leapfrog scheme and reduces to standard second-order finite difference time stepping in the physical domain Ω where τx=τy=0. With our spatial step size h=1/100 and our temporal step size Δt=1/1000, we clearly satisfy the CFL stability condition
which is necessary for the explicit finite difference scheme we implemented to be stable.

To solve the adjoint wave equation (7), we first perform the change of variable t=Tt 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 Ω. 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,x)|(0,T]×Ω. To add noise to the synthetic data, we use the following strategy. At each point xΩ where the ultrasound signal is measured, we generate an independent Gaussian “white-noise” process wx(t), t(0,T], that satisfies

where E denotes the operation of taking expectation. We then scale the white noise wx(t) according to the power of the signal p(t,x), to generate noisy data p˜(t,x) with a specified noise-to-signal ratio (NSR) η

In our numerical simulations in the upcoming section, we set the tolerance level ϵ=106 in Algorithm 1 and run this algorithm for a maximum number of K=1000 iterations. The quantity AHkp*L2((0,T]×Ω)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 A and A 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 AHkp*L2((0,T]×Ω)2 is encountered.


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 Γ=1 and the diffusion coefficient γ(x)=0.02+0.01sin(2π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,

The second to fourth sources, g2, g3, and g4, are obtained by rotating g1 by π/2, π, and 3π/2, respectively, along the boundary. Note again that g1 is normalized against U0 already.

To measure the quality of the reconstructions, we use relative L2 errors. Let f be a quantity to be reconstructed, ft its true value, and fr the reconstructed value. Then, the relative L2 error of the reconstruction, denoted by EL2(f), is the ratio between the size of the error in the reconstruction and the size of the true quantity. That is,



Experiment I

In the first group of numerical simulations, we attempt to reconstruct the absorption coefficients (σ,μ) shown in Fig. 1. We first perform reconstructions, using Algorithm 1, on the initial pressure field H generated by the four illumination sources {gi}i=14 from acoustic data of different noise levels. The quality of the reconstructions, in terms of the relative L2 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 g1 and g2, respectively. Shown are the true initial pressure field H, the reconstructed H using clean data (NSR η=0.0) and noisy data with NSR η=0.1. The relative L2 errors of the reconstruction in Fig. 2 are 0.04% for the case of η=0.0 and 2.7% for the case of η=0.1, whereas those for the reconstructions in Fig. 3 are 0.06% for the case of η=0.0 and 2.6% for the case of η=0.1.

Fig. 1

The true absorption coefficients, (a) σ (cm1) and (b) μ×U0 (cm1), 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 η).


Fig. 2

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


Fig. 3

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


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 (σ,μ) 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 η=0.00, η=0.05, and η=0.10, respectively. The quality of the reconstructions is very high with relative L2 errors (1.71,4.08)×102, (3.80,6.33)×102, and (5.15,8.01)×102, respectively; see the fourth column of Table 1 for the reconstruction result using data with η=0.01, which we did not show here since it is too similar to the case with η=0.00.

Fig. 4

The absorption coefficients (a) σ and (b) μ×U0 reconstructed using noisy data with NSR η=0.00, η=0.05, and η=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.


Experiment II

In the second group of numerical simulations, we study the reconstruction of the absorption coefficients (σ,μ) shown in Fig. 5. In Fig. 6, we present the true initial pressure filed H computed with illumination source g1, and the reconstructions of this H with clean ultrasound data (NSR η=0.00) and noisy data (NSR η=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 L2 errors 9.30×104 and 2.51×102 for η=0.00 and η=0.10, respectively, even when H is this complicated. We have also performed similar reconstructions for H generated from the other three illumination sources g2, g3, and g4. 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) σ (cm1) and (b) μ×U0 (cm1), used to generate synthetic data in experiment II.


Fig. 6

(a) The initial pressure field H(x) generated from illumination g1 using the true absorption coefficients in Fig. 5, (b) as well as the reconstructions of H using ultrasound data with NSR η=0.00 (clean data) and (c) NSR η=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.


To separate μ from σ 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 η=0.00, η=0.05, and η=0.10, respectively. The relative L2 errors in the reconstructions are (2.19,5.26)×102, (4.31,7.73)×102, and (5.60,9.34)×102, 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) σ and (b) μ×U0 in Fig. 5 using data at different noise levels (NSR η=0.00, η=0.05, and η=0.10 from left to right).



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.53

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 γ 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.


The authors have no relevant financial interests in this article and no potential conflicts of interest to disclose.


We would like to thank Prof. John C. Schotland (University of Michigan) for useful discussions on the two-photon absorption photoacoustic imaging. This work was partially supported by the National Science Foundation through grants DMS-1321018 and DMS-1620473.



P. Beard, “Biomedical photoacoustic imaging,” Interface Focus, 1 602 –631 (2011). Google Scholar


P. Burgholzer, H. Grun and A. Sonnleitner, “Photoacoustic tomography: sounding out fluorescent proteins,” Nat. Photonics, 3 378 –379 (2009). NPAHBY 1749-4885 Google Scholar


B. T. Cox, J. G. Laufer and P. C. Beard, “The challenges for quantitative photoacoustic imaging,” Proc. SPIE, 7177 717713 (2009). PSISDG 0277-786X Google Scholar


P. Kuchment and L. Kunyansky, “Mathematics of thermoacoustic tomography,” Eur. J. Appl. Math., 19 191 –224 (2008). 0956-7925 Google Scholar


C. Li and L. Wang, “Photoacoustic tomography and sensing in biomedicine,” Phys. Med. Biol., 54 R59 –R97 (2009). PHMBA7 0031-9155 Google Scholar


A. A. Oraevsky, S. L. Jacques and F. K. Tittel, “Measurement of tissue optical properties by time-resolved detection of laser-induced transient stress,” Appl. Opt., 36 402 –415 (1997). APOPAI 0003-6935 Google Scholar


A. A. Oraevsky, A. A. Karabutov, “Optoacoustic tomography,” Handbook of Optical Biomedical Diagonstics, SPIE Press, Bellingham, Washington (2002). Google Scholar


S. K. Patch and O. Scherzer, “Photo- and thermo-acoustic imaging,” Inverse Prob., 23 S1 –S10 (2007). INPEEY 0266-5611 Google Scholar


O. Scherzer, Handbook of Mathematical Methods in Imaging, Springer-Verlag, Berlin (2010). Google Scholar


E. M. Strohm, M. J. Moore and M. C. Kolios, “Single cell photoacoustic microscopy: a review,” IEEE J. Sel. Top. Quantum Electron., 23 137 –151 (2016). Google Scholar


Photoacoustic Imaging and Spectroscopy, Taylor & Francis, Oxford (2009). Google Scholar


L. V. Wang, “Photoacoustic tomography,” Scholarpedia, 9 10278 (2014). 1941-6016 Google Scholar


J. Xia, J. Yao and L.V. Wang, “Photoacoustic tomography: principles and advances,” Electromagn. Waves, 147 1 –22 (2014). PELREX 1043-626X Google Scholar


S. Arridge et al., “Accelerated high-resolution photoacoustic tomography via compressed sensing,” Phys. Med. Biol., 61 8908 –8940 (2016). PHMBA7 0031-9155 Google Scholar


J. R. Cook, W. Frey and S. Emelianov, “Quantitative photoacoustic imaging of nanoparticles in cells and tissues,” ACS Nano, 7 1272 –1280 (2013). ANCAC3 1936-0851 Google Scholar


B. Cox et al., “Quantitative spectroscopic photoacoustic imaging: a review,” J. Biomed. Opt., 17 061202 (2012). JBOPFO 1083-3668 Google Scholar


B. T. Cox et al., “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). APOPAI 0003-6935 Google Scholar


X. Gao et al., “Quantitative imaging of microvasculature in deep tissue with a spectrum-based photo-acoustic microscopy,” Opt. Lett., 40 970 –973 (2015). OPLEDP 0146-9592 Google Scholar


W. Poon et al., “Determination of biodistribution of ultrasmall, near-infrared emitting gold nanoparticles by photoacoustic and fluorescence imaging,” J. Biomed. Opt., 20 066007 (2015). JBOPFO 1083-3668 Google Scholar


D. Razansky and V. Ntziachristos, “Hybrid photoacoustic fluorescence molecular tomography using finite-element-based inversion,” Med. Phys., 34 4293 –4301 (2007). MPHYA6 0094-2405 Google Scholar


D. Razansky, C. Vinegoni and V. Ntziachristos, “Multispectral photoacoustic imaging of fluorochromes in small animals,” Opt. Lett., 32 2891 –2893 (2007). OPLEDP 0146-9592 Google Scholar


J. Ripoll and V. Ntziachristos, “Quantitative point source photoacoustic inversion formulas for scattering and absorbing media,” Phys. Rev. E, 71 031912 (2005). Google Scholar


P. Shao, T. Harrison and R. J. Zemp, “Iterative algorithm for multiple illumination photoacoustic tomography (MIPAT) using ultrasound channel data,” Biomed. Opt. Express, 3 3240 –3249 (2012). BOEICL 2156-7085 Google Scholar


T. Tarvainen et al., “Reconstructing absorption and scattering distributions in quantitative photoacoustic tomography,” Inverse Prob., 28 084009 (2012). Google Scholar


C. Tian et al., “Imaging and sensing based on dual-pulse nonlinear photoacoustic contrast: a preliminary study on fatty liver,” Opt. Lett., 40 2253 –2256 (2015). OPLEDP 0146-9592 Google Scholar


P. K. Upputuri, M. Krisnan and M. Pramanik, “Microsphere enabled subdiffraction-limited optical-resolution photoacoustic microscopy: a simulation study,” J. Biomed. Opt., 23 045001 (2017). JBOPFO 1083-3668 Google Scholar


K. Wang and M. A. Anastasio, “A simple Fourier transform-based reconstruction formula for photoacoustic computed tomography with a circular or spherical measurement geometry,” Phys. Med. Biol., 57 N493 (2012). PHMBA7 0031-9155 Google Scholar


K. Wang et al., “Investigation of iterative image reconstruction in three-dimensional optoacoustic tomography,” Phys. Med. Biol., 57 5399 –5423 (2012). PHMBA7 0031-9155 Google Scholar


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


R. J. Zemp, “Quantitative photoacoustic tomography with multiple optical sources,” Appl. Opt., 49 3566 –3572 (2010). APOPAI 0003-6935 Google Scholar


B. Zhu and E. M. Sevick-Muraca, “Reconstruction of sectional images in frequency-domain based photoacoustic imaging,” Opt. Express, 19 23286 –23297 (2011). OPEXFF 1094-4087 Google Scholar


P. Chantharasupawong, R. Philip and J. Thomas, “Simultaneous optical and photoacoustic measurement of nonlinear absorption,” Appl. Phys. Lett., 102 041116 (2013). APPLAB 0003-6951 Google Scholar


Y.-H. Lai et al., “Nonlinear photoacoustic microscopy via a loss modulation technique: from detection to imaging,” Opt. Express, 23 525 –536 (2014). OPEXFF 1094-4087 Google Scholar


G. Langer et al., “Two-photon absorption-induced photoacoustic imaging of rhodamine B dyed polyethylene spheres using a femtosecond laser,” Opt. Express, 21 22410 –22422 (2013). OPEXFF 1094-4087 Google Scholar


G. Langer et al., “Two-photon absorption-induced photoacoustic and luminescence imaging employing a femtosecond laser,” Proc. SPIE, 8943 89431V (2014). PSISDG 0277-786X Google Scholar


G. J. Tserevelakis et al., “Hybrid multiphoton and optoacoustic microscope,” Opt. Lett., 39 1819 –1822 (2014). OPLEDP 0146-9592 Google Scholar


B. E. Urban et al., “Investigating femtosecond-laser-induced two-photon photoacoustic generation,” J. Biomed. Opt., 19 085001 (2014). JBOPFO 1083-3668 Google Scholar


P. W. Winter et al., “Two-photon instant structured illumination microscopy improves the depth penetration of super-resolution imaging in thick scattering samples,” Optica, 1 181 –191 (2014). Google Scholar


Y. Yamaoka, M. Nambu and T. Takamatsu, “Frequency-selective multiphoton-excitation-induced photoacoustic microscopy (MEPAM) to visualize the cross sections of dense objects,” Proc. SPIE, 7524 75642O (2010). PSISDG 0277-786X Google Scholar


Y. Yamaoka, M. Nambu and T. Takamatsu, “Fine depth resolution of two-photon absorption-induced photoacoustic microscopy using low-frequency bandpass filtering,” Opt. Express, 19 13365 –13377 (2011). OPEXFF 1094-4087 Google Scholar


Y. Yamaoka and T. Takamatsu, “Enhancement of multiphoton excitation-induced photoacoustic signals by using gold nanoparticles surrounded by fluorescent dyes,” Proc. SPIE, 7177 71772A (2009). PSISDG 0277-786X Google Scholar


C. S. Yelleswarapu and S. R. Kothapalli, “Nonlinear photoacoustics for measuring the nonlinear optical absorption coefficient,” Opt. Express, 18 9020 –9025 (2010). OPEXFF 1094-4087 Google Scholar


H. Mahr, “Two-photon absorption spectroscopy,” Quantum Electronics, Academic Press, Cambridge, Massachusetts (2012). Google Scholar


M. Rumi and J. W. Perry, “Two-photon absorption: an overview of measurements and principles,” Adv. Opt. Photonics, 2 451 –518 (2010). AOPAC7 1943-8206 Google Scholar


A. C. Tam and C. K. Patel, “Two-photon absorption spectra and cross-section measurements in liquids,” Nature, 280 304 –306 (1979). Google Scholar


W. Denk, J. Strickler and W. Webb, “Two photon laser scanning fluorescence microscopy,” Science, 248 73 –76 (1990). SCIEAS 0036-8075 Google Scholar


P. T. C. So, “Two-photon fluorescence light microscopy,” Encyclopedia of Life Sciences, Nature Publishing Group, London (2002). Google Scholar


J. Ying, F. Liu and R. R. Alfano, “Spatial distribution of two-photon-excited fluorescence in scattering media,” Appl. Opt., 38 224 –229 (1999). APOPAI 0003-6935 Google Scholar


W. R. Zipfel, R. M. Williams and W. Webb, “Nonlinear magic: multiphoton microscopy in the biosciences,” Nat. Biotechnol., 21 1369 –1377 (2003). NABIF9 1087-0156 Google Scholar


Y. Bae, J. J. Song and Y. B. Kim, “Photoacoustic study of two-photon absorption in hexagonal ZnS,” J. Appl. Phys., 53 615 –619 (1982). JAPIAU 0021-8979 Google Scholar


Y. Yamaoka et al., “Photoacoustic microscopy using ultrashort pulses with two different pulse durations,” Opt. Express, 23 17063 –17072 (2014). OPEXFF 1094-4087 Google Scholar


G. Bal and K. Ren, “Multi-source quantitative PAT in diffusive regime,” Inverse Prob., 27 075003 (2011). INPEEY 0266-5611 Google Scholar


K. Ren and R. Zhang, “Nonlinear quantitative photoacoustic tomography with two-photon absorption,” SIAM J. Appl. Math., 78 (2018). SMJMAP 0036-1399 Google Scholar


R. Dautray and J.-L. Lions, Mathematical Analysis and Numerical Methods for Science and Technology, 6 Springer-Verlag, Berlin (1993). Google Scholar


J. Laufer et al., “Quantitative determination of chromophore concentrations from 2d photoacoustic images using a nonlinear model-based inversion scheme,” Appl. Opt., 49 1219 –1233 (2010). APOPAI 0003-6935 Google Scholar


A. Pulkkinen et al., “A Bayesian approach to spectral quantitative photoacoustic tomography,” Inverse Prob., 30 065012 (2014). INPEEY 0266-5611 Google Scholar


K. Ren, H. Gao and H. Zhao, “A hybrid reconstruction method for quantitative photoacoustic imaging,” SIAM J. Imaging Sci., 6 32 –55 (2013). Google Scholar


K. Ren, R. Zhang and Y. Zhong, “Inverse transport problems in quantitative PAT for molecular imaging,” Inverse Prob., 31 125012 (2015). INPEEY 0266-5611 Google Scholar


T. Saratoon et al., “A gradient-based method for quantitative photoacoustic tomography using the radiative transfer equation,” Inverse Prob., 29 075006 (2013). INPEEY 0266-5611 Google Scholar


J. Zhang and M. A. Anastasio, “Reconstruction of speed-of-sound and electromagnetic absorption distributions in photoacoustic tomography,” Proc. SPIE, 6086 608619 (2006). PSISDG 0277-786X Google Scholar


X. Zhang et al., “Forward-backward splitting method for quantitative photoacoustic tomography,” Inverse Prob., 30 125012 (2014). INPEEY 0266-5611 Google Scholar


G. Bal and K. Ren, “On multi-spectral quantitative photoacoustic tomography in diffusive regime,” Inverse Prob., 28 025010 (2012). INPEEY 0266-5611 Google Scholar


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). JOAOD6 0740-3232 Google Scholar


D. Razansky et al., “Multispectral opto-acoustic tomography of deep-seated fluorescent proteins in vivo,” Nat. Photonics, 3 412 –417 (2009). NPAHBY 1749-4885 Google Scholar


P. Shao, B. Cox and R. J. Zemp, “Estimating optical absorption, scattering, and Grueneisen distributions with multiple-illumination photoacoustic tomography,” Appl. Opt., 50 3145 –3154 (2011). APOPAI 0003-6935 Google Scholar


Z. Yuan and H. Jiang, “Simultaneous recovery of tissue physiological and acoustic properties and the criteria for wavelength selection in multispectral photoacoustic tomography,” Opt. Lett., 34 1714 –1716 (2009). OPLEDP 0146-9592 Google Scholar


T. Ding, K. Ren and S. Vallelian, “A one-step reconstruction algorithm for quantitative photoacoustic imaging,” Inverse Prob., 31 095005 (2015). INPEEY 0266-5611 Google Scholar


A. Pulkkinen et al., “Direct estimation of optical parameters from photoacoustic time series in quantitative photoacoustic tomography,” IEEE Trans. Med. Imag., 35 2497 –2508 (2016). ITMID4 0278-0062 Google Scholar


K. Ren and F. Triki, “A global stability estimate for the photoacoustic inverse problem in layered media,” (2017). Google Scholar


P. Burgholzer et al., “Exact and approximative imaging methods for photoacoustic tomography using an arbitrary detection surface,” Phys. Rev. E, 75 046706 (2007). Google Scholar


B. T. Cox, S. R. Arridge and P. C. Beard, “Photoacoustic tomography with a limited-aperture planar sensor and a reverberant cavity,” Inverse Prob., 23 S95 –S112 (2007). INPEEY 0266-5611 Google Scholar


D. Finch, M. Haltmeier and Rakesh, “Inversion of spherical means and the wave equation in even dimensions,” SIAM J. Appl. Math., 68 392 –412 (2007). SMJMAP 0036-1399 Google Scholar


M. Haltmeier, T. Schuster and O. Scherzer, “Filtered backprojection for thermoacoustic computed tomography in spherical geometry,” Math. Methods Appl. Sci., 28 1919 –1937 (2005). Google Scholar


M. Haltmeier and G. Zangerl, “Spatial resolution in photoacoustic tomography: effects of detector size and detector bandwidth,” Inverse Prob., 26 125002 (2010). INPEEY 0266-5611 Google Scholar


Y. Hristova, “Time reversal in thermoacoustic tomography-an error estimate,” Inverse Prob., 25 055008 (2009). INPEEY 0266-5611 Google Scholar


B. Huang et al., “Improving limited-view photoacoustic tomography with an acoustic reflector,” J. Biomed. Opt., 18 110505 (2013). JBOPFO 1083-3668 Google Scholar


L. Kunyansky, “Thermoacoustic tomography with detectors on an open curve: an efficient reconstruction algorithm,” Inverse Prob., 24 055021 (2008). INPEEY 0266-5611 Google Scholar


L. V. Nguyen, “A family of inversion formulas in thermoacoustic tomography,” Inverse Probl. Imaging, 3 649 –675 (2009). Google Scholar


G. Paltauf et al., “Photoacoustic tomography with integrating area and line detectors,” Photoacoustic Imaging and Spectroscopy, 251 –263 CRC Press, Boca Raton, Florida (2009). Google Scholar


J. Qian et al., “An efficient Neumann-series based algorithm for thermoacoustic and photoacoustic tomography with variable sound speed,” SIAM J. Imaging Sci., 4 850 –883 (2011). Google Scholar


R. W. Schoonover and M. A. Anastasio, “Compensation of shear waves in photoacoustic tomography with layered acoustic media,” J. Opt. Soc. Am. A, 28 2091 –2099 (2011). JOAOD6 0740-3232 Google Scholar


B. E. Treeby, “Acoustic attenuation compensation in photoacoustic tomography using time-variant filtering,” J. Biomed. Opt., 18 036008 (2013). JBOPFO 1083-3668 Google Scholar


A. Björck, Numerical Methods for Least Squares Problems, SIAM, Philadelphia, Pennsylvania (1996). Google Scholar


H. Ammari et al., “Mathematical modelling in photo-acoustic imaging of small absorbers,” SIAM Rev., 52 677 –695 (2010). SIREAD 0036-1445 Google Scholar


G. Bal, K. Ren, “Non-uniqueness result for a hybrid inverse problem,” Tomography and Inverse Transport Theory, 559 29 –38 American Mathematical Society, Providence, Rhode Island (2011). Google Scholar


A. V. Mamonov and K. Ren, “Quantitative photoacoustic imaging in radiative transport regime,” Comm. Math. Sci., 12 201 –234 (2014). Google Scholar


J.-P. Bérenger, “A perfectly matched layer for the absorption of electromagnetic waves,” J Comput. Phys., 114 185 –200 (1994). Google Scholar


F. Collino and C. Tsogka, “Application of the perfectly matched absorbing layer model to the linear elastodynamic problem in anisotropic heterogeneous media,” Geophysics, 66 294 –307 (2001). GPYSA7 0016-8033 Google Scholar

Biographies for the authors are not available.

© 2018 Society of Photo-Optical Instrumentation Engineers (SPIE) 1083-3668/2018/$25.00 © 2018 SPIE
Patrick Bardsley, Kui Ren, and Rongting Zhang "Quantitative photoacoustic imaging of two-photon absorption," Journal of Biomedical Optics 23(1), 016002 (2 January 2018).
Received: 5 July 2017; Accepted: 5 December 2017; Published: 2 January 2018

Back to Top