Acoustic impact of the human skull on transcranial photoacoustic imaging

: With balanced spatial resolution, imaging depth, and functional sensitivity, photoacoustic tomography (PAT) hold great promise for human brain imaging. However, the strong acoustic attenuation and aberration of the human skull ( ∼ 8 mm thick) are longstanding technical challenges for PAT of the human brain. In this work, we numerically investigated the impacts of the stratified human skull on photoacoustic wave propagation ( i.e. , the forward model) and PAT image formation ( i.e. , the inverse model). We simulated two representative transcranial PAT implementations: photoacoustic computed tomography (PACT) and photoacoustic macroscopy (PAMac). In the forward model, we simulated the detailed photoacoustic wave propagation from a point or line source through a digital human skull. The wave attenuation, refraction, mode conversation, and reverberation were thoroughly investigated. In the inverse model, we reconstructed the transcranial PACT and PAMac images of a point or line target enclosed by the human skull. Our results demonstrate that transcranial PAMac suffers mainly from wave reverberation within the skull, leading to prolonged signal duration and reduced axial resolution. Transcranial PACT is more susceptible to the skull’s acoustic distortion, mode conversion, and reverberation, which collectively lead to strong image artifacts and deteriorated spatial resolutions. We also found that PACT with a ring-shaped transducer array shows more tolerance of the skull’s adverse impacts and can provide more accurate image reconstruction. Our results suggest that incorporating the skull’s geometry and acoustic properties can improve transcranial PAT image reconstruction. We expect that our results have provided a more comprehensive understanding of the acoustic impact of the human skull on transcranial PAT.

strong acoustic attenuation, and large acoustic impedance mismatch, the skull is a dispersive barrier for acoustic wave propagation [12]. Interestingly, unlike traditional ultrasound imaging, photoacoustic waves are inherently broadband, further complicating the wave distortion through the skull. With only one-way acoustic travel, PAT has great promise for human brain imaging and should be less sensitive to the skull's acoustic distortion than ultrasonic imaging [13]. Nevertheless, the strong acoustic attenuation of the human skull remains a technical challenge for transcranial PAT, which substantially reduces the signal to noise ratio [14]. While the signal strength can be improved by optimizing the light delivery, acoustic detection, and contrast enhancement [15], the image quality of transcranial PAT is still severely deteriorated by the skull's aberration, reverberation, and mode conversion of the acoustic waves [16][17][18][19]. Chemical debridement [10,20], surgical thinning [21], and cranial window techniques [22] have been explored to reduce or mitigate the impact of the skull; however, these methods are invasive and not preferred for functional brain imaging.
Researchers have devoted many efforts toward better understanding the mechanism of photoacoustic wave propagation through the skull. To understand the relationship between the human skull bone density and acoustic parameters, Pichardo et al. measured the longitudinal acoustic wave velocity and attenuation coefficient of freshly excised human skulls [23]. Mohammadi et al. confirmed that the thickness of the skull is the main factor affecting PA signal attenuation and phase distortion [24]. Later, Lee et al. confirmed that the sound pressure was attenuated by 89% through the skull and most of the high frequency energy was lost [25]. The above studies have collectively demonstrated that the human skull is a significant challenge for transcranial photoacoustic imaging. However, none were able to provide a detailed understanding of acoustic wave propagation in a stratified human skull in the context of acoustic attenuation, refraction, mode conversion, and reverberation.
In this paper, we present the first numerical study on transient acoustic waveform propagation through the stratified bone model (Section 2) and human skull (Section 3), especially the mode conversion between longitudinal waves and shear waves within the skull. We also demonstrate the different impacts of the human skull on imaging formation in different PAT implementations (Section 4). For this work, due to the low-frequency ultrasound waves simulated for the human skull, we consider only PACT and PAMac with 1 MHz acoustic waves.

Photoacoustic wave propagation on a digital human skull model
To accurately understand the propagation of photoacoustic waves in the human skull, we employed Wavenology EL (Wave Computation Technologies, Inc., Durham, North Carolina), a commercial simulation platform. Wavenology EL is based on the finite difference time domain (FDTD) method and capable of simulating the acoustic characteristics in complex structures. Using Wavenology EL, we calculated the following first-order partial differential equations in a conservative form and obtained acoustic pressure and velocity fields: where q is the quasi-velocity-strain, f , g and h represents the flux variables in three dimensions, and b is the space-time dependent source vector. Specifically, where ρ is the mass density of the material, v x , v y and v z are the three-dimensional cartesian velocity vectors, ε is the strain tensor, τ is the stress tensor, and b x , b y and b z are the body force vectors. Details about the FDTD method can be found in our previous publications [26][27][28]. Figure 1(a) shows the 3D digital model of a human skull (approximate size: 200 mm×210 mm ×200 mm) used in our simulation. The human skull consists of three layers with different acoustic properties: the inner ivory, the marrow, and the outer ivory. To clearly describe the wave propagation, we put the photoacoustic signal source inside the skull cavity. We selected 1 MHz as the central ultrasound frequency as a balance between the desired spatial resolution and the acoustic attenuation by the skull, based on the previous study [29]. Our FDTD algorithm requires 20 sampling points per wavelength (PPW) in the media to ensure the simulation accuracy. Therefore, the 3D spatial sampling points were 2600×2800×2600 voxels at 1 MHz. To reduce computation load and memory expense, we assumed the skull was relatively homogeneous within the cross-sectional imaging plane of PAT and performed a 2D forward simulation instead of 3D simulations ( Fig. 1(b)). We first tested our simulation platform on a simplified five-layered plane model, as shown in Fig. 2(A). The first and fifth layers are semi-infinite water, with a mass density of 1000 kg/m 3 and a longitudinal wave speed of 1500 m/s. The second layer is outer ivory, with a thickness of 3 mm, a mass density of 1900kg/m 3 , a longitudinal wave speed of 2900 m/s, and a shear wave speed of 1450 m/s. The third layer is bone marrow, with a thickness of 4 mm, a mass density of 1700kg/m 3 , a longitudinal wave speed of 2500 m/s, and a shear wave speed of 1250 m/s. The fourth layer is inner ivory, which has the same parameters as the outer ivory. The detailed acoustic parameters of the skull and water are listed in Table 1 [30]. The conversion of the acoustic attenuation factor is described in the literature [31,32]. The 1 MHz monopole photoacoustic point source (a Ricker profile) is in the first water layer, 80 mm from the boundary.
Because only longitudinal waves are supported in water, wave conversion happens at the interfaces of delamination. To be more specific, at the liquid-solid boundary, the refracted longitudinal waves are partially converted into shear waves; the refracted longitudinal (shear) waves are again partially converted into shear (longitudinal) waves at each solid-solid boundary due to reflection and refraction; at the solid-liquid boundary, all refracted waves are converted to longitudinal waves propagating into water. To analyze the effect of skull attenuation on the  photoacoustic signal, we simulated the acoustic wave propagation in three different scenarios: 1) water without the skull, denoted as 'water-only scenario'; 2) the layered model described in Fig. 2(A), but with lossless bone layers (i.e., zero attenuation factor), denoted as 'bone-loss-free scenario'; and 3) the same layered model as the second scenario, except that the bone layers are lossy, denoted as 'bone-loss scenario'. In the simulation, the origin is set at the center of the x-axis. An omnidirectional point transducer is placed in the top water layer. The photoacoustic point source adopts the Ricker wave function with a center frequency of 1 MHz and a -6 dB bandwidth of 100%, as shown in Fig. 2

(B).
To understand the wave distortions in this simple layered model, we compared the waveforms received by three representative detector points, as shown in Fig. 3. The emission source is located at the origin of coordinates. The coordinates of the three detector points are (50, −50) mm, (50, −26) mm, (50, 0) mm. In the water-only scenario, the acoustic waveform remains unchanged and serves as the reference signal (the solid black lines in Figs. 3(B, E, H)). For both the bone-loss-free and bone-loss scenarios, as the detector gets closer to the direct top of the source, the amplitude of the earliest wave through direct longitudinal transmission becomes stronger, while the later-arrived waves through refraction, reflection, and mode conversion become weaker. This is mainly because both refraction and mode conversion are less prominent as the incident angle decreases. The direct-transmission signal arrives earlier in both bone scenarios than in the water-only scenario, because the longitudinal wave travels faster in the bone. Due to the influence of the loss parameters on the wave propagation speed, the first signal arrives earlier in the lossy model than the loss-free model [34][35][36]. We also observed that the mode-converted signal (the double arrows in Fig. 3(E)) arrives at a similar time in the bone-loss-free scenario as in the water-only scenario, because the shear wave in the bone travels at a similar speed as the longitudinal wave in water. In the bone-loss-free scenario, acoustic waves generated by refraction and mode conversion can oscillate in the solid layers for a long time and thus lead to prolonged signals received by the detector (the dashed blue lines in Figs. 3(B, E, H)). The frequency spectra of the received signals also become more fluctuated due to constructive and destructive interferences that are closely relevant to the layers' dimensions (the dashed blue lines in Figs. 3(C, F, I)). In the bone-loss scenario, the amplitude of the received signals, especially at relatively high frequencies, is reduced by at least 75% (the dot red lines in Figs. 3(C, F, I)). The strong acoustic attenuation by the bone layers also quickly diminishes the reverberations inside the bone and thus shortens the duration of the received signal. To better understand the skull's distortion on different target shapes, we compared the waveforms in the bone-loss scenario using a point source and a line source, detected by three representative individual detectors (Figs. 4(A, C, E)). The line source is located along the x-axis from -5 mm to 5 mm. The emitted acoustic pressure of the line source is the same as that of the point source. Although photoacoustic signals from a point source are spherical waves and thus diverge quickly, photoacoustic signals from a line source are cylindrical or elliptical waves and thus diverge slowly. When the incident angle is small, the signals received by the detector directly above the source are highly consistent between the two sources ( Fig. 4(F)). As the incident angle increases, the received signal amplitude from the line source drops much more rapidly than that from the point source (Figs. 4(B-D)). The results demonstrate that the impact of the skull is highly dependent on the target shape.  Figure 5 and Visualization 1, Visualization 2, and Visualization 3 further demonstrate the complete photoacoustic wave propagation in the planar bone model, which are also briefly described here. When the point source in water is excited, the acoustic wave spreads around as a spherical wave until it encounters a bone boundary. Due to the mismatch in acoustic impedance between the bone and water, there is reflection, refraction, and mode conversion at the water-bone interface. The refracted longitudinal waves are partially converted to slower shear waves into the bone. The shear waves, however, have a smaller refraction angle and thus less wavefront distortion than the longitudinal waves. The longitudinal waves and shear waves continue propagating inside the bone, and again there is reflection, refraction, and mode conversion at the bone-water interface. The refracted shear waves are completely converted into longitudinal waves because water does not support shear waves. Without attenuation in the bone, the waves can oscillate in the bone and transmit into the water for a prolonged period (Fig. 5(A), Bone-loss-free and point source scenario). With attenuation inside the bone, the residual waves are quickly diminished (Fig. 5(B), Bone-loss and point source scenario). With a line source, the acoustic wave resembles a cylindrical wave and is better detected just above the source (Fig. 5(C), Bone-loss and line source scenario). After validating our simulation platform on a planar bone model, we simulated the photoacoustic wave propagation through a stratified human skull model ( Fig. 1(B)). In the simulation, the human skull is immersed in water. The soft brain tissue in the skull is replaced by water even though the acoustic properties of brain tissue are slightly different [17,30,37]. The coordinate origin is set at the center of the model. The length of the model along the x-axis is 205 mm and the width along the y-axis is 200 mm. We set the point source at the origin of coordinates. Three omnidirectional detectors, located at (−95, 90) mm, (−95, 50) mm, and (−95, 0) mm, are used to receive acoustic waves in space. The acoustic parameters of the human skull are summarized in Table 1.
To quantitatively analyze the impact of the human skull on the acoustic wave propagation, we again compared three scenarios: water-only, bone-loss-free, and bone-loss, as shown in Fig. 6. We used the waveforms received by the detectors in the water-only scenario as the ground truth, as shown by the solid black lines in Figs. 6(B, E, H). In the bone-loss-free and bone-loss scenarios, the early arrived signals via direct longitudinal transmission have greatly varied amplitudes, heavily dependent upon the relative position of the detector and the source as well as the curvature of the bone. The skull's irregular geometry results in complex local boundary conditions that largely determine the refraction and mode conversion activities. For example, detector point 3 is directly above the source, and the early arrived waves in the bone-loss-free and bone-loss scenarios are 41% and 78% weaker, respectively, than in the water-only scenario ( Fig. 6(H)). When the detectors are farther away from the source, such as at detection point 2, the signal amplitudes are reduced by 58% and 95%, respectively (Fig. 6(E)). Interestingly, when the detector is even farther away, such as at detection point 1, the signal drop becomes less significant, most likely because the local bone curvature better matches the incoming acoustic wavefront. Such complex dependence on the local bone geometry highlights the necessity of accurately modeling the acoustic wave propagation in transcranial PAT. The bone-loss scenario exhibits quickly diminished reverberations inside the bone, which shortens the received signal duration. Figures 6(C, F, I) show the frequency spectra of the received photoacoustic signals. The high frequency components are attenuated more in the lossy model. Similar to the planar bone model, the signal frequency spectrum through the human skull has drastic fluctuations due to the constructive and destructive interferences at the boundaries. Such interferences result in the strengthening or weakening of certain wavelengths, depending on the local bone dimensions and acoustic properties. We also simulated the acoustic wave propagation from a line source and compared the results with the point source. The line source is located on the x-axis from -5 mm to 5 mm. The signals received by the detector directly above the line source are almost the same as the point source, as shown in Figs. 7(H, I). This observation is consistent with the planar bone model (Fig. 4(F)).
However, the signal from the line source drops quickly as the detector moves away from the source, as shown in Figs. 7(C-F), mainly due to the cylindrical signal profile of the line source. With wave refraction, reverberation, and mode conversion in the skull, the received signal profile is rather complex. For example, at detection point 2, we observed that there are at least three major groups of signals that arrive relatively early (Fig. 7(E), blue box). The first group of signals, the longitudinal waves that transmit directly from the source to the detector through the skull, arrives the earliest due to the high speed of longitudinal waves inside the bone. The second group of signals is the direct transmission of the longitudinal waves that are first converted into shear waves and then back to longitudinal waves. The third group of signals are mainly reverberation waves. We observed the same phenomenon in our previous work in mouse skulls at much higher frequencies [12]. Figure 8 and Visualization 4, Visualization 5, and Visualization 6 show the complete acoustic wave propagation from a point source and a line source in the human skull model with and without loss.

The human skull's effects on PAT image formation
The impact of the human skull on transcranial PAMac and PACT images were further investigated. For PAMac, a single-element focused ultrasonic transducer with a central frequency of 1 MHz and a focal length of 90 mm is used to detect the time-resolved PA signals. 1D scanning is performed along the x-axis to acquire PA signals at each location. The PAMac image is then formed by directly converting the received time-resolved signals back into the spatial domain, assuming a uniform speed of sound of 1500 m/s (Figs. 9(A-B)). The skull induces error in the reconstructed source position along the y-axis, i.e., the acoustic axis of the focused transducer (Figs. 9(C-D)). The fast longitudinal wave inside the skull is the main reason for the position error.
For PACT, a linear array of 191 detectors is evenly distributed above the skull surface, between (-95, 90) mm and (95, 90) mm. The traditional delay-and-sum method and time-reversal method are used for image reconstruction. The reconstructed PACT images of a point target and a line target demonstrate that the skull not only results in a large error in the target position, but also induces severe image artifacts (Fig. 10). Because the signals received by all detectors are used to reconstruct the final image, PACT is more vulnerable than PAMac to wave distortions due to refraction, mode conversion, and reverberation. More importantly, the wavefront distortions perceived by each of the detectors are highly dependent on the local skull geometry. Without knowing the exact distortion pattern and assuming the homogenous speed of sound, the image reconstruction by the traditional delay-and-sum method essentially amplifies the errors contributed by each detector.
We also used the time-reversal algorithm to reconstruct the PACT image, which is more adaptive than the delay-and-sum method. When the skull's structure and acoustic properties are accounted for in the time-reversal algorithm, the reconstructed image clearly has fewer artifacts than the delay-and-sum result, as shown in Figs. 11(C-D). The reconstructed position of a line target also shows better accuracy than the delay-and-sum result, as shown in Figs. 11(E-F). These results demonstrate that combining the time-reversal reconstruction method with accurate acoustic properties of the skull can substantially improve the transcranial PACT image quality.
To further investigate the skull's impact on different detection modes, we simulated a ringshaped ultrasonic transducer array with a total of 360 elements uniformly distributed around the skull (radius: 95 mm). Here, we also used the time-reversal method to reconstruct the PACT images. When the skull's acoustic properties are not accounted for in the time-reversal reconstruction process (Figs. 12(A-B)), the image has strong artifacts, like with the linear array result. However, using the same reconstruction method, the reconstructed target position with the ring-shaped array is more accurate than that with the linear array (Figs. 10(C-D-E-F)).
We believe that the more accurate target location is mainly due to the large acceptance angle of the ring-shaped array, which helps to cancel out the impact of the irregular geometry of the skull. Furthermore, when the skull's structure and acoustic properties are accounted for during the reconstruction process, there are far fewer artifacts in the ring-shaped array result, as shown in Figs. 12(C-D-E-F) than in the linear array result (Figs. 11(C-D-E-F). These results demonstrate that the ring-shaped array is advantageous over the linear array for transcranial PACT.

Conclusion and discussion
In this work, we simulated photoacoustic wave propagation through a stratified human skull model and numerically investigated the skull's impact on transcranial PAT. The numerical simulation results reveal the wave reflection, refraction, mode conversion, and reverberation process at the skull-tissue interface and within the skull. The reconstructed PACT and PAMac images confirm that phase aberration, mode conversion, and reverberation are the main causes of image deterioration. We also demonstrate that, even with similar ultrasound frequencies and effective detection apertures, transcranial PAMac results in better image quality than transcranial PACT. Transcranial PACT using a linear transducer array is heavily affected by the skull's acoustic distortion and the limited detection angles, which can cause severe reconstruction artifacts. However, it must be noted that PACT still has advantages over PAMac in penetration depth and imaging speed for functional applications.
While our simulation work provides insights on the engineering challenges associated with transcranial PAT, our results have led to several potential solutions to improve the transcranial PAT image quality. First, we confirmed that the main reason for artifacts in the transcranial PACT image is the phase aberration caused by the inconsistent wave speed in the skull and soft tissues [12]. Accounting for the structural and acoustic parameters of the skull, the time-reversal reconstruction method can obtain relatively accurate results, as shown in Figs. 10(E, F). Therefore, one solution is to use the full-wave simulation method to model the complex and inhomogeneous structures and incorporate the wave propagation in the image reconstruction. Several relevant methods have been proposed to reduce the effect of the human skull on transcranial photoacoustic imaging. Because shear waves in the skull have less impedance mismatch than longitudinal waves, G. T. Clement et al. proposed a shear wave induction technique to reduce the wavefront distortion [33]. Javaherian et al. presented an optimization inverse framework based on coupled equations [38]. Mohammadi et al. proposed a vector space similarity model to improve the image quality [39]. The Wang group investigated the shear wave propagation in the solid planar geometry [37] and proposed a de-aberration layered universal back-projection (L-UBP) method [29]. Nevertheless, all these methods can benefit from accurate skull structure and acoustic parameters, which may be often difficult to obtain in practical medical applications.
Our results have shown that wave reverberation inside the skull induces artifacts in both transcranial PACT and PAMacs. We demonstrate that the detector just above the target receives the largest signal amplitude with the least phase aberration and mode conversion ( Fig. 6(H) and Fig. 7(H)). Therefore, it may be beneficial to utilize only a subset of transducer detectors in PACT according to the target location, which will lead to a trade-off between spatial resolution and image artifacts. Generally speaking, PAMac images are less affected by the wave distortion by the skull because of the strong skull loss and focused detection mode, but PAMac is usually too slow for functional brain imaging [40]. Therefore, improving the imaging speed of PAMac by advanced scanning mechanisms may provide a practical solution for transcranial brain imaging. For transcranial PACT, it is also important to increase the total detection angle coverage. In our models, the ring-shaped ultrasonic detector array improved the signal strength and reduced the impact of the skull's acoustic distortion [41][42][43]. Half-time image reconstruction, in which the reverberation signals are separated from the direct transmission signal, can also be explored to reduce the effect of reverberation [44].
In summary, we expect this work can lend understanding to the impact of the human skull on transcranial PAT and inform the development of innovative methods to mitigate such impacts.