Comparison of compressional and elastic wave simulations for patient-specific planning prior to transcranial photoacoustic-guided neurosurgery

Abstract. Significance: Simulations have the potential to be a powerful tool when planning the placement of photoacoustic imaging system components for surgical guidance. While elastic simulations (which include both compressional and shear waves) are expected to more accurately represent the physical transcranial acoustic wave propagation process, these simulations are more time-consuming and memory-intensive than the compressional-wave-only simulations that our group previously used to identify optimal acoustic windows for transcranial photoacoustic imaging. Aim: We present qualitative and quantitative comparisons of compressional and elastic wave simulations to determine which option is more suitable for preoperative surgical planning. Approach: Compressional and elastic photoacoustic k-Wave simulations were performed based on a computed tomography volume of a human cadaver head. Photoacoustic sources were placed in the locations of the internal carotid arteries and likely positions of neurosurgical instrument tips. Transducers received signals from three previously identified optimal acoustic windows (i.e., the ocular, nasal, and temporal regions). Target detectability, image-based target size estimates, and target-to-instrument distances were measured using the generalized contrast-to-noise ratio (gCNR), resolution, and relative source distances, respectively, for each simulation method. Results: The gCNR was equivalent between compressional and elastic simulations. The areas of the −6  dB contours of point spread functions utilized to measure resolution differed by 0.33 to 3.35  mm2. Target-to-instrument distance measurements were within 1.24 mm of the true distances. Conclusions: These results indicate that it is likely sufficient to utilize the less time-consuming, less memory-intensive compressional wave simulations for presurgical planning.

bone and underlying pituitary tumors. 1 Although the procedure is generally safe, 2 morbidity and mortality rates rise to 14% to 23% and 24% to 26%, respectively, if iatrogenic injury to the internal carotid arteries (ICAs) occurs. [3][4][5] Current intraoperative guidance techniques, such as stereotactic guidance and endoscopy, enable monitoring of the ICAs in close proximity to the surgical site; however, they suffer from two primary limitations. First, stereotactic guidance is subject to registration errors, which can become increasingly large as patient anatomy is disrupted during surgery and deviates from the anatomy in preoperative x-ray computed tomography (CT) or magnetic resonance images. Second, endoscopy is unable to identify the ICAs when they are obscured by bone or other tissues in the operative path. Our group is investigating the use of transcranial photoacoustic imaging as an intraoperative imaging technique for realtime visualization of the ICAs to address these two limitations, as detailed in our original research papers [6][7][8][9][10][11] and summarized in the literature surveys on this topic. 12,13 To achieve photoacoustic imaging for guidance of endonasal transsphenoidal surgery, we propose the insertion of light-transmitting fiber optic devices in the nasal cavity, similar to other surgical tools. [6][7][8][9][10][11][12] This light source will then excite the hemoglobin within the ICAs, converting the absorbed optical energy to acoustic energy that is received by an externally placed ultrasound receiver. The external ultrasound receiver placement results in a transcranial photoacoustic imaging scenario, which is challenged by acoustic interactions with bone and is known to degrade image quality. [14][15][16] Previous work from our group developed and demonstrated a simulation method to identify naturally occurring acoustic windows in the adult human skull to minimize acoustic interactions with bone and to provide high-contrast photoacoustic images of the ICAs. 10,11 We demonstrated that patient-specific simulations have the potential to enable preoperative planning to determine appropriate placement of imaging system components. 10,11 Our vision for a presurgical workflow with the use of simulations is shown in Fig. 1, in direct comparison with a workflow without access to presurgical simulations. Without simulations, a surgeon may need to search and find optimal locations to place an ultrasound receiver to best receive the acoustic signals within the skull. This intraoperative process may be prone to the expenditure of time that could instead be spent operating on the patient and would increase the time that the patient is under anesthesia and on the operating table. In addition, the surgeon may have to abandon the use of intraoperative photoacoustic image guidance during the procedure if a viable receiver location is not determined in time, thereby sacrificing the benefits of real-time photoacoustic-based guidance during the surgery. 10 In contrast, the simulations that we propose have the potential to preoperatively identify patient-specific ultrasound receiver locations. 10 The surgeon may then use the simulation results and the patient's preoperative CT scan to construct a time-efficient surgical plan for each operation. With this step completed, the neurosurgery may Fig. 1 Surgical workflow with and without the use of patient-specific simulations. Without photoacoustic simulations, a surgeon will likely have to test multiple ultrasound receiver locations and their feasibility. Simulations have the potential to efficiently identify these locations before the surgical procedure. then be safely executed with the benefit of real-time, intraoperative photoacoustic image guidance.
Our group previously identified the ocular cavity, temporal region, and nasal cavity as three optimal ultrasound receiver locations, using compressional-wave-only simulations. 10,11 However, the shear waves known to propagate within dense media such as bone were excluded from our initial demonstrations. As an alternative, elastic wave simulations, which include both compressional and shear waves, are expected to more accurately represent the physical acoustic process. However, elastic wave simulations, such as the k-Wave elastic simulations based on the classical Kelvin-Voigt absorption model, are time consuming and memory intensive. 17 Time and memory costs reduce the ease and likelihood of clinical translation. This paper presents a comparison of compressional and elastic photoacoustic k-Wave simulations to investigate which simulation type (i.e., compressional or elastic) is required for preoperative surgical planning. The remainder of this paper is organized as follows. Section 2 details the simulation methods and quantitative metrics used for the comparison. Section 3 presents the resulting comparisons. Section 4 discusses the implications of the results with respect to the compressional wave simulations that were used in our experimental validation studies, with regard to the vision shown in Fig. 1 and with an eye toward reducing barriers to clinical translation. Finally, Sec. 5 concludes the paper.

Simulation Configuration
Three-dimensional (3D) photoacoustic k-Wave 17-19 simulations were performed after converting the CT volume of a human cadaver skull into heterogeneous, volumetric maps of corresponding density, compressional and shear wave sound speeds, and compressional and shear wave absorption prefactors. Figure 2 shows an example axial slice from the CT volume of a human cadaver skull. Homogeneous volumes were additionally modeled as the average density, sound speed, and absorption values of brain tissue to provide baseline photoacoustic images that do not contain heterogeneous tissue effects, such as aberration, attenuation, scattering, and reverberation. The simulated tissue properties for the heterogeneous and homogeneous cases are reported in Table 2. A two-dimensional (2D) cross-section of the propagating wave at various time points for the homogeneous and heterogeneous cases is shown in Fig. 3, with more details available in Video 1. The computational grid was defined with a symmetric voxel size of 0.3 mm × 0.3 mm × 0.3 mm. Acoustic wave propagation was simulated with a sampling frequency of 82 MHz (i.e., time increment of 1.12 −8 s) on a NVIDIA Quadro RTX 6000 GPU (Table 1). Fig. 2 Axial slice from the CT volume of the human cadaver skull demonstrating a 2D crosssection of the 3D simulation configuration. Spherical photoacoustic sources were placed within the LCA and at distances of 6 to 13 mm from the LCA to represent the tip of a surgical instrument (a distance of 6 mm is shown). Green lines illustrate locations of independently placed ultrasound transducers.
Phased array ultrasound transducers with 0.3 mm pitch, 13.5 mm height, 0 mm kerf, and 64 elements (as reported in Table 1) were positioned to receive transcranial photoacoustic signals from three acoustic windows: (1) the left ocular cavity, (2) the left temple region, and (3) the nasal cavity, as shown in Fig. 2. A spherical photoacoustic source with a diameter of 0.3 mm (i.e., a point source) or 4 mm (i.e., the diameter of a typical adult carotid artery) was placed in the  location of the left carotid artery (LCA), as shown in Fig. 2. Background absorption was not modeled in the photoacoustic source distributions for the following two reasons. First, background absorption can be mitigated by carefully selecting the illumination wavelength to preferentially excite the oxygenated hemoglobin over surrounding tissues. 28 Second, our previous observations in transcranial photoacoustic imaging of the cadaveric ICAs did not display evidence of background absorption. 10 To investigate photoacoustic image characteristics when the tip of a surgical instrument is also in the imaging plane, a 2-mm-diameter spherical photoacoustic source was positioned within the surgical site at distances of 6 to 13 mm from the 4-mm-diameter LCA target. These distances were measured from the centers of the LCA and instrument targets. Considering this measurement and the best possible system resolution described in Sec. 2.2, the minimum distance simulated was 6 mm. The distances of the LCA source relative to the center of the transducers located on the ocular, nasal, and temple regions were 7.40, 8.02, and 5.03 cm, respectively.

Image Formation and Analysis
Received transducer channel data were bandpass filtered to contain −6 dB frequencies in the range of 1 to 5 MHz. Randomly distributed Gaussian noise was added to the received channel data obtained with the transducer located in the temporal region to model the electronic noise of an imaging system, resulting in a 15-dB channel signal-to-noise ratio (SNR). The same noise distribution was then added to the received channel data obtained with the remaining transducer locations (i.e., nasal and temple regions) to simulate the same noise floor for the three transducer locations, each viewing photoacoustic signals of the same 4 mm target.
Photoacoustic delay-and-sum (DAS) images were generated from the filtered channel data with additive noise. Although advanced beamforming techniques exist to compensate for acoustic heterogeneity, 29-31 the DAS beamformer is a more standard choice that was selected to compare acoustic heterogeneity effects between the two simulation types. DAS image quality (i.e., resolution, target visibility, and target detectability) was measured for each transducer location. Resolution was assessed by calculating the area of the −6 dB contour of the point spread function (PSF), measured from images of the point target. The best possible system resolution was also measured as the distance of the minimum cross-section through the center of the −6 dB contour with the smallest area. Target visibility was assessed using images of the the 4-mmdiameter LCA target to measure contrast, SNR, and contrast-to-noise ratio (CNR), while target detectability was assessed using the generalized contrast-to-noise ratio (gCNR). 32,33 These image quality metrics are defined as follows: 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 ; 1 1 6 ; 9 5 Contrast ¼ 20 log 10  18,[20][21][22][23][24][25] The two shear wave properties of sound speed in brain and absorption power law prefactors in bone and brain are not explicitly known. However, these properties can be estimated as approximately half the corresponding compressional speed of sound and approximately double the corresponding compressional absorption power law prefactor. 24,26 Speed of sound (m/s) Absorption power law prefactor 19 (7) E Q -T A R G E T ; t e m p : i n t r a l i n k -; e 0 0 2 ; 1 1 6 ; 7 2 3 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 ; 1 1 6 ; 6 9 1 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 ; 1 1 6 ; 6 5 4 where μ t and μ b are the means, σ t and σ b are the standard deviations, and h t and h b are the histograms of the signal amplitudes within ellipsoidal regions of interest (ROIs) placed within the photoacoustic target (denoted by subscript t) or within the background of the photoacoustic image (denoted by subscript b), N is the number of bins in the histogram, and k is the index of the bin. A total of N ¼ 145 bins were used to create the histograms for the gCNR measurements. For each image, a singular target ROI was centered on the brightest pixel in the image. For the images obtained with the ocular and nasal transducer locations, six background ROIs were located around the target ROI to calculate means and standard deviations of image quality metrics. For the images obtained with the temple region transducer location, only five background ROIs were used because the sixth was outside the field of view of the transducer. The areas of the target ROI and each background ROI were equivalent. Target-to-instrument tip distances were measured from DAS photoacoustic images of the LCA and the instrument tip, calculated as the Euclidean distance between the centroid and each source. The distance error was calculated as the absolute difference between the measured targetto-instrument tip distance and the known target-to-instrument distance. Table 3 compares the execution time and memory usage of the 3D compressional and elastic simulations, performed with the three independently placed transducers shown in Fig. 2. The last row of this table reports the mean time and memory usage across the three simulations to demonstrate that elastic simulations required ∼4.6× more time and 2.7× more memory than compressional simulations.

Results
Figures 4(a) and 4(b) show simulated photoacoustic channel data from the heterogeneous compressional and elastic simulations of the 0.3-mm-diameter LCA source obtained from the temple acoustic window. Although the initial wavefront shapes across elements look similar between the two simulation types, the compressional simulation has greater reverberations in the photoacoustic signal subsequent to the initial wavefront. To further appreciate the differences in the wavefronts between the simulation types, Fig. 4(c) shows the difference image of the elastic channel data subtracted from the compressional channel data. The wavefront in the compressional simulation has a greater amplitude than the wavefront in the elastic simulation, as indicated by the positive values along the wavefront in the difference image. Specifically, the maximum amplitudes of the wavefronts were 9e −5 and 6e −5 Pa for the compressional and elastic simulations, respectively. Figure 5(a) shows simulated photoacoustic images generated from the heterogeneous simulations of the 4-mm-diameter LCA source. The shape, location, and visibility of the targets generally agree between the compressional and elastic simulation pairs for each transducer location. Figure 5(b) quantifies the target visibility and detectability of the heterogeneous simulated images displayed in Fig. 5(a). There were minimal differences in contrast, SNR, and CNR measurements, and the gCNR measurements were equivalent. Therefore, the visibility and detectability of these photoacoustic image targets ranges from similar to equivalent for compressional and elastic wave simulations. As noted by Kempski et al., 33 there is no benefit to improving  particular methods when the gCNR is already near its maximum value of 1.0, and in this case, the method is the type of simulation used for presurgical planning.
Figures 6(a) and 6(b) show the −6 dB contours of the point target images from the homogeneous and heterogeneous simulations, respectively. The best resolution of the system was 1.4 mm, obtained from the shortest cross-section of the −6 dB contour of the homogeneous compressional simulation image obtained from the temporal region. Figure 6(c) shows the calculated areas of the contours. The contour areas of the baseline homogeneous simulated images were 9.33, 9.80, and 6.60 mm 2 for the ocular, nasal, and temple acoustic windows, respectively, for compressional simulations and 10.27, 11.05, 7.08 mm 2 , respectively, for elastic simulations. The areas of the corresponding simulated images obtained with the heterogeneous skull model were 29.78, 38.26, and 24.22 mm 2 , respectively, for compressional simulations and 28.76, 35.29, and 24.21 mm 2 , respectively, for elastic simulations. Therefore, the resolution of these photoacoustic images is 0.45 to 1.24 mm 2 better for the homogeneous compressional wave simulations than the homogeneous elastic wave simulations. However, in the presence of bone, the elastic wave simulations have 0.33 to 3.35 mm 2 better resolution than compressional wave simulations for the heterogeneous case. Figure 7(a) shows an annotated photoacoustic image from a heterogeneous simulation of the LCA and a surgical instrument tip, obtained with the ocular acoustic window. The ground truth target-to-instrument distance is 11.27 mm, and the measured target-to-instrument distance is 11.18 mm in this image. Figure 7(b) shows image-based measurements of multiple target-to-instrument distances as a function of ground truth distances. The largest distance error between measurements and ground truth was obtained with the nasal cavity images, measuring 1.24 and 1.16 mm for the compressional and elastic wave simulations, respectively.
The error ranges for the ocular and temporal regions were 0.09 to 0.91 mm and 0.14 to 0.88 mm for the compressional and elastic simulations, respectively. When comparing the target-to-instrument distance errors between the compressional and elastic simulations, the differences in errors were 0.03 to 0.09, 0.07 to 0.45, and 0.03 to 0.22 mm for the ocular, nasal, and temporal regions, respectively. Therefore, the relative target-to-instrument distances can be measured with accuracy within the 1.4 mm best possible resolution of the system from both compressional and elastic wave simulations.

Discussion
This paper details our investigations of k-Wave elastic wave simulations as a method for presurgical planning of transducer placement during transcranial photoacoustic-guided neurosurgery. In particular, we assessed comparative differences in wave propagation, image quality (i.e., target visibility, target detectability, and target size estimation), and target-to-instrument tip distances (i.e., distances between a critical blood vessel and surgical instrument tip) using photoacoustic channel data and images generated from compressional and elastic transcranial simulations, resulting in three key observations. First, wavefronts in the channel data were similar in terms of general wavefront shape, but they differed in amplitude and reverberations between the compressional and elastic simulations, as shown in Fig. 4. Second, target visibility and detectability were qualitatively and quantitatively similar between the compressional and elastic simulations for the three acoustic windows tested, as shown in Fig. 5. Third, the accuracy of target-to-instrument tip distances measured using either compressional or elastic simulations was generally similar between elastic and compressional simulations and was within the best possible resolution of the system, as shown in Fig. 7. Although the latter two observations (based on results presented in Figs. 5 and 7) may seem counterintuitive with respect to shear waves being the primary acoustic component during acoustic propagation in bone, we hypothesize that the similarity of target visibility, target detectability, and target-to-instrument distances between the two simulation types is due to the strategic choice of transducer location 10 (i.e., the length of acoustic interaction with bone is minimized along the pathway from the LCA source to transducers) and the relatively minimal thickness of bone with respect to the distance between the source and each transducer.
Our three key observations indicate that it is likely sufficient to utilize the less timeconsuming, less memory-intensive k-Wave compressional simulations for presurgical planning, confirming that appropriate complexity was included when comparing simulation to experimental results from a human cadaver. 10 The arguably minimal 0.33 to 3.35 mm 2 area differences in PSFs between k-Wave compressional and elastic simulations indicate that future investigations surrounding the effects of PSF differences could be useful prior to completely eliminating elastic simulations from the presurgical planning process. In particular, for the three acoustic windows tested, the areas of the −6 dB contours of the PSFs achieved with the compressional wave simulations was up to 1.24 mm 2 better than that of the elastic wave simulations for the homogeneous case and up to 3.35 mm 2 worse than that of the elastic simulations for the heterogeneous Fig. 7 (a) Heterogeneous simulated photoacoustic images of the two sources representing the LCA and instrument tip, obtained with the ocular acoustic window at a relative distance of 11.27 mm. The target-to-instrument distance is the length of the line between the source centers in the image. (b) Comparison of measured target-to-instrument distances as a function of the true target-to-instrument distances for compressional (○) and elastic (Δ) measurements in the ocular, nasal, and temporal acoustic windows. The ideal 1:1 relationship is shown as the dashed black line, with ± 1.4 mm boundaries, representing the best possible resolution of the system, indicated by the dashed gray lines.
case, as shown in Fig. 6. A neurosurgeon user study could possibly be performed to determine if these differences in image resolution would substantially alter a neurosurgeon's presurgical plan.
Note that our observations are specifically catered to simulation results obtained with the k-Wave Toolbox software package. A comparison with multiple other simulation packages [34][35][36] was beyond the scope of this paper. Using the k-Wave simulation package, we demonstrated for the first time that it is likely sufficient to sacrifice the accuracy of elastic wave simulations by solely relying on compressional wave simulations for presurgical planning. Although this sacrifice affects the photoacoustic waveforms sensed by the simulated receivers (as shown in Fig. 4) and image resolution (as shown in Fig. 6), it may not impact the image quality needed to make surgical decisions based on target detectability and target-to-instrument distance. Thus, we now have a better understanding of this pathway to advance intraoperative transcranial photoacoustic imaging of the ICAs toward surgical use.
Regarding the benefit of shorter simulation times, though minimally invasive neurosurgeries are rarely emergency procedures, surgeons waiting on simulations to complete before operating on a patient might create bottlenecks in presurgical planning. This issue could be a barrier to entry for this technology in hospitals that might not have access to or financial resources for powerful GPUs. Without the GPUs that we used to accelerate our simulations, elastic wave simulations could take several hours to complete and ultimately disrupt the presurgical workflow. When utilizing compressional wave simulations, which require ∼4.6× less time to simulate, surgeons can avoid these presurgical bottlenecks, regardless of their accessibility to GPUs. In addition, reduced computational memory places less strain on computational resource requirements for presurgical simulations and can enable surgeons to perform multiple patient-specific simulations at once. Thus, when implemented with compressional wave simulations, our proposed presurgical planning step has the potential to rapidly identify patient-specific optimal transducer locations for incorporation into the surgical plan, in support of our vision outlined in Fig. 1.
In summary, presurgical identification of optimal photoacoustic imaging system component locations with compressional wave simulations offers three primary benefits. First, this presurgical identification would eliminate wasting valuable operating room time to search for and find a suitable transducer location, thereby reducing the time that the patient is under anesthesia, the total procedure duration, and the medical cost. The second benefit is removal of the potential barrier that the surgeon may be unable to identify a viable transducer location and therefore unable to use photoacoustic image-guidance during the procedure. The third benefit is identification of transducer locations that would minimize image quality degradation from the presence of bone and thereby produce photoacoustic images of the ICAs with the best image quality possible for the patient. Possible future work includes a comparison of simulation types when modeling more complex source distributions, such as the possibility for optical absorption from the deoxygenated blood in the endothelial tubes of the cavernous sinus, and followup neurosurgeon user studies.

Conclusions
The work presented in this paper is the first to reveal that the less time-consuming and less memory-intensive compressional k-Wave simulations are likely sufficient for identification of optimal transducer locations for transcranial photoacoustic-guided surgery. This assessment is based on target visibility and detectability (e.g., target contrast, SNR, CNR, and gCNR) and relative source-to-instrument distance producing similar or identical quantitative measurements for compressional and elastic simulations. These results have multiple implications for reducing barriers for clinical translation of simulation-based photoacoustic-guided surgery.

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