2.1.

Forward Model

The complete forward model combines the model of light propagation through the scattering medium and the model of photoacoustic wave propagation. The photon diffusion equation is used to describe the fluence distribution inside the medium. The simple analytical expressions for fluence inside the target and background medium²³ are implemented in the model. The fluence inside the target, ${\mathrm{\Phi }}_{\text{in}}$, and background medium, ${\mathrm{\Phi }}_{\text{out}}$, are only related to background optical properties, target optical properties, and target position as:

In Eqs. (1a) and (1b), $A_{l,m}$, $B_{l,m}$, $C_{l,m}$, and $D_{l,m}$ are coefficients. $j_l(x)$ and $n_l(x)$ are the spherical Bessel and Neumann functions, respectively, and $Y_{l,m}(\theta ,\varphi )$ is the spherical harmonics. $k_{\mathrm{out}}$ and $k_{\mathrm{in}}$ are the complex wave numbers of the background medium and target, respectively, and are related to their optical properties. As a result, the absorbed optical energy distribution ($Q$) can be determined and related to the background optical properties, target optical properties, and target position as well. The initial photoacoustic pressure wave induced by a short-pulsed laser is proportional to the absorbed optical energy density distribution and Gruneisen coefficient ($\mathrm{\Gamma }$). The initial pressure, $p_0(r,{\mu }_a)$, and the subsequent acoustic wave, $p(r,t)$, can be expressed by the following equations:²⁴

2.2.

Fitting Procedure

The complete fitting procedure, described by a three-step process, is illustrated in Fig. 1. The first step is to obtain the background optical properties using the diffuse optical measurement system.²⁵ Since a probe with a flat surface was used for housing the ultrasound transducer, optical sources, and detectors, a semi-infinite geometry was a reasonable approximation for the optical measurements. The following set of linear equations were used to relate the calibrated amplitude $I_{\alpha \beta }$ and phase ${\varphi }_{\alpha \beta }$ measured at the source-detector pair $\alpha \beta$, to the source-detector separation ${\rho }_{\alpha \beta }$ as,²⁶

Fig. 1

Fitting procedure to quantitatively recover the target optical properties.

The optical properties of the target can be fitted using an iterative method to minimize the object function. The nonlinear regression function fminsearch, used in the minimization procedure, is available in MATLAB. If the optical properties of the background are unknown, then a four-parameter fitting procedure (${\mu }_a$ and ${\mu }_s^{\prime }$ of the background medium, and ${\mu }_a$ and ${\mu }_s^{\prime }$ of the target) is needed. In this case, the two-parameter fitting procedure can be modified to a four-parameter one as follows: (a) the first step in the three-step-procedure is skipped; (b) the target location information is determined in the second step as before; and (c) the object function in the third step is modified to incorporate ${\mu }_a$ and ${\mu }_s^{\prime }$ of the background medium. This modified procedure hence allows the parameters, ${\mu }_a$ and ${\mu }_s^{\prime }$ of the background medium, and ${\mu }_a$ and ${\mu }_s^{\prime }$ of the target, to be fitted.

2.3.

Experimental System and Phantom

Figure 2(a) shows the experimental configuration of the diffuse optical measurements-assisted (DOM) PAT system. A Ti:Sapphire (Symphotics TII, LS-2134) laser optically pumped by a $Q$-switched Nd:YAG laser (Symphotics-TII, LS-2122), delivered 10-ns pulses at 15 Hz and 750 nm. The laser output was coupled into two 1000-micron step-index input fibers (Thorlabs, BFL48-1000) using a convex lens and circular beamsplitter arrangement as shown in the figure. The convex lens, which has 99% transmittance at 750-nm wavelength, has a focal length of 20 cm and focuses the light into the pair of fibers that are placed at its focal point. The beam splitter was used to split the incoming beam into two beams, one for each fiber. It was measured to have a 60% transmitting and 36% reflecting split ratio for horizontally polarized light at the 750 nm that was used for the experiment. The overall coupling efficiency of the setup (including the losses in the lens and beam splitter) was about 85%, and the total output energy was about 18 mJ per pulse. The two output ends of the fibers were then mounted on opposite (longitudinal) sides of a low-frequency linear ultrasound transducer described below and used for photoacoustic imaging.

Fig. 2

(a) Experimental configuration of diffuse-optical-measurement-assisted PAT system; (b) probe geometry used for simulation and phantom experiments.

A low-frequency ultrasound transducer was employed to maximize the sensitivity for photoacoustic signal detection. The ultrasound transducer manufactured by Vermon (France) had 64 elements with 0.85-mm pitch. The center frequency of the transducer was 1.3 MHz with a 6-dB response from 500 to 1800 kHz. An integrated acoustic lens with 2.5-cm focal length increased the sensitivity at imaging depths for which the optical fluences were low. The photoacoustic signals from the 64-array elements were individually amplified by 60 dB using receiver electronics designed and constructed by our group and multiplexed into eight data acquisition channels with a sampling frequency of 40 MHz and 12-bit precision.¹³ Due to the multiplexing, eight laser firings were required to generate a single 64-channel capture. Because the laser had 15-Hz repetition rate, the acquisition rate was about two frames per second.

The diffuse optical measurement system was a frequency-domain system and consisted of light source and detection subsystems, as well as a hand-held probe.²⁵ The source had four laser diodes of wavelengths 740, 780, 808, and 830 nm. Each laser diode was modulated at 140 MHz and sequentially coupled into six optical source fibers on the probe through optical switches. On the receiving end, the reflected light from the turbid medium associated with each of the illumination sources was detected simultaneously with 12 optical fiber bundles of 3 mm diameter. The detection fibers then coupled the light into 12 parallel photomultiplier tubes (PMTs) and electronics channels. An aluminum foil-covered plastic plate with dimensions of $9\times 10\mathrm{cm}$ was used for the hand-held probe as shown in Fig. 2(b). The middle slot was used for the ultrasound transducer, and the light source and detector fibers were deployed around the probe. The reflective aluminum foil prevented light from being absorbed by the probe and generating image artifacts.²⁸ Semi-infinite partial reflection boundary condition was used for the estimation of the background optical properties, and the aluminum cover provided an effective refraction coefficient of 0.6. In this study, the background optical properties obtained from the 780-nm laser diode only were used.

In the phantom experiments, 0.6% Intralipid solution was used to emulate the background breast tissue. The calibrated optical absorption coefficient (${\mu }_a$) and reduced scattering coefficient (${\mu }_s^{\prime }$) of the Intralipid solution were in the range of 0.01 to 0.02 and 4.0 to $8.0{\mathrm{cm}}^{-1}$, respectively, for the experiments. A $1{\mathrm{cm}}^{-3}$-cubical soft absorber of higher (${\mu }_a=0.28{\mathrm{cm}}^{-1}$), and another of lower (${\mu }_a=0.08{\mathrm{cm}}^{-1}$) optical contrast was made from polyvinyl chloride plastisol and used to mimic a malignant and benign lesion respectively. The calibration was performed on a larger piece phantom of dimensions of $10\times 10\times 5{\mathrm{cm}}^3$ using the DOM system. In general, a malignant lesion has higher contrast with optical absorption coefficient ranging from 0.18 to $0.4{\mathrm{cm}}^{-1}$, and a benign lesion has the lower contrast with absorption coefficient from 0.007 to $\text{0.15}{\mathrm{cm}}^{-1}$.^29–31

3.1.

Simulation Results with/without Background Optical Properties

Photoacoustic signals with 2% zero-mean Gaussian noise of a target of 1 cm diameter was simulated. In order to mimic breast legions of different contrasts, the simulation was performed for different ${\mu }_a$ of the target ranging from 0.30 to $0.05{\mathrm{cm}}^{-1}$ in $0.05{\mathrm{cm}}^{-1}$ decrements and ${\mu }_s^{\prime }=7.0{\mathrm{cm}}^{-1}$. The target was submerged in a scattering medium, having optical properties of ${\mu }_a=0.02{\mathrm{cm}}^{-1}$ and ${\mu }_s^{\prime }=7.0{\mathrm{cm}}^{-1}$, and was located at the position ($x=0$, $y=0$, $z=2.0$) cm relative to the center of the probe. If one assumes that the target location and size are already known from information obtained from the photoacoustic image, then there remain four unknowns: ${\mu }_a$ and ${\mu }_s^{\prime }$ of the target, and ${\mu }_a$ and ${\mu }_s^{\prime }$ of the background medium. Figure 3(a) and 3(b) are the results from the four-parameter fitting. The solid lines in Fig. 3(a) represent the true ${\mu }_a$ value of the target (blue line) and background (black line), and the red stars and green squares are the fitted ${\mu }_a$ values for the target and background, respectively. The $x$-axis represents the expected target ${\mu }_a$ and background ${\mu }_a$ in ${\mathrm{cm}}^{-1}$, and $y$-axis is the fitted ${\mu }_a$ in ${\mathrm{cm}}^{-1}$. When the target has low contrast, such as $0.05{\mathrm{cm}}^{-1}$, the error is largest: 102% for the fitted target ${\mu }_a$, and 399% for the fitted background ${\mu }_a$. On the other hand, as the target contrast gets higher, the errors gradually decrease to a minimum and then increase again as the contrast increases. The errors in the fitted target ${\mu }_a$ were 102%, 34%, 3.1%, 16.5%, 28.4%, and 36% corresponding to the target ${\mu }_a$ of 0.05, 0.10, 0.15, 0.20, 0.25, and $0.30{\mathrm{cm}}^{-1}$, respectively. The fitted background ${\mu }_a$ decreased from 0.099 to $0.022{\mathrm{cm}}^{-1}$ when the target contrast increased from 0.05 to $0.30{\mathrm{cm}}^{-1}$ with the error ranging from 400% to 7%. Clearly, four-parameter fitting is not accurate enough to determine the target ${\mu }_a$. This is because many unknowns can make the nonlinear fitting algorithm settle at the local minimum of the object function. Figure 3(b) shows the fitted results for the target ${\mu }_s^{\prime }$ and background ${\mu }_s^{\prime }$ with the error ranging from 12% to 27% for target ${\mu }_s^{\prime }$, and from 20% to 70% for background ${\mu }_s^{\prime }$. The fitted results for the case where the background ${\mu }_a$ and ${\mu }_s^{\prime }$ are already known from DOM are shown in Fig. 3(c) and 3(d). It is seen from Fig. 3(c) that the errors for the fitted target ${\mu }_a$ have decreased to 4.6%, 1.2%, 0.7%, 0.55%, 0.5%, and 0.3% corresponding to the true target ${\mu }_a$ of 0.05, 0.10, 0.15, 0.20, 0.25, and $0.30{\mathrm{cm}}^{-1}$, respectively. The errors for the fitted target ${\mu }_s^{\prime }$ also ranged from 27% to 36%. Overall, the fitting accuracy for both the high and low contrast targets can be greatly improved if the total number of unknown parameters is decreased. Besides, the run-time of the fitting algorithm is also greatly reduced. This improvement was realized by using background optical properties obtained from the diffused optical measurements to reduce the four-parameter fitting to a two-parameter one. Figure 3(b) and 3(d) shows that the fitted target ${\mu }_s^{\prime }$ is not very sensitive to the change in the background optical properties. As a result, the target ${\mu }_a$ is the only fitted parameter in the phantom experiments.

Fig. 3

Fitting results from the simulation data with 2% noise: (a) fitted background ${\mu }_a$ and target ${\mu }_a$ with four unknown parameters; (b) fitted background ${\mu }_s^{\prime }$ and target ${\mu }_s^{\prime }$ with four unknown parameters; (c) fitted target ${\mu }_a$ with two unknown parameters; (d) fitted target ${\mu }_s^{\prime }$ with two unknown parameters.

3.2.

Phantom Experiments

The optical properties of the background medium can be calculated using Eqs. (3a) and (3b). Using the frequency domain diffuse optical measurement system described in Sec. 2.3, a set of amplitudes ($I_{\alpha \beta }$) and phases (${\varphi }_{\alpha \beta }$) is measured at different distances (${\rho }_{\alpha \beta }$), where $\alpha$ represents the source number and $\beta$ represents the detector number. Because the parameters for each PMT and individual laser diode varied considerably from one to another, we had to calibrate the gain and phase shift for each channel. Considering the gain difference and phase shift, Eq. (3) can be modified as

Fig. 4

Measured amplitude and phase profiles of diffuse waves: (a) $\log ({\rho }_{\alpha \beta }^2{I}_{\alpha \beta })$ versus source-detector distances; (b) phases versus source-detector distances.

In the phantom experiments, the following calibration process was introduced to scale the simulation data to the measurements from the U.S. transducer array: we started with a phantom object with known optical properties, which was a $1\text{-}{\mathrm{cm}}^3$ cubical soft absorber made of polyvinyl chloride plastisol with the calibrated optical properties, ${\mu }_a=0.30{\mathrm{cm}}^{-1}$ and ${\mu }_s^{\prime }=6.2{\mathrm{cm}}^{-1}$. An empirical relation, shown in Eq. (7), was established from PAT simulation trials of the phantom. The equation considered the scan angle of each transducer element in the simulation, and weighted the simulated pressure according to the distance of each pixel from the transducer elements as

Fig. 5

Comparison of the experimental and simulated signals: (a) experimental signal and simulated signal obtained from one U.S. transducer element; (b) the total intensities obtained from experimental signals and simulation signals from 64 channels.

A $1{\mathrm{cm}}^3$-cubical absorber of higher (${\mu }_a=0.28{\mathrm{cm}}^{-1}$) and another of lower (${\mu }_a=0.08{\mathrm{cm}}^{-1}$) optical contrast were used as soft tissue targets to mimic a high-contrast malignant lesion and low-contrast benign lesion, respectively. The absorber was placed inside Intralipid with calibrated optical properties of ${\mu }_a=0.019{\mathrm{cm}}^{-1}$ and ${\mu }_s^{\prime }=5.24{\mathrm{cm}}^{-1}$, and imaged at depths of 1.5, 2.0, 2.5, and 3.0 cm measured from the transducer face to the center of the absorber. Figure 6(a) shows the PAT images of the higher-contrast target with its center located at 2.0 cm from the transducer. The horizontal axis of the PAT images is the lateral dimension (in mm) in the $x$-direction, and the vertical axis is the depth (in mm) in the $z$-direction. The high-resolution PAT images clearly delineate the front and back face of the target, as shown in the circled regions in the figure. The location of the target can be read as (0.2, 2.0) cm in $x$ and $z$ directions. Because the photoacoustic images provided only a two-dimensional (2-D) image, the center of the target in the $y$ dimension was set to $y=0$, which is the center of the transducer [as shown in Fig. 2(b)]. The fitted target ${\mu }_a$ converged to $0.255{\mathrm{cm}}^{-1}$ after 18 iterations, which was 91% of the true value. Table 1 lists the fitted results and errors for the higher-contrast absorber at different depths. The errors in the fitted data range from 7.1% to 16.8%; these are 18% to 37% better than the reconstructed DOT results with guidance from PAT image in Ref. 14. In the referenced paper, the maximum reconstructed absorption coefficients for the higher-contrast target have errors about 25% to 54%.

Fig. 6

(a) PAT images of the high contrast target with center located at 2.0 cm depth; (b) PAT images of the low contrast target with center located at 2.0 cm depth. The x-axis of the PAT images is lateral dimension and vertical $z$-axis is the depth. The units are mm. The circled region shows the cubic target.

Table 1

Phantom experiment for single high-contrast target (calibrated μa=0.28  cm−1) and low-contrast target (calibrated μa=0.08  cm−1).

Similar experiments were performed for the lower-contrast target. Figure 6(b) shows the PAT images of this target whose center location may be read as $(-0.5,2.0)\mathrm{cm}$. The contrast of the PAT images was not as good as that of the higher-contrast one. However, the front and back faces of the target were fairly resolved. After 14 iterations, the fitted target ${\mu }_a$ had a final converged value of $0.086{\mathrm{cm}}^{-1}$, which is 108% of the true value. The fitting results at other depths are listed in Table 1, and the errors at different depths range from 3.8% to 11.2%. These values are either comparable or 20% better than the DOT results for low contrast target in Ref. 14.

3.3.

Sensitivity Test

The sensitivity of the two-parameter fitting procedure to the target parameters and background optical properties was evaluated. We found that accurate target information, namely the location ($x$, $y$, $z$) and size, are critical for the fitting procedure. As shown in Fig. 7(a), where the $x$-axis is the measured radius of target and $y$-axis is the fitted ${\mu }_a$ in ${\mathrm{cm}}^{-1}$, if the target’s radius is not estimated accurately, varying from 0.4 to 0.6 cm, the fitted ${\mu }_a$ is greatly influenced by this error. This is easy to understand because the total signal is related to the volume of the target. When the estimated target size is smaller than that of true size, the fitted absorption coefficient increases to match the measured signal; in this case, the errors increase from 15% to 74%. An example of this is shown in Fig. 7(a) for a target centered at 2.5 and 3.0 cm, and in general is true for all depths. The fitted ${\mu }_a$ is also strongly related to the target depth as shown in Fig. 7(b). In this figure, the $x$-axis is the assumed target depth, and we still used the data in which the target was centered at 2.5 and 3.0 cm depths. When the target depth input is not accurate, the fitted error also increases from 15% to 100%. On the other hand, the $x\text{-y}$ location of target has small influence on the fitting accuracy as shown in Fig. 7(c). The $x$-axis represents the target location in the $x$ direction, and $y$-axis represents the fitted target ${\mu }_a$ values. When the location of the target was moving in the $x$-direction away from the true location up to 0.2 cm, the fitted ${\mu }_a$ varied less than 7%.

Fig. 7

The influence of structural information: (a) radius; (b) depth; (c) $x\text{-}y$ plane.

The background optical properties are also critical to the two-parameter fitting procedure. Figure 8(a) shows that the target fitted ${\mu }_a$ changes as the input background ${\mu }_a$ increases. In the figure, the $x$-axis represents the background ${\mu }_a$ in ${\mathrm{cm}}^{-1}$, and $y$-axis is the fitted ${\mu }_a$. The error ranges from 16% to 45%. The background ${\mu }_s^{\prime }$ also has similar influence on the fitting procedure as shown in Fig. 8(b). In this figure, $x$-axis represents the background ${\mu }_s^{\prime }$, and $y$-axis is the fitted target ${\mu }_a$. When the background ${\mu }_s^{\prime }$ is not accurate, the fitting error also increases from 16% to 60%. On the other hand, the target ${\mu }_s^{\prime }$ has a small influence on the fitting accuracy as shown in Fig. 8(c). The $x$-axis represents the target ${\mu }_s^{\prime }$ and ranges from 1 to 11 in $2{\mathrm{cm}}^{-1}$ increments, and $y$-axis represents the fitted target ${\mu }_a$. When the target ${\mu }_s^{\prime }$ was changed, the fitted ${\mu }_a$ varied less than 20%.

Fig. 8

The influence of optical properties: (a) background ${\mu }_a$ change; (b) background ${\mu }_s^{\prime }$ change; (c) target ${\mu }_s^{\prime }$ change.

3.4.

Blood Tube Experiments

To test the robustness of the fitting procedure, we also performed experiments with blood collected from a sacrificed mouse. The blood was diluted and injected into a thin tube with 0.58-mm inner diameter and 0.9-mm outer diameter. To estimate the absorption coefficient of the blood, the formula from Ref. 33 was used.

Fig. 9

PAT image of the blood tube.