1.1.

Photoacoustic Imaging

Breast cancer is the second leading cause of cancer-related death for American women,¹ and early detection and accurate diagnosis are critical for reducing its mortality rate.² The current standard of care includes mammography screening as well as adjunct supplemental ultrasound for diagnosing suspicious lesions.³ However, mammography requires ionizing radiation and has low sensitivity in dense breast tissue while ultrasound has low specificity for breast cancer.^1–4 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 (${\mathrm{HbO}}_2$, Hb), PAI is capable of visualizing vasculature and, through multispectral measurements of tissue absorption, mapping tissue blood oxygen saturation (${\mathrm{SO}}_2$).^7,8 As some hallmarks of malignant cancer include tissue hypoxia and abnormal vasculature,⁹ PAI may potentially improve cancer detection by providing information on tissue function to complement structural information provided by other breast imaging modalities.^10–13 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.^13–20 PAI device performance assessment currently relies on combinations of bench testing (including tissue-mimicking phantoms),^21–25 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.

1.2.

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,³⁰ 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).^31–33 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,34–41 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,⁴⁵ 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)⁵² 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.^{38–40,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.

1.3.

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.

2.1.

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 model^28,53,54 with k-Wave^48,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, $P_0$, is then computed based on the well-known relationship

Fig. 1

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

2.2.

Model Verification and Validation

2.2.1.

Verification against literature

To verify the accuracy of our model, results from a previous modeling study using a similar modeling approach⁴² 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 ${\mu }_{a\_\text{breast}}=0.03{\mathrm{mm}}^{-1}$, ${\mu }_s^{\prime }=0.5{\mathrm{mm}}^{-1}$, $g=0.9$, $n=1.42$). Three different cyst contrast scenarios were considered: cyst absorption lower than the background ($0.0015{\mathrm{mm}}^{-1}$, negative contrast), equal to the background ($0.03{\mathrm{mm}}^{-1}$, zero contrast), or higher than the background ($0.09{\mathrm{mm}}^{-1}$, 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).

2.3.

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\mathrm{mJ}/{\mathrm{cm}}^2$ over all simulation cases. The MC grid/domain was set to $45\times 45\times 45\mathrm{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.

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, $\alpha$, was chosen to mimic breast tissue with $\alpha =0.75f^N\mathrm{dB}/\mathrm{cm}$, where $f$ is the frequency and $N=1.5$.⁵⁷ We then measured axial and lateral spatial resolution in resulting images to determine how image quality varies with transducer frequency response.

3.1.

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.

3.2.

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.²⁵ Computed lateral resolution values for simulated images are similar to those from experimental images (DAS and Verasonics algorithms) as well as our previous reports.^23–25 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.²⁵ 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.

3.3.

Parametric Study

3.3.1.

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\mathrm{mm}\times 2.5\mathrm{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 $<\sim 1.3\mathrm{dB}$ on average for the largest beam cases (12.6 mm circle, $2\mathrm{mm}\times 20\mathrm{mm}$ ellipse, and $4\mathrm{mm}\times 40\mathrm{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\mathrm{mJ}/{\mathrm{cm}}^2$, 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 $\sim 0.2$ to $5.0{\mathrm{cm}}^2$, 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\mathrm{mm}\times 20\mathrm{mm}$ and $4\mathrm{mm}\times 40\mathrm{mm}$) beams, the lateral intensity uniformity becomes less depth-dependent, with on average $<4.35\mathrm{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\mathrm{mm}\times 2.5\mathrm{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.⁶² 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\mathrm{mm}\times 2.5\mathrm{mm}$ and $4\mathrm{mm}\times 40\mathrm{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\mathrm{mJ}/{\mathrm{cm}}^2$. Energy deposition colorbars in $\mathrm{mJ}/{\mathrm{cm}}^3$, 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\mathrm{mm}\times 2.5\mathrm{mm}$, $1\mathrm{mm}\times 10\mathrm{mm}$, $2\mathrm{mm}\times 20\mathrm{mm}$, and $4\mathrm{mm}\times 40\mathrm{mm}$.

3.3.2.

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.²⁵

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.⁶³ 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.