Translator Disclaimer
8 November 2019 Multidomain computational modeling of photoacoustic imaging: verification, validation, and image quality prediction
Author Affiliations +

As photoacoustic imaging (PAI) technology matures, computational modeling will increasingly represent a critical tool for facilitating clinical translation through predictive simulation of real-world performance under a wide range of device and biological conditions. While modeling currently offers a rapid, inexpensive tool for device development and prediction of fundamental image quality metrics (e.g., spatial resolution and contrast ratio), rigorous verification and validation will be required of models used to provide regulatory-grade data that effectively complements and/or replaces in vivo testing. To address methods for establishing model credibility, we developed an integrated computational model of PAI by coupling a previously developed three-dimensional Monte Carlo model of tissue light transport with a two-dimensional (2D) acoustic wave propagation model implemented in the well-known k-Wave toolbox. We then evaluated ability of the model to predict basic image quality metrics by applying standardized verification and validation principles for computational models. The model was verified against published simulation data and validated against phantom experiments using a custom PAI system. Furthermore, we used the model to conduct a parametric study of optical and acoustic design parameters. Results suggest that computationally economical 2D acoustic models can adequately predict spatial resolution, but metrics such as signal-to-noise ratio and penetration depth were difficult to replicate due to challenges in modeling strong clutter observed in experimental images. Parametric studies provided quantitative insight into complex relationships between transducer characteristics and image quality as well as optimal selection of optical beam geometry to ensure adequate image uniformity. Multidomain PAI simulation tools provide high-quality tools to aid device development and prediction of real-world performance, but further work is needed to improve model fidelity, especially in reproducing image noise and clutter.




Photoacoustic Imaging

Breast cancer is the second leading cause of cancer-related death for American women,1 and early detection and accurate diagnosis are critical for reducing its mortality rate.2 The current standard of care includes mammography screening as well as adjunct supplemental ultrasound for diagnosing suspicious lesions.3 However, mammography requires ionizing radiation and has low sensitivity in dense breast tissue while ultrasound has low specificity for breast cancer.14 Photoacoustic imaging (PAI) is a rapidly emerging modality combining pulsed optical excitation and acoustic detection for deep mapping of light-absorbing chromophores to depths of several centimeters.5,6 Because optical absorption in tissue is dominated by oxy- and deoxy-hemoglobin (HbO2, Hb), PAI is capable of visualizing vasculature and, through multispectral measurements of tissue absorption, mapping tissue blood oxygen saturation (SO2).7,8 As some hallmarks of malignant cancer include tissue hypoxia and abnormal vasculature,9 PAI may potentially improve cancer detection by providing information on tissue function to complement structural information provided by other breast imaging modalities.1013 A wide variety of PAI system designs are described in the literature, and differences in system design parameters such as illumination geometry, optical wavelength(s), transducer array geometry, and transducer acoustic frequency response are expected to produce systems with different levels of performance.1320 PAI device performance assessment currently relies on combinations of bench testing (including tissue-mimicking phantoms),2125 preclinical animal studies,26,27 and clinical trials,17,18 each of which possesses varying degrees of burden, cost, and ability to adequately represent in vivo imaging scenarios in complex tissues. Therefore, there is a need for additional nonclinical tools to more rapidly and accurately predict real-world performance of PAI devices.


Multidomain Simulation of PAI

Computational modeling is an increasingly common approach for PAI device evaluation and performance assessment, including simulation of optical, acoustic, and thermal transport phenomena.28,29 Models with sufficient accuracy and real-world fidelity can be used to predict device performance and thus supplement or eventually replace costlier phantom, animal, and/or clinical testing. Fast, cheap prediction of fundamental photoacoustic image quality metrics often measured using phantom tests (e.g., spatial resolution, penetration depth) is especially appealing as such metrics can inform iterative optimization of device design. However, it is unclear how to evaluate photoacoustic computational model performance or establish model credibility for a given intended application. Due to increased use of computational models in medical device development and evaluation,30 a standardized approach for assessing model credibility through verification and validation (V&V) procedures has been developed through the American Society of Mechanical Engineers (ASME).3133 Within this framework, verification involves determining that a computational model accurately solves the underlying set of equations. Validation involves determining how accurately a model represents a real-world system through comparison to experimental data.

Various approaches for multidomain PAI simulations using different optical and acoustic components have been proposed. Here, we provide a brief review of proposed models used for PAI and discuss observed limitations in their performance assessment. Monte Carlo (MC) simulation is the gold standard method for modeling light–tissue interactions in turbid media and has been previously applied to PAI. MC has been used to compare performances of different PAI device designs,20,3441 to evaluate target lesion visualization and detectability,39,42 and to enable quantitative PAI.43,44 Common tools for modeling acoustic wave propagation in tissue include Field II,45 which has been used to simulate photoacoustic response and quantify spatial resolution of a proposed PAI system,46,47 and k-Wave,48,49 a popular open-source photoacoustic simulation toolbox for MATLAB used by several groups to study PAI systems.42,43

For PAI simulation, several groups have proposed multidomain finite element models based on commercial software (e.g., COMSOL)50,51 or open-source packages (e.g., ONELAB)52 to simulate photoacoustic processes by explicitly modeling heat transfer, solid mechanics, and acoustic wave propagation. While these models show promising results, their greater complexity due to a larger number of parameters and modeling assumptions implies greater difficulty in model verification and validation. These studies also relied on the diffusion approximation of light transport, rather than using gold-standard MC simulations. Several studies have demonstrated combining MC simulations with various acoustic propagation models,29,38,40 and several groups have coupled MC with the k-space forward acoustic model implemented in k-Wave.42,43 Our literature review of previously proposed integrated photoacoustic models showed that a majority of studies did not validate model outputs against experiments. Some studies compared simulated output images to input images based on numerical phantoms, but such activities fall under verification, as validation against experiments is an explicit requirement of the ASME V&V framework and is essential for model accuracy evaluation. In addition, several studies only performed qualitative verification of simulated images, rather than providing quantitative assessment. Lastly, while each individual modeling component (e.g., MC or k-Wave) may be well validated in a general sense, within the ASME framework this does not necessarily imply sufficient V&V of the entire modeling chain for a particular “context of use,” such as predicting photoacoustic image quality metrics in silico. Therefore, despite growing interest in PAI modeling, there is apparent lack of consensus regarding model assessment methodology and the level of rigor needed to establish usefulness of such models.

PAI models have been previously applied toward prediction of PAI device performance.3840,50,51 In addition to deficiencies in verification and validation described above, studies that performed quantitative image quality assessment often characterized only one image quality metric (e.g., spatial resolution) or only considered detectability of a single target, rather than an array of targets throughout the field of view. Evaluation of multiple targets is important as several image quality metrics have been shown to vary within the image.25,36 Furthermore, the impact of several key design parameters, including optical beam geometry and transducer frequency response, on photoacoustic image quality has not been adequately evaluated through parametric study using a well-validated model. Thus, it is unclear whether proposed models are truly suitable for predicting fundamental image quality metrics of PAI devices.


Study Objectives

Our goals were to apply ASME V&V principles toward rigorous performance assessment of an integrated PAI model representative of approaches seen in the literature and to evaluate the model’s capability for predicting real-world device performance in terms of fundamental image quality metrics. To this end, our aims were to: (1) couple a previously developed three-dimensional (3D) MC model with k-Wave, (2) verify the model by comparing outputs against published simulation data, (3) validate the model against experimental images acquired with a custom PAI system, and (4) demonstrate model utility through parametric study of how key system design parameters influence image quality.


Materials and Methods


Computational Model of PAI

To develop a multidomain PAI model representative of common approaches from the literature, we coupled our previously developed 3D voxel-based MC model28,53,54 with k-Wave48,49 (Fig. 1). Tissue optical properties and inclusion geometry are defined in the MC model, which enables computation of the spatial fluence and energy deposition distributions. The initial acoustic pressure distribution, P0, is then computed based on the well-known relationship

Eq. (1)

where r is the position, Γ(r) is the Gruneisen parameter describing the conversion efficiency of light to sound, μa(r) is the absorption coefficient, and Φ(r) is the optical fluence. Γ(r) is known to vary with tissue type and temperature,42 while fluence distribution changes with tissue type and wavelength-dependent optical absorption and scattering coefficients. In this study, Γ(r) was assumed to be homogeneous at a fixed value of 0.2.18 k-Wave was used to solve coupled first-order acoustic equations using a k-space pseudospectral time-domain approach to compute wave propagation in tissue. The simulated radio-frequency (RF) signals then were used to reconstruct in silico photoacoustic images using a delay-and-sum (DAS) beamforming algorithm with dynamic aperture.

Fig. 1

PAI model flowchart. The small 2D images represent the output of each step.



Monte Carlo simulation

Light transport was simulated using our previously developed and validated 3D voxel-based MC model.28,53,54 Briefly, the model is based on a 3D voxel grid capable of simulating heterogeneous geometries such as multi-layer structures containing cylindrical inclusions, with tissue optical properties (absorption coefficient μa, scattering coefficient μs, anisotropy factor g, and refractive index n) assigned to each voxel based on labeled tissue type. Based on convergence analysis, for all simulations except the verification study (grid size matched to the published data) and the validation of spatial resolution (grid size of 0.1 mm), we used a voxel size of 0.15 mm, with 200 million launched photons. All MC simulations were performed using a supercomputing cluster comprising 94 computing nodes with 8 CPUs and 24 GB of RAM per node.55 A two-dimensional (2D) xz plane in the computed 3D energy deposition map at 9  mm offset from the beam center in y direction was extracted to calculate initial pressure within the image plane. The offset plane was used to mimic the real-world experimental setups, in which the ultrasound transducer is offset from the optical source.17,20,56


Acoustic simulation

One-way acoustic wave propagation was simulated using k-Wave,48,49 which uses a k-space pseudospectral time-domain solution of the coupled first-order wave equations describing conservation of momentum [Eq. (2)], conservation of mass [Eq. (3)], and the pressure–density relation [Eq. (4)]:

Eq. (2)


Eq. (3)


Eq. (4)

where u, p, ρ, ρ0, and c0 corresponding to acoustic particle velocity, acoustic pressure, acoustic density, ambient density, and isotropic sound speed, respectively. k-Wave offers broad capabilities ideal for simulating acoustic physical processes during PAI, including simulation of acoustically heterogeneous media, arbitrary source geometry, and simulated acoustic signal measurement by discrete transducers (or multiple-element transducer arrays) with specified geometry and frequency response. k-Wave model inputs include: (1) discretized medium geometry, (2) medium acoustic properties, (3) initial acoustic sources, and (4) acoustic detector array geometry and frequency response. In this study, the frequency response of the transducers was implemented as a frequency domain Gaussian filter utilizing center frequency and fractional bandwidth of the simulated transducer using k-Wave features. Based on convergence analysis, simulations were created using an axial grid size of one-fourth of the transducer acoustic wavelength and a lateral grid size one-third of the array element pitch. Tissue was assumed to be homogeneous with speed of sound equal to 1480 m/s. Tissue acoustic attenuation was assumed to be negligible for model verification and validation studies while parametric studies assumed attenuation similar to breast tissue (α=0.75f1.5  dB/cm/MHz1.5 where f is the frequency).57 All acoustic simulations were performed in 2D using a personal computer with i7-5600 CPU processor at 2.6 GHz processing frequency and 16 GB of RAM.


Model Verification and Validation


Verification against literature

To verify the accuracy of our model, results from a previous modeling study using a similar modeling approach42 were replicated and compared. This model was selected as the authors provided sufficiently detailed descriptions of optical and acoustic model parameters as well as phantom geometries to enable replication of reported model outputs. Heijblom et al. modeled a 1064 nm, transmission-mode PAI system to investigate visualization of a benign breast cyst within a breast mimicking layer (with μa_breast=0.03  mm1, μs=0.5  mm1, g=0.9, n=1.42). Three different cyst contrast scenarios were considered: cyst absorption lower than the background (0.0015  mm1, negative contrast), equal to the background (0.03  mm1, zero contrast), or higher than the background (0.09  mm1, positive contrast). For the purposes of verifying our model outputs for these scenarios, all optical properties, acoustic properties, laser beam characteristics, and phantom geometry [phantom 1, Fig. 2(a)] were matched. While Heijblom et al. performed 3D acoustic simulations and used a 3D backprojection algorithm for reconstruction, we performed 2D acoustic simulations and employed a DAS with dynamic aperture algorithm for reconstruction. Furthermore, Heijblom et al. modeled the cyst as a 10-mm diameter at the depth of 13 mm (from surface to the upper boundary of the cylinder) and reported that cylindrical and spherical targets produced only minor differences in received signals. We thus modeled the cyst as a 10-mm-diameter cylinder at a depth of 13 mm instead of a 10-mm-diameter sphere as our MC model already supports cylinders, and we are generally interested in cylindrical target geometries. Our simulation results were compared quantitatively with the results of Heijblom et al. by computing the percent difference between computed energy deposition profiles, decay rate of energy deposition, and cyst thickness from reconstructed photoacoustic images.

Fig. 2

Phantom geometries used for simulations. (a) Tissue-mimicking phantom with embedded cyst structure of diameter 10 mm (phantom 1). (b) Filament array of diameter 0.05 mm (phantom 2). (c) Penetration depth phantom (phantom 3) with embedded targets of 0.5 mm. (d) 2D array phantom with embedded targets of 0.5 mm (phantom 4).



Experimental validation in phantoms

To validate our PAI model, model outputs were compared to experimental image data acquired in physical phantoms (phantom 2 and phantom 3) with well-characterized optical and acoustic properties using a custom PAI system allowing for switching of ultrasound transducers as described elsewhere.24,25 The system consisted of a tunable near-infrared optical parametric oscillator (Phocus Mobile, Opotek, Inc., Carlsbad, California) capable of emitting pulses with a 5-ns duration over 690 to 950 nm at a rate of 10 Hz and a research-grade ultrasound acquisition system with 128 channels (Vantage 128, Verasonics, Inc., Kirkland, Washington). An elliptical beam (4×40  mm) was aligned directly adjacent to the transducer for in-plane illumination with transducer in contact with the imaging medium. Phantoms were imaged at 800 nm with a radiant exposure of 10  mJ/cm2, which is 32% of the ANSI maximum permissible exposure.58 Images were acquired with an 8 MHz, 128 element linear array with 109% bandwidth and 0.3 mm pitch (L11-4v, Verasonics).

The transducer was gel-coupled with a protective aluminum foil layer to reduce unwanted photoacoustic signal generation at the tissue surface. Thirty frames of experimental data were collected for each validation study and then averaged for postprocessing and analysis.

Experimental and simulated data were both used to reconstruct photoacoustic images. To approximate the background noise/clutter, 2D uniform random noise was generated and then scaled by a depth-dependent exponential curve empirically derived from fitted experimental background RF signals. The generated background noise/clutter map was then added to each simulated RF data channel to approximate 2D random noise in resultant simulated images. A theoretical approach to simulating speckle noise has also been proposed.59 The beamforming processing was a DAS algorithm with dynamic aperture assuming a homogeneous speed of sound of 1480 m/s. Image magnitude was computed using the Hilbert transform. Background images were subtracted from both simulated and experimental images pixel wise. Experimental images were background-subtracted by subtracting the average of 30 image frames acquired with no laser emission. Simulated background was generated as the simulated image of a phantom containing no embedded targets. Background-subtracted images were normalized to the greatest target inclusion intensity and log compressed for display.

Two 7×7×5  cm phantoms comprised of diluted Intralipid (I141, Sigma-Aldrich, Inc.) and India ink (3080-F, Chartpak, Inc., Leeds, Massachusetts) were imaged to provide experimental data for model validation. First, a spatial resolution phantom was constructed containing a seven-by-five grid of 50-μm-diameter black nylon monofilament suture wires with 5-mm horizontal spacing and 7.5-mm vertical spacing [phantom 2, Fig. 2(b)]. Phantom optical properties were measured using integrating sphere spectrophotometry and the inverse adding-doubling method60 (see Table 1). Filaments were assumed to have an absorption coefficient of 4  cm1 to produce high-contrast line sources. Axial and lateral spatial resolution were computed as the full-width at half-maximum of the linear (not log compressed) profile of each target intensity. Second, a penetration depth phantom was constructed containing a diagonal array of 0.5-mm-inner diameter polytetrafluoroethylene (PTFE) tubes at depths from 5 to 35 mm with 5-mm vertical and horizontal spacing (STT-24, Component Supply Company, Fort Meade, Florida). Tubes were filled with India ink solutions with an optical absorption coefficient of 4  cm1.

Table 1

Simulation tissue optical properties at 800 nm (a. low-absorbing background; b. medium-absorbing background).

μa (cm−1)μs′ (cm−1)gnμa (cm−1)μs′ (cm−1)gn
Validation (filament array phantom)(a) 0.02100.91.334.00.0010.631.47 61
(b) 0.1
Validation (penetration array phantom)(a) 0.02100.91.334.00.0010.631.4761
(b) 0.1
Parametric study (2D vessel array phantom)280.1100.61.45.0130.9871.4

Penetration depth and target contrast were quantified by computing tube contrast ratio (CR) and signal-to-noise ratio (SNR) as

Eq. (5)


Eq. (6)

where I is the mean target ROI intensity, B is the mean background ROI intensity, and σ is the background ROI standard deviation. Target ROIs (1.8×2.5  mm) were analyzed by removing pixels with intensity <50% of maximum intensity, and local background ROIs of equal size were selected at the same depths of each target ROI. Penetration depth was computed as the interpolated depth for which SNR reached 6 dB, a threshold that was identified through qualitative image assessment and is consistent with the threshold reported for this system in our previous work.25 However, as detectability is a function of SNR and target size, the 6-dB threshold may only be appropriate for small size target evaluated in this study.


Parametric Study

We used the model to conduct a parametric study of how key system optical design parameters (optical illumination geometry) and acoustic design parameters (detector frequency and bandwidth) impact PAI device performance. To investigate the effect of different optical beam geometries, we used a numerical phantom (phantom 4) containing a 2D grid of 0.5-mm diameter cylindrical targets at depths from 5 to 35 mm with 10-mm vertical and 7.5-mm horizontal spacing. We considered five circular and elliptical beams of varying size, with each circle and ellipse matched approximately in beam area (Table 2). Small differences in beam area of <3% are due to limitations in MC grid resolution. The total energy in each beam was adjusted to ensure a constant radiant exposure of 10  mJ/cm2 over all simulation cases. The MC grid/domain was set to 45×45×45  mm to minimize computation time; we verified that extending the bottom and vertical simulation boundaries to the dimensions of validation phantoms did not significantly affect energy deposition distributions. The lateral uniformity was investigated for different cases by computing the ratio between the intensity of the center targets to periphery. Further, the target lateral detectability versus depth was characterized using the CR of the targets.

Table 2

Optical beam dimensions used for parametric study.

Circular beam
Circle 1Circle 2Circle 3Circle 4Circle 5
Radius (mm)
Area (mm2)2.018.0432.16128.67498.75
Elliptical beam
Ellipse 1Ellipse 2Ellipse 3Ellipse 4Ellipse 5
Minor and major radii (mm×mm)0.25×2.50.5×51×102×204×40
Area (mm2)1.967.8531.41125.66502.65

To further evaluate the effect of ultrasound detector center frequency and bandwidth on the system performance, we simulated imaging of the numerical filaments phantom [phantom 2, Fig. 2(b)] using the L11-4v transducer but varying center frequency (2.5, 5.0, 7.5, and 10 MHz) and bandwidth (50% and 100%). Tissue acoustic attenuation, α, was chosen to mimic breast tissue with α=0.75fN  dB/cm, where f is the frequency and N=1.5.57 We then measured axial and lateral spatial resolution in resulting images to determine how image quality varies with transducer frequency response.


Results and Discussion


Model Verification

Outputs of our model, including depth profiles of energy deposition and RF amplitude through the cyst center, were compared to those reported by Heijblom et al. for negative, zero, and positive cyst contrast conditions (Fig. 3). Energy deposition profiles were normalized to area under the curve to facilitate visualization, while RF amplitude profiles were normalized to maximum amplitude in the negative cyst contrast case. Qualitative inspection shows consistent energy deposition and RF amplitude trends with depth between our model and data reported by Heijblom et al., as well as visualization of cyst upper and lower boundaries. There are small discrepancies in terms of RF signal trace shape and precise cyst boundary location that may be attributed to our use of 2D acoustic simulations as opposed to 3D, as well as a different choice of photoacoustic image reconstruction algorithm. Quantitatively, the difference in signal amplitude between the models was within 9% to 22% at different depth segments, which is sufficient considering the level of noise apparent in the published data. Further, our model predicted an energy distribution decay rate of 0.27/mm averaged across all three different absorption conditions, compared to the 0.20 to 0.22/mm range reported by Heijblom et al. Discrepancies between quantitative metrics from each model may be due to the higher number of photons used in our simulation which led to smoother, more convergent energy deposition profiles. Overall, verification results indicate that our model adequately reproduces outputs consistent with a published model despite some differences in model components and implementation.

Fig. 3

Comparison of energy deposition (top row) and acoustic RF profiles (bottom row) with depth from our model outputs (black solid lines) and those of Heijblom et al. (blue dashed lines) for negative, zero, and positive cyst contrast. The black dashed lines denote upper and lower cyst boundaries.



Model Validation

Figure 4 shows photoacoustic images reconstructed from simulated and experimental RF data in the filament array phantom (phantom 2), as well as computed axial and lateral resolution. Target appearance is similar in simulated and experimental images, and spatial resolution values computed from simulated DAS images were consistent with those computed from experimental images, including images processed using a proprietary Verasonics reconstruction algorithm rather than our DAS algorithm (images not shown). These results are consistent with our previous work, which used the same transducer but employed only the Verasonics algorithm.25 Computed lateral resolution values for simulated images are similar to those from experimental images (DAS and Verasonics algorithms) as well as our previous reports.2325 Axial and lateral resolution were found to slightly vary with depth in simulated images, possibly due to the use of nonoptimized dynamic aperture focusing in deeper tissue regions. Discrepancies in lateral resolution between experimental DAS and Verasonics images may also be due to unknown processing steps applied to the RF data by the Verasonics acquisition system. These results indicate that our model can predict image quality metrics such as spatial resolution in sufficient agreement with experimental results to suggest predictive capability.

Fig. 4

Reconstructed photoacoustic images of filament phantom (phantom 2) for (a) simulated and (b) experimental RF data and (c) computed axial and (d) lateral resolution from simulated and experimental data. The color bar is in dB.


Figure 5 shows simulated and experimental reconstructed photoacoustic images of the penetration depth phantom (phantom 3). Qualitative inspection of these images indicates general agreement, but experimental images present clutter and reflection artifacts not present in simulated images. The differences between simulated and experimental data can be observed from line plots across the second target (Fig. 5, lower row). From this, it is clear that experimental data include clutter and artifact absent in the simulated data. These artifacts are caused by diffusely reflected light from the phantom being absorbed at the transducer surface, an effect we have observed with our system previously.25 Simulation of this effect would require knowledge of the optical properties and geometry of the aluminum foil shielding, which has an irregular, irreproducible shape as a fresh piece must be applied in each experiment.

Fig. 5

Upper row: Reconstructed photoacoustic images from penetration depth phantom (phantom 3) for (a) and (b) low-absorbing and (c) and (d) medium-absorbing background, using (a) and (c) experimental and (b)–(d) simulated data. Data are normalized to the intensity of the shallowest target intensity. The color bar is in dB. Lower row: line plot across second target (white line in a) for depth of 5 to 20 mm.


Simulated and experimental images produced consistent trends in CR and SNR as functions of depth (Fig. 6). An SNR threshold of 6 dB was found to be correlated with detectability limits based on qualitative inspection. Table 3 shows the computed penetration depth using interpolated target SNR for simulation and experiment. Simulated images predicted a slightly greater penetration depth for the case of low-absorption case with penetration depth up to 33.2 mm compared to 27.7 mm computed for experimental images. This may be due to the absence of clutter and limited fidelity of background noise injected into the simulation data, or due to discrepancies in modeling transducer element sensitivity and directivity. However, penetrations depths were very close for the medium-absorbing background. Further model improvements are needed to adequately reproduce image characteristics that influence contrast-based image quality metrics.

Fig. 6

Comparison of (a) target CR and (b) SNR in simulated and experimental images. Here, low abs. and med. abs refer to low absorption and medium absorption phantoms, respectively.


Table 3

Maximum penetration depth for simulated and experimental images based on SNR threshold of 6 dB.

Low-absorbing backgroundMedium-absorbing background
SNR threshold-based limit (mm)33.227.726.926.9


Parametric Study


Optical beam geometry

In addition to predicting experimental system performance, we evaluated the utility of our model for investigating the effect of optical beam geometry on system performance. As shown in Fig. 7, energy deposition in targets strongly depends on spatial location, beam geometry, and beam size, with greater deposited energy for larger beam sizes. Since radiant exposure was held constant, this trend is likely caused by enhanced fluence under large beams due to multiple scattering in turbid media. We also observed that small circular and elliptical beams result in poor lateral intensity uniformity, with targets in the center of the field of view appearing with much higher intensity than those on the periphery. For small beam cases (0.8- to 6.4-mm diameter circles, 0.25  mm×2.5  mm ellipse), the range in intensity between the center and peripheral targets in each row was up to an order of magnitude (20 dB), while it was <1.3  dB on average for the largest beam cases (12.6 mm circle, 2  mm×20  mm ellipse, and 4  mm×40  mm ellipse). On the other hand, lateral uniformity was less depth-dependent for the large circular and elliptical beams, which is likely due to increased fluence as a result to multiple scattering typical of large-area beams. Furthermore, elliptical beams generally provide superior lateral intensity uniformity to circular beams of equal area. However, as all beam cases used the same radiant exposure of 10  mJ/cm2, a greater energy per pulse would be required to achieve the improved performance observed with larger beams. This an important consideration in selecting an optical source during PAI device design. Since the beam areas studied varied from 0.2 to 5.0  cm2, the required energy varied from 2 to 50 mJ, which would require high-energy sources such as second-harmonic-pumped OPOs. As shown in Fig. 8(a), lateral intensity uniformity was found to substantially improve with depth for all but the largest circular and elliptical beams, which is expected due to light diffusion in a turbid medium. For larger circular (12.6 mm) and elliptical (2  mm×20  mm and 4  mm×40  mm) beams, the lateral intensity uniformity becomes less depth-dependent, with on average <4.35  dB variation down to 25 mm [Fig. 8(a)]. This is likely due to the increased similarity of optical beam size with the transducer dimensions, providing more uniform diffuse illumination above deep, peripheral targets and thus increasing local energy deposition in those regions. Figures 8(b) and 8(c) show the computed CR as a function of depth for different beam size cases. It can be observed that smaller beams (0.8 to 6.4 mm diameter circles and 0.25  mm×2.5  mm ellipses) result in large difference in CR values (about 15 dB) for center and peripheral targets as a function of depth, which leads to poor peripheral target detectability. On the other hand, using larger beam sizes, the CR values of center and peripheral targets become closer (with maximum difference of about 3 dB) resulting in better peripheral target visibility as a function of depth. It is worth noting that a recently released commercial-grade PAI system (LAZR-X, VisualSonics, Inc.) allows users to swap in linear array fiberoptic bundles of different lengths to optimize light delivery for specific imaging applications.62 This indicates that design of optimal beam geometry for a given imaging scenario is perhaps more complex than anticipated, and our results highlight the critical importance of beam geometry as a design parameter affecting PAI device performance.

Fig. 7

Energy deposition maps and corresponding simulated photoacoustic images for (a) and (b) 0.8- and 12.6-mm circular beams and (c) and (d) elliptical beams of size 0.25  mm×2.5  mm and 4  mm×40  mm. The small lower-right figure in each energy deposition map is an en face view of beam fluence at the phantom surface, which were self-normalized for visualization purposes. All beam cases used a fixed uniform radiant exposure of 10  mJ/cm2. Energy deposition colorbars in mJ/cm3, photoacoustic image colorbars in dB.


Fig. 8

(a) Lateral intensity variation in dB. These values were computed as a ratio of the center target intensity to the one on most lateral position. (b) CR of center (dashed line) and periphery targets using circular and (c) elliptical beams. Due to overlap between the first and second beam size cases, the second beam size for both circular and elliptical was eliminated for clarity. Circles 1, 3, 4, and 5 correspond to circular beams with 0.8, 3.2, 6.4, and 12.6 mm diameter. Ellipses 1, 3, 4, and 5 refer to elliptical beams of size 0.25  mm×2.5  mm, 1  mm×10  mm, 2  mm×20  mm, and 4  mm×40  mm.



Acoustic detector parameters

Figure 9 shows the reconstructed photoacoustic images from simulated data, using four different center frequencies (columns) and two different bandwidths (rows) for the ultrasound transducer. Since medium acoustic attenuation generally increases with frequency, penetration depth decreases with center frequency. Lateral resolution improves as center frequency increases while higher bandwidth (Fig. 9, bottom row) results in worse lateral resolution as opposed to the lower bandwidth (Fig. 9, top row). Figure 10 shows computed axial and lateral resolution from simulated data as functions of depth for detectors with different center frequency and bandwidth. Axial resolution [Fig. 10(a)] varies approximately linearly with depth, which supports our validation and previous results.25

Fig. 9

Reconstructed photoacoustic images of filament phantom (phantom 2) using ultrasound transducer arrays with varying center frequency (columns) as well as fractional bandwidth of 50% (top row) and 100% (bottom row). Each image was normalized to its maximum target intensity.


Fig. 10

(a) Computed axial and (b) lateral resolution using ultrasound detector of four different center frequencies with 50% and 100% bandwidth (solid versus dashed lines).


Results demonstrate that axial resolution is improved when using higher center frequencies, with diminishing returns for higher frequencies (7.5 to 10 MHz), but fractional bandwidth has a significant effect on both axial and lateral resolution (Fig. 10). Lateral resolution was also found to improve with increasing center frequency but decrease with increasing bandwidth. This effect is caused by center frequency downshifting in attenuating media, which is more pronounced with greater signal bandwidth. Axial and lateral resolution results are both consistent with those reported in a previous study of plane-wave ultrasound imaging.63 These results demonstrate the effect of ultrasound detector characteristics on system performance and show the potential of the proposed model for improving understanding of PAI device design consequences for particular imaging applications.



We have developed a multidomain model of PAI physical processes and conducted model verification and validation to ensure modeling results are representative of real-world PAI system performance. Further, we have investigated the impact of key system design consideration on image quality using standard metrics (e.g., spatial resolution, CR, and SNR) by conducting parametric studies. Our optical parametric study results indicated that for fixed radiant exposure, large, elliptical beams offer better peripheral target detectability deeper in the tissue to circular beams or short elliptical beams. Our acoustic parametric study demonstrated that both axial and lateral resolution improve using higher center frequencies, while lateral resolution was found to improve with smaller bandwidth. In spite of these accomplishments, there is room to improve the realism of our model in the future by incorporating 3-D acoustic propagation—which will enable simulation of out-of-plane generated signals—and more accurately representing the true transducer impulse response. These steps and others will allow models to better simulate image quality and reliably predict the impact of noise and clutter. By improving our understanding of complex optical and acoustic physical processes, modeling tools can facilitate device design optimization and performance evaluation of PAI systems.


The authors have no conflicts of interest, financial or otherwise, to disclose. The mention of commercial products, their sources, or their use in connection with material reported herein is not to be construed as either an actual or implied endorsement of such products by the U.S. Department of Health and Human Services.


We gratefully acknowledge funding support from the FDA Medical Countermeasures Initiative (MCM259 and MCM282) and the FDA Office of Women’s Health. We would also like to acknowledge the FDA high-performance computing team for their assistance.



W. A. Berg et al., “Ultrasound as the primary screening test for breast cancer: analysis from ACRIN 6666,” J. Natl. Cancer Inst., 108 (4), (2015). JNCIEQ Google Scholar


D. A. Berry et al., “Effect of screening and adjuvant therapy on mortality from breast cancer,” N. Engl. J. Med., 353 (17), 1784 –1792 (2005). NEJMAG 0028-4793 Google Scholar


R. Shah, K. Rosso and S. D. Nathanson, “Pathogenesis, prevention, diagnosis and treatment of breast cancer,” World J. Clin. Oncol., 5 (3), 283 –298 (2014). Google Scholar


W. Xia, W. Steenbergen and S. Manohar, “Photoacoustic mammography: prospects and promises,” Breast Cancer Manag., 3 (5), 387 –390 (2014). Google Scholar


L. V. Wang and S. Hu, “Photoacoustic tomography: in vivo imaging from organelles to organs,” Science, 335 (6075), 1458 –1462 (2012). SCIEAS 0036-8075 Google Scholar


J. L. Su et al., “Advances in clinical and biomedical applications of photoacoustic imaging,” Expert Opin. Med. Diagn., 4 (6), 497 –510 (2010). Google Scholar


M. Xu and L. V. Wang, “Photoacoustic imaging in biomedicine,” Rev. Sci. Instrum., 77 (4), 041101 (2006). RSINAK 0034-6748 Google Scholar


M. Toi et al., “Visualization of tumor-related blood vessels in human breast by photoacoustic imaging system with a hemispherical detector array,” Sci. Rep., 7 41970 (2017). SRCEC3 2045-2322 Google Scholar


J. C. Forster et al., “A review of the development of tumor vasculature and its effects on the tumor microenvironment,” Hypoxia, 5 21 –32 (2017). Google Scholar


K. S. Valluru and J. K. Willmann, “Clinical photoacoustic imaging of cancer,” Ultrasonography, 35 (4), 267 –280 (2016). Google Scholar


M. Mehrmohammadi et al., “Photoacoustic imaging for cancer detection and staging,” Curr. Mol. Imaging, 2 (1), 89 –105 (2013). Google Scholar


M. Heijblom et al., “Photoacoustic image patterns of breast carcinoma and comparisons with magnetic resonance imaging and vascular stained histopathology,” Sci. Rep., 5 11778 (2015). SRCEC3 2045-2322 Google Scholar


E. I. Neuschler et al., “A pivotal study of optoacoustic imaging to diagnose benign and malignant breast masses: a new evaluation tool for radiologists,” Radiology, 287 (2), 398 –412 (2017). RADLAX 0033-8419 Google Scholar


M. Pramanik et al., “Design and evaluation of a novel breast cancer detection system combining both thermoacoustic (TA) and photoacoustic (PA) tomography,” Med. Phys., 35 (6), 2218 –2223 (2008). MPHYA6 0094-2405 Google Scholar


S. Manohar et al., “The twente photoacoustic mammoscope: system overview and performance,” Phys. Med. Biol., 50 (11), 2543 –2557 (2005). PHMBA7 0031-9155 Google Scholar


X. Li et al., “High resolution functional photoacoustic tomography of breast cancer,” Med. Phys., 42 (9), 5321 –5328 (2015). MPHYA6 0094-2405 Google Scholar


A. A. Oraevsky et al., “Clinical optoacoustic imaging combined with ultrasound for coregistered functional and anatomical mapping of breast tumors,” Photoacoustics, 12 30 –45 (2018). Google Scholar


E. Fakhrejahani et al., “Clinical report on the first prototype of a photoacoustic tomography system with dual illumination for breast cancer imaging,” PLoS One, 10 (10), e0139113 (2015). POLNCL 1932-6203 Google Scholar


P. K. Upputuri and M. Pramanik, “Recent advances toward preclinical and clinical translation of photoacoustic tomography: a review,” J. Biomed. Opt., 22 (4), 041006 (2016). JBOPFO 1083-3668 Google Scholar


M. Li et al., “Linear array-based real-time photoacoustic imaging system with a compact coaxial excitation handheld probe for noninvasive sentinel lymph node mapping,” Biomed. Opt. Express, 9 (4), 1408 –1422 (2018). BOEICL 2156-7085 Google Scholar


M. Fonseca et al., “Characterisation of a phantom for multiwavelength quantitative photoacoustic imaging,” Phys. Med. Biol., 61 (13), 4950 –4973 (2016). PHMBA7 0031-9155 Google Scholar


J. R. Cook, R. R. Bouchard and S. Y. Emelianov, “Tissue-mimicking phantoms for photoacoustic and ultrasonic imaging,” Biomed. Opt. Express, 2 (11), 3193 –3206 (2011). BOEICL 2156-7085 Google Scholar


C. Jia et al., “Two-layer heterogeneous breast phantom for photoacoustic imaging,” J. Biomed. Opt., 22 (10), 106011 (2017). JBOPFO 1083-3668 Google Scholar


W. C. Vogt et al., “Biologically relevant photoacoustic imaging phantoms with tunable optical and acoustic properties,” J. Biomed. Opt., 21 (10), 101405 (2016). JBOPFO 1083-3668 Google Scholar


W. C. Vogt et al., “Phantom-based image quality test methods for photoacoustic imaging systems,” J. Biomed. Opt., 22 (9), 095002 (2017). JBOPFO 1083-3668 Google Scholar


K. E. Wilson et al., “Multiparametric spectroscopic photoacoustic imaging of breast cancer development in a transgenic mouse model,” Theranostics, 4 (11), 1062 –1071 (2014). Google Scholar


A. Hariri et al., “The characterization of an economic and portable LED-based photoacoustic imaging system to facilitate molecular imaging,” Photoacoustics, 9 10 –20 (2018). Google Scholar


T. Gould, Q. Wang and T. J. Pfefer, “Optical-thermal light-tissue interactions during photoacoustic breast imaging,” Biomed. Opt. Express, 5 (3), 832 –847 (2014). BOEICL 2156-7085 Google Scholar


S. L. Jacques, “Coupling 3D Monte Carlo light transport in optically heterogeneous tissues to photoacoustic signal generation,” Photoacoustics, 2 (4), 137 –142 (2014). Google Scholar


T. M. Morrison et al., “Advancing regulatory science with computational modeling for medical devices at the FDA’s Office of Science and Engineering Laboratories,” Front. Med., 5 241 (2018). FMBEEQ Google Scholar


ASME V&V 40–2018, Assessing Credibility of Computational Modeling Through Verification and Validation, American Society of Mechanical Engineers, New York (2018). Google Scholar


ASME V&V 10– 2006, Guide for Verification and Validation in Computational Solid Mechanic, American Society of Mechanical Engineers, New York (2006). Google Scholar


ASME V&V 20– 2009, Standard for Verification and Validation in Computational Fluid Dynamics & Heat Transfer, American Society of Mechanical Engineers, New York (2009). Google Scholar


V. Periyasamy and M. Pramanik, “Advances in Monte Carlo simulation for light propagation in tissue,” IEEE Rev. Biomed. Eng., 10 122 –135 (2017). Google Scholar


K. Sivasubramanian et al., “Optimizing light delivery through fiber bundle in photoacoustic imaging with clinical ultrasound system: Monte Carlo simulation and experimental validation,” J. Biomed. Opt., 22 (4), 41008 (2017). JBOPFO 1083-3668 Google Scholar


G. S. Sangha, N. J. Hale and C. J. Goergen, “Adjustable photoacoustic tomography probe improves light delivery and image quality,” Photoacoustics, 12 6 –13 (2018). Google Scholar


S. H. El-Gohary et al., “Design study on photoacoustic probe to detect prostate cancer using 3D Monte Carlo simulation and finite element method,” Biomed. Eng. Lett., 4 (3), 250 –257 (2014). Google Scholar


G. Paltauf, P. R. Torke and R. Nuster, “Modeling photoacoustic imaging with a scanning focused detector using Monte Carlo simulation of energy deposition,” J. Biomed. Opt., 23 (12), 121607 (2018). JBOPFO 1083-3668 Google Scholar


V. Periyasamy and M. Pramanik, “Monte Carlo simulation of light transport in tissue for optimizing light delivery in photoacoustic imaging of the sentinel lymph node,” J. Biomed. Opt., 18 (10), 106008 (2013). JBOPFO 1083-3668 Google Scholar


Y. Lou et al., “Impact of nonstationary optical illumination on image reconstruction in optoacoustic tomography,” J. Opt. Soc. Am. A, 33 (12), 2333 –2347 (2016). JOAOD6 0740-3232 Google Scholar


Y. Lou, A. Oraevsky and M. A. Anastasio, “Application of signal detection theory to assess optoacoustic imaging systems,” Proc. SPIE, 9708 97083Z (2016). PSISDG 0277-786X Google Scholar


M. Heijblom et al., “Appearance of breast cysts in planar geometry photoacoustic mammography using 1064-nm excitation,” J. Biomed. Opt., 18 (12), 126009 (2013). JBOPFO 1083-3668 Google Scholar


B. A. Kaplan et al., “Monte-Carlo-based inversion scheme for 3D quantitative photoacoustic tomography,” Proc. SPIE, 10064 100645J (2017). PSISDG 0277-786X Google Scholar


J. Laufer et al., “In vitro measurements of absolute blood oxygen saturation using pulsed near-infrared photoacoustic spectroscopy: accuracy and resolution,” Phys. Med. Biol., 50 (18), 4409 –4428 (2005). PHMBA7 0031-9155 Google Scholar


J. A. Jensen, “FIELD: a program for simulating ultrasound systems,” in 10th Nordicbaltic Conf. Biomed. Imaging, 351 –353 (1996). Google Scholar


J. Gamelin et al., “A real-time photoacoustic tomography system for small animals,” Opt. Express, 17 10489 –10498 (2009). OPEXFF 1094-4087 Google Scholar


S. Vaithilingam et al., “Three-dimensional photoacoustic imaging using a two-dimensional CMUT array,” IEEE Trans. Ultrason. Ferroelectr. Freq. Control, 56 (11), 2411 –2419 (2009). ITUCER 0885-3010 Google Scholar


B. T. Cox et al., “k-space propagation models for acoustically heterogeneous media: application to biomedical photoacoustics,” J. Acoust. Soc. Am., 121 (6), 3453 –3464 (2007). JASMAN 0001-4966 Google Scholar


B. E. Treeby and B. T. Cox, “k-Wave: MATLAB toolbox for the simulation and reconstruction of photoacoustic wave fields,” J. Biomed. Opt., 15 (2), 021314 (2010). JBOPFO 1083-3668 Google Scholar


Z. Wang, S. Ha and K. Kim, “A new design of light illumination scheme for deep tissue photoacoustic imaging,” Opt. Express, 20 (20), 22649 –22659 (2012). OPEXFF 1094-4087 Google Scholar


C. Sowmiya and A. K. Thittai, “Simulation of photoacoustic tomography (PAT) system in COMSOL(R) and comparison of two popular reconstruction techniques,” Proc. SPIE, 10137 101371O (2017). PSISDG 0277-786X Google Scholar


C. Fadden and S.-R. Kothapalli, “A single simulation platform for hybrid photoacoustic and RF-acoustic computed tomography,” Appl. Sci., 8 (9), 1568 (2018). Google Scholar


T. J. Pfefer et al., “A three-dimensional modular adaptable grid numerical model for light propagation during laser irradiation of skin tissue,” IEEE J. Sel. Top. Quantum Electron., 2 (4), 934 –942 (1996). IJSQEN 1077-260X Google Scholar


T. J. Pfefer et al., “Modeling laser treatment of port wine stains with a computer-reconstructed biopsy,” Lasers Surg. Med., 24 (2), 151 –166 (1999). LSMEDI 0196-8092 Google Scholar


M. Mikailov et al., “Scaling bioinformatics applications on HPC,” BMC Bioinf., 18 (14), 501 (2017). BBMIC4 1471-2105 Google Scholar


P. K. Upputuri and M. Pramanik, “Photoacoustic imaging in the second near-infrared window: a review,” J. Biomed. Opt., 24 (4), 040901 (2019). JBOPFO 1083-3668 Google Scholar


F. S. Foster and J. W. Hunt, “Transmission of ultrasound beams through human tissue—focusing and attenuation studies,” Ultrasound Med. Biol., 5 (3), 257 –268 (1979). USMBA3 0301-5629 Google Scholar


ANSI Z136.1:2014, American National Standard for Safe Use of Lasers, Laser Institute of America, Orlando, Florida (2014). Google Scholar


E. Hysi et al., “Insights into photoacoustic speckle and applications in tumor characterization,” Photoacoustics, 14 37 –48 (2019). Google Scholar


G. Yoon, D. N. G. Roy and R. C. Straight, “Coherent backscattering in biological media: measurement and estimation of optical properties,” Appl. Opt., 32 (4), 580 –585 (1993). APOPAI 0003-6935 Google Scholar


R. Michels, F. Foschum and A. Kienle, “Optical properties of fat emulsions,” Opt. Express, 16 (8), 5907 –5925 (2008). OPEXFF 1094-4087 Google Scholar


FUJIFILM VisualSonics, Inc., “Vevo LAZR-X multimodal imaging | FUJIFILM VisualSonics,” Google Scholar


Z. Alomari et al., “The effect of the transducer parameters on spatial resolution in plane-wave imaging,” in IEEE Int. Ultrason. Symp. (IUS), 1 –4 (2015). Google Scholar


Nima Akhlaghi received his BS and MS degrees in physical and biomedical electronics from Kyiv Polytechnic Institute, Ukraine, and Sumy State University, Ukraine, in 2009 and 2010. He received his PhD in electrical and computer engineering from George Mason University in 2017. He joined the FDA in 2017 as a visiting scientist where he is conducting ultrasound imaging research in the Office of Science and Engineering Laboratories (OSEL). His research interests include ultrasound and photoacoustic imaging, modeling, and elastography.

T. Joshua Pfefer received his BS degree in mechanical engineering from Northwestern University, his MS degree in mechanical engineering, and his PhD in biomedical engineering from the University of Texas at Austin. He then trained as a research fellow at the Wellman Laboratories of Photomedicine. In 2000, he joined the FDA, where he is the leader of the Optical Diagnostic Devices Laboratory. In 2018, he was named a fellow of SPIE for contributions to biophotonics.

Keith A. Wear received his BA degree in applied physics from the University of California, San Diego and his MS and PhD degrees in applied physics from Stanford University. He is a research physicist at the FDA. He is an associate editor of the Journal of the Acoustical Society of America, Ultrasonic Imaging, and IEEE Transactions on Ultrasonics, Ferroelectrics, and Frequency Control. He is a fellow of the Acoustical Society of America, the American Institute for Medical and Biological Engineering, and the American Institute of Ultrasound in Medicine.

Brian S. Garra received his MD from the University of Washington followed by residency in diagnostic radiology at the University of Utah, becoming board certified in 1980. He is a professor of radiology at the University of Vermont and at George Washington University. A fellow of the American Institute of Ultrasound in Medicine and the Society of Radiologists in Ultrasound, he does device evaluation and research in elastography, complex phantom development and global ultrasound outreach at the FDA.

William C. Vogt received his BS degree in mechanical engineering from the University of Massachusetts Amherst in 2009 and his PhD in biomedical engineering from Virginia Polytechnic Institute and State University in 2013. Since joining the FDA in 2013, he has been conducting regulatory science to develop tools and test methods for evaluating the safety and effectiveness of photoacoustic imaging devices. His research interests include photoacoustic imaging, tissue phantoms, nanoparticles, standardization, and biophotonic medical device characterization and evaluation.

© The Authors. Published by SPIE under a Creative Commons Attribution 4.0 Unported License. Distribution or reproduction of this work in whole or in part requires full attribution of the original publication, including its DOI.
Nima Akhlaghi, T. Joshua Pfefer, Keith A. Wear, Brian S. Garra, and William C. Vogt "Multidomain computational modeling of photoacoustic imaging: verification, validation, and image quality prediction," Journal of Biomedical Optics 24(12), 121910 (8 November 2019).
Received: 31 May 2019; Accepted: 14 October 2019; Published: 8 November 2019

Back to Top