Multidomain computational modeling of photoacoustic imaging: verification, validation, and image quality prediction

Abstract. 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. [1][2][3][4] Photoacoustic imaging (PAI) is a rapidly emerging modality combining pulsed optical excitation and acoustic detection for deep mapping of lightabsorbing chromophores to depths of several centimeters. 5,6 Because optical absorption in tissue is dominated by oxy-and deoxy-hemoglobin (HbO 2 , Hb), PAI is capable of visualizing vasculature and, through multispectral measurements of tissue absorption, mapping tissue blood oxygen saturation (SO 2 ). 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. [10][11][12][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][14][15][16][17][18][19][20] PAI device performance assessment currently relies on combinations of bench testing (including tissue-mimicking phantoms), [21][22][23][24][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.

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). [31][32][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, 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. [38][39][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.

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.

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 E Q -T A R G E T ; t e m p : i n t r a l i n k -; e 0 0 1 ; 3 2 6 ; 3 1 6 P 0 ðrÞ ¼ ΓðrÞμ a ðrÞΦðrÞ; (1) wherer 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 wavelengthdependent 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 radiofrequency (RF) signals then were used to reconstruct in silico photoacoustic images using a delay-and-sum (DAS) beamforming algorithm with dynamic aperture.

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 Journal of Biomedical Optics 121910-2 December 2019 • Vol. 24 (12) 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) x − z 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 E Q -T A R G E T ; t e m p : i n t r a l i n k -; e 0 0 3 ; 6 3 ; 2 4 6 E Q -T A R G E T ; t e m p : i n t r a l i n k -; e 0 0 4 ; 6 3 ; 2 whereũ, p, ρ, ρ 0 , and c 0 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.75 f 1.5 dB∕cm∕MHz 1.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.

Verification against literature
To verify the accuracy of our model, results from a previous modeling study using a similar modeling approach 42 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 mm −1 , μ 0 s ¼ 0.5 mm −1 , g ¼ 0.9, n ¼ 1.42). Three different cyst contrast scenarios were considered: cyst absorption lower than the background (0.0015 mm −1 , negative contrast), equal to the background (0.03 mm −1 , zero contrast), or higher than the background (0.09 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 Journal of Biomedical Optics 121910-3 December 2019 • Vol. 24 (12) 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.

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∕cm 2 , 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 method 60 (see Table 1). Filaments were assumed to have an absorption coefficient of 4 cm −1 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 cm −1 . Penetration depth and target contrast were quantified by computing tube contrast ratio (CR) and signal-to-noise ratio (SNR) as 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 Journal of Biomedical Optics 121910-4 December 2019 • Vol. 24 (12) 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∕cm 2 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.
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.75 f N 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.

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. Table 1 Simulation tissue optical properties at 800 nm (a. low-absorbing background; b. medium-absorbing background).

Background
Inclusions   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.

Model Validation
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  Journal of Biomedical Optics 121910-6 December 2019 • Vol. 24 (12) 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. [23][24][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. 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.
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.

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∕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 ∼0.2 to 5.0 cm 2 , the required energy varied from 2 to 50 mJ, which would require high-energy sources such as secondharmonic-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 depthdependent, 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)     Journal of Biomedical Optics 121910-8 December 2019 • Vol. 24 (12) 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. 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 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.

Conclusion
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 propagationwhich will enable simulation of out-of-plane generated signalsand 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.

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