Numerical evaluation of the skull for human neuromodulation with transcranial focused ultrasound

Objective. Transcranial focused ultrasound is an emerging field for human non-invasive neuromodulation, but its dosing in humans is difficult to know due to the skull. The objective of the present study was to establish modeling methods based on medical images to assess skull differences between individuals on the wave propagation of ultrasound. Approach. Computational models of transcranial focused ultrasound were constructed using CT and MR scans to solve for intracranial pressure. We explored the effect of including the skull base in models, different transducer placements on the head, and differences between 250 kHz or 500 kHz acoustic frequency for both female and male models. We further tested these features using linear, nonlinear, and elastic simulations. To better understand inter-subject skull thickness and composition effects we evaluated the intracranial pressure maps between twelve individuals at two different skull sites. Main results. Nonlinear acoustic simulations resulted in virtually identical intracranial pressure maps with linear acoustic simulations. Elastic simulations showed a difference in max pressures and full width half maximum volumes of 15% at most. Ultrasound at an acoustic frequency of 250 kHz resulted in the creation of more prominent intracranial standing waves compared to 500 kHz. Finally, across twelve model human skulls, a significant linear relationship to characterize intracranial pressure maps was not found. Significance. Despite its appeal, an inherent problem with the use of a noninvasive transcranial ultrasound method is the difficulty of knowing intracranial effects because of the skull. Here we develop detailed computational models derived from medical images of individuals to simulate the propagation of neuromodulatory ultrasound across the skull and solve for intracranial pressure maps. These methods allow for a much better understanding of the intracranial effects of ultrasound for an individual in order to ensure proper targeting and more tightly control dosing.


Introduction
Transcranial focused ultrasound (tFUS) is an appealing approach for noninvasive neuromodulation of cortical tissue for a wide variety of applications, including those with deeper cortical targets as well as deeper structures in the brain beyond the cortex. Other non-invasive technologies, such as transcranial magnetic stimulation (TMS), which uses magnetic fields that pass unimpeded through the skull, have proven an effective form of noninvasive neuromodulation [1]. However, despite the advantage of passing through the skull unimpeded, TMS produces electric fields in the brain with a spatial scale of centimeters [2] and with peak effects limited to the cortical surface [3], limiting the spatial and depth specificity of neural stimulation. In contrast, tFUS can offer a higher spatial resolution on the scale of millimeters, and has been used safely and effectively for cortical neuromodulation in mouse [4], rabbit [5], monkey [6], and humans [7][8][9][10]. However, unlike TMS, ultrasound does not pass unimpeded through the skull. The skull represents the primary barrier to ultrasound due to its high frequency dependent attenuation, dispersion, and refraction of ultrasound waves which results in a significant loss of energy and distortion of the transmitted ultrasound beam.
Our previous work using computational modelling of neuromodulatory transcranial focused ultrasound using a simplified skull found the inclusion of brain tissues and geometry to have little impact on transcranial pressure maps and to result in safe levels of heating in the skull and brain tissue [11]. However, changing the simplified skull geometry from flat to curved reduced peak intensity by 40%, affirming the skull as the major determinant of intracranial pressure and thus recommending the use of detailed skull models to more accurately determine intracranial pressure. It is the focus of this work to develop a model of transcranial focused ultrasound using detailed representations of the skull while regarding all other tissues as water, as the modeling of soft tissues has already been shown to have negligible effects [11]. Detailed computational models of the skull have been used previously in applications including aberration correction in high intensity focused ultrasound for thermal ablation [12] and non-thermal ablation with microbubble enhanced focused ultrasound [13]. In the present study, we developed computational models of the skull based on medical images to evaluate simulation methods for tFUS for noninvasive neuromodulation and to explore the differences in intracranial pressure maps between individuals. Evaluations of linear and nonlinear simulations, the effects of the skull base and size, elastic simulations, and low frequency ultrasound simulations provide further insight on the transcranial wave propagation of ultrasound for human neuromodulation. These issues are explored in general using one female and one male skull model in order to provide a more balanced discussion of modeling methods, while a larger dataset is later used to discuss the variability of intracranial pressure maps between individuals.

Methods
Computational models were developed using images from magnetic resonance (MR) imaging and computerized tomography (CT) of individuals to evaluate the wave propagation of tFUS across the skull and the resultant intracranial pressure maps. Simulations were performed using the k-Wave MATLAB toolbox [14], which uses a pseudospectral time domain method to solve discretized wave equations on a spatial grid. k-Wave is one of a few readily available numerical solvers that has been extensively used in the literature [15,16]. Additionally, k-Wave has undergone thorough numerical testing for highly heterogeneous transcranial simulations and has been shown to give equal or better accuracy compared to finite-difference time domain methods [17,18]. Acoustic simulations using k-Wave were used on baseline models to evaluate simulation methods before investigating the variation in intracranial pressure maps due to individual differences.

Medical images
The baseline acoustic models utilized in this study to establish modeling methods were constructed using female and male datasets from the Visible Human Project ® (VHP) [19]. Each dataset consists of MR, CT and cryosection images taken in the axial plane of the head at various slice thicknesses (4 mm MR, 1 mm CT). CT images were used to construct the acoustic model of the skull, while MR images were used to target tFUS at the primary motor hand representation (M1) based on brain morphology. The precentral gyrus was identified for targeting according to morphological landmarks outlined in previous literature, where the 'hand knob' was used on transverse sections to identify M1, and verified on sagittal sections using the 'hook sign' [20,21]. In addition to M1 placement, we also tested the effect of transducer placement at the temporal window. Placement of the transducer in simulations of ultrasound through the temporal window was determined by identifying the pterion based on the skull morphology observed in CT scans, approximately 2.6 cm behind and 1.3 cm above the posterolateral margin of the frontozygomatic suture [22]. The pterion is where the frontal, parietal, temporal, and sphenoid bones of the skull unite, and is regarded as the thinnest and weakest part of the skull. Following identification of the target, the ultrasound source was placed coaxially with the target, 1.26 mm from the outer surface of the skull. The skull was then rotated to orient the outer skull surface to be approximately tangential to the computational source as determined by visual inspection of orthogonal planes intersecting the acoustic axis (centerline) of the source. For comparisons between models, skull thicknesses below the transducer were determined along the acoustic axis of the source.
To evaluate differences in intracranial pressures between individuals, an additional dataset of CT and MR images (N = 10) that were originally acquired for surgical planning and verification of deep brain stimulator lead implantation at the Center for Magnetic Resonance Research (CMRR) at the University of Minnesota was used. CT and MR image sets with 1.5 mm and 1 mm slice thickness respectively were used for modeling tFUS similarly as described for the previous baseline medical images. The only additional prep aration needed for these image sets was to remove the deep brain stimulator leads within the skull from CT scans, which was done manually using ImageJ [23]. Each image slice of the medical scans was inspected visually, and if a DBS electrode was present, an encapsulating circle was drawn over the DBS lead and imaging artifact, and the enclosed region specified as water. This process was repeated up to the region of the burr hole in the skull. Once at the burr hole, there was no further manipulation of the image to either fill in the burr hole in the skull or remove further wires and leads on the outside of the skull. If the placement of the ultrasound transducer on the outside of the skull to target intracranial sites of interest were to require placement of the ultrasound transducer in the vicinity of either the burr hole or tunneled extension wires of the DBS leads, that image set would be excluded from use for simulations of transcranial ultrasound. Altogether, images of 12 individuals (3 females and 9 males) with an average age and standard deviation of 62.2 ± 14.8 years at the time of medical imaging were utilized to construct acoustic models. The variation of intracranial pressure maps across individuals was evaluated using linear regression models and statistically tested with a significance level of 0.05.

Computer simulations
CT and MR images were first co-registered in MATLAB [24] to guide targeting of ultrasound to the 'hand knob' representation of the precentral gyrus based on cortical geometry. The images were then resampled using linear interpolation for acoustic simulations at a finer resolution of 0.42 mm and the acoustic parameters for simulation calculated from the CT images. The skull was extracted manually using a threshold intensity value of 150 HU and the intracranial space was assumed to be homogenous as ultrasound reflections between soft tissues are small [11]. By considering only the skull for our computational models, we reduced model complexity and computation time while including the main determinant (skull bone) of transcranial ultrasound wave propagation behavior [11,25]. Acoustic parameters were calculated from CT data assuming a linear relationship between skull porosity and the acoustic parameters [12,26]. For elastic simulations requiring shear acoustic parameters, they were set proportional to the compressional parameters as relationships for the shear parameters as a function were not found. Simulations were then carried out using software and services provided by the Minnesota Supercomputing Institute at the University of Minnesota. While these modeling methods may introduce uncertainties in results, such methods have been used successfully for aberration correction and forward problem simulations in previous literature [12,13,26,27].
The computational model of the ultrasound transducer used in simulations was constructed to recreate empirical acoustic pressure maps of focused ultrasound transmitted in an acoustic test tank, similar to previous work [11]. The ultrasound transducer modeled was a custom-designed singleelement focused transducer (Blatek, Inc., State College, PA) with a center frequency of 0.5 MHz, a diameter of 30 mm, and a focal length of 30 mm. The transcranial ultrasonic neuromodulation waveform that the simulations were based on has been previously described [7], and has an acoustic frequency of 500 kHz with a pulse duration of 360 µs (180 cycles) pulsed at a repetition frequency of 1.0 kHz over half a second. Simulation of the transducer model in free water produced a max pressure of 1030 kPa.
Skull models were considered to be immersed in water in order to couple the model to the ultrasound source. The acoustic properties of the skull model were calculated from CT Hounsfield units (H) based on a porosity (ψ) calculated using equation (1) [12,26]. Following calculation of porosity, the compressional speed of sound (c skull,c ), density (ρ skull ), and attenuation (α skull,c ) of the skull were calculated using equations (2)-(4) [12,26] with the parameters given in table 1. The minimum and maximum compressional attenuation of skull was taken from Connor 2003 [28] and corresponds to the extreme values of attenuation listed at 500 kHz. The attenuation of water was calculated for a frequency of 500 kHz at 37 °C, corresponding to body temperature, using the built in function in k-Wave based on previous literature [29]. For elastic wave simulations, shear wave speed (c skull,s ) was set as (1400/2550)c skull,c and shear wave attenuation (α skull,s ) was set as (90/85)α skull,c [30]. Though acoustic wave simulations do not include the effects of shear mode conversion and attenuation that elastic wave simulations do, they impose less of a computational burden than elastic wave simulations due to the reduced number of field quantities and parameters (stress tensors and shear mode parameters). However, significant shear modes can be generated when the incident wave angle is greater than a critical angle [31], thus we evaluated both simulation methods to determine the amount of influence shear waves have on intracranial pressure maps.
While the accuracy of the pseudospectral time domain method decreases in inhomogeneous media, the error reduces with increasingly finer spatial discretization [32]. A spatial discretization of 7 points per wavelength in water (0.42 mm) was chosen for our simulations, which results in simulation  [29] c bone,c = 3100 [12] ρ bone = 2200 [12] α min,skull,c = 21.5, α max,skull,c = 208.9 [28] c bulk = 2850 [31] ρ bulk = 1732 [31] α bulk = 85.0 [31] domains for the baseline VHP female and male skulls of approximately 500 × 500 × 500 grid points. The full skull models of this size take about 30 h to solve on the high performance computing cluster. Previous convergence testing in 3D transcranial ultrasound models found that a discretization of 6 points per wavelength was sufficient to reduce the max pressure error to less than 5%, while the full width at half maximum (FWHM) volume achieves convergence at approximately 6 points per wavelength [18]. The maximum and minimum values of pressure were recorded for each point in the simulation domain, and were used for comparison between models and calculations of intracranial FWHM volumes. A Courant-Friedrichs-Lewy stability criterion (CFL number) of 0.1 was used for all simulations, except for elastic wave simulations which required a smaller CFL number of 0.05 for a stable solution. The default simulation length in k-Wave is determined by the amount of time it takes the ultrasound frequency of interest to travel to geometrically opposite corners of the simulation grid at the minimum speed of sound contained within the simulation domain. This resulted in a default simulation including ~115 cycles of 500 kHz ultrasound in our baseline models. Our typical transcranial ultrasonic neuromodulation waveform has 180 cycles, so simulations including 180 cycles were run for comparison. For the baseline VHP female and male skulls, less than 0.2% change was found in the peak intracranial pressure, along with less than 2% change in FWHM volume between the two simulation time lengths. These changes were regarded as not significant enough to warrant the increased simulation time associated with forcing 180 cycles, thus the default simulation lengths were used for model analyses.
To evaluate the generation of standing waves by ultrasound, standing waves were estimated using a spatial filter designed to quantify the pattern of modulation of the pressure map by constant interference as outlined in previous literature [33]. Overall, the fast oscillations of the pressure distribution are isolated along each dimension with a high pass 11th-order Butterworth filter with a cut-off of half the wavelength of ultrasound in water at the acoustic frequency of interest. A Hilbert transform is then applied to obtain the envelope of the isolated high frequency modulations, and the maximum of the x, y, or z component is taken to produce a single map estimating the spatial modulation pattern. This pattern is largely associated with the presence of standing waves, as some artifacts are produced by the high spatial frequencies of oscillations in the pressure map of a focused wave at the location of the focus [33].

Baseline linear acoustic simulations
Simulations were first run on the VHP female and male skulls targeting the hand knob of M1 to establish modeling methods (figures 1 and 2). The intracranial max pressure produced in the female and male skulls respectively were 177.32 kPa and 137.79 kPa. The presence of the skull in the model reduced the peak pressure produced by the transducer in free water by 6 to 7 fold. The skull also distorted the beam from the axis of the transducer by 3.17° and 1.08° in the female and male skulls respectively, as measured by the major centroidal principal axis of the FWHM volume. As the skulls were manually rotated in an effort to orient the transducer tangentially to the skull surface over the M1 location, the differences in the dist ortions of the beam from the axis of the transducer are dependent on the degree of tangential misalignment of the transducer and skull surfaces to some degree, with the individual geometry of the skulls also contributing to the dist ortions of the ultrasound beam. Since the orientations of the female and male skulls with respect to their transducers were constant between simulation models, comparisons between cases for each of the respective skull models are still valid and informative.

Nonlinear acoustic simulations
To investigate the contribution of nonlinear effects to the intracranial pressure maps, simulations were run with a nonlinear acoustic model and compared to the baseline linear acoustic models. The B/A nonlinearity parameter was set to 5 for water and 374 for bone [34]. While the value of the B/A nonlinearity parameter has not yet been satisfactorily determined in the skull, as the Goldberg number is expected to be small due to the high attenuation of the skull [35], we used a relatively high number throughout the skull compared to the rest of the model to get a conservative estimate of the impact of model nonlinearity. Between the linear and nonlinear models, there was less than a 1% change in both the maximal pressures and FWHM volume, with virtually no difference in pressures along the main beam path (maximum pressure difference of 1.3% located in the far field). Additionally, the location of intracranial max pressure was identical between linear and nonlinear models. These results reinforce that at low source intensities and pressures, wave propagation is well described by linear acoustic theory.

Skull base effects
To reduce model complexity and simulation time, it may be attractive to simulate a small subspace that encompasses the transducer source, focal region, and a subsection of the skull along the acoustic beam path. However, simulation of this subspace would neglect the interaction of reflected ultrasound from the skull base and other intracranial surfaces with pressures at the acoustic focus. To explore the effect of the skull base on intracranial pressure maps, linear acoustic simulations were run on the VHP female and male skulls excluding their bottom halves (figures 3 and 4). The peak pressures increased in the baseless models by 2.0% in the female and 2.4% in the male, while the FWHM volume increased by 12.9% in the female and decreased by 0.8% in the male. The location of peak pressure shifted by 2.6 mm and 1.0 mm in the female and male skulls respectively. These results suggest that skull base reflections can slightly influence the intracranial pressure map, dependent on the skull, as the larger volume change and translation of maximum pressure location noted in the baseless female skull is likely due to the smaller size of the female skull facilitating effects from reflections compared to the male skull.
To explore the role of skull size further, we constructed a model of transcranial focused ultrasound using CT scans from a rhesus macaque with 0.6 mm slice thickness. The focal distance of the modeled transducer was halved to 15 mm to similarly keep the main beam profile near the inside surface of the skull by the transducer without extending to the opposite side, as in the human models. This was accomplished by decreasing the radius of curvature of the bowl shaped transducer by 43% while keeping the aperture constant to change solely the focal distance of the transducer. The strength of the modified transducer was then scaled to produce an identical max pressure as the previous, unmodified transducer. Linear acoustic simulations were then run on the macaque skull with and without the bottom half of the skull (figure 5). Compared to the full macaque skull, the baseless macaque skull had a 26.4% higher max pressure, a 1.5% decrease in FWHM volume, and a translation of max pressure location of 3.4 mm. These results suggest that skull base reflections can play a larger role in animal models with smaller skulls compared to humans. Greater wave confinement will occur in smaller skulls, resulting in stronger and more numerous interactions, however, predicting whether wave interactions will be constructive or destructive is difficult to predict for unideal geometries.

Elastic simulations
Due to the increased complexity of accounting for elastic waves in modeling, simulations of elastic wave propagation were conducted in a restricted subsection of the skull directly under the transducer as shown in figure 6. A linear acoustic simulation using an identical setup with the same subsection of the skull as in the elastic model was constructed to facilitate comparisons between the linear acoustic and elastic simulations. For the subsection, elastic simulations took approximately 14 h to solve on the high performance computing cluster, while acoustic simulations took half an hour. Between the two models, the female elastic skull had an increase of 3.5% in max pressure, decrease of 8.9% in FWHM volume, and a translation of max pressure location by 0.6 mm. The male skull had an increase of 14.5% in max pressure, increase of 15.0% in FWHM volume, and a translation of max pressure location by 0.9 mm. While shear waves do not propagate in liquid media and the conversion from compressional to shear waves inside the skull is generally not considered significant at incidence angles lower than 20° [31], these results from a subsection of skull suggest that the inclusion of shear waves in elastic models for more accurate intracranial pressure maps may be warranted, though the simulation of elastic wave models is considerably more resource intense than that of acoustic wave models, especially for full skull models.

Homogenous skull properties
The effect of considering the skull to be a homogenous medium with consistent bulk properties was investigated to gauge how much of an influence the mapping of material properties using medical images influences intracranial pressure maps, all while keeping the realistic skull geometries constant. This reduces the complexity of the model and the amount of time required for simulation. The skull was assigned bulk properties as listed in table 1 and the comparison of the intracranial pressure maps to the baseline models are shown in figures 7 and 8. In the homogenous female skull, the max pressure increased by 18.3%, the FWHM volume decreased by 2.1%, and the point of maximal pressure translated by 0.6 mm compared to the baseline model. In the homogenous male skull, the max pressure increased by 58.7%, the FWHM volume increased by 16.5%, and the point of maximal pressure translated by 5.1 mm. These results seem to indicate that simplifying the skull by assigning it homogenous properties may be sufficient for some skull models with regard to the position of peak pressures and FWHM volume, such as with the female skull model, but does not necessarily hold for all skull geometries, as in the case of the male skull model which had larger changes in the FWHM volume and location of peak pressure.

Transmission through the temporal window
The temporal window is considered an 'acoustic window' in that it is one of the areas of the skull that offers the least resistance to the transmission of ultrasound due to it being the thinnest part of the skull. We modeled the transmission of ultrasound using the same source transducer as in previous female and male VHP skull models through the temporal window to investigate the effect on the intracranial pressure maps ( figure  9). The skull thickness below the transducer and deviation of peak pressure location from their free field pres sure maps is reported in table 2. The intracranial pressures were greater for transmission through the temporal window compared to the ultrasound pressure maps targeting the hand knob by 23.8% in the female model and 84.6% in the male model, while FWHM volumes were smaller by 71.0% in the female model and 30.4% in the male model. These results indicate that the temporal window can provide a transcranial access site for peak pressures and optimal focusing using ultrasound.

Low frequency ultrasound transmission
Among the factors that contribute to the choice of acoustic frequency, the frequency dependent transmission of ultrasound through the skull is one key factor, where lower frequencies generally have a higher transmission coefficient through the skull [36,37]. To test how much a lower acoustic frequency may change the intracranial pressure maps of transcranial focused ultrasound, we ran simulations of the baseline linear acoustic models of the VHP female and male skull at 250 kHz; half of the original 500 kHz (figure 10). The model ultrasound transducer was first modified to ensure a focal length of 30 mm and identical pressures at the lower frequency. We were unable to match the FWHM volume of the original 500 kHz transducer without substantially changing the aperture of the 250 kHz model transducer, so the FWHM volumes are not presented for comparison.
In the female skull, the intracranial max pressure decreased by 5.9% targeting the hand knob and 39.2% through the temporal window compared to 500 kHz simulations. In the male skull, the intracranial max pressure decreased by 12.8% targeting the hand knob and 14.6% through the temporal window compared to 500 kHz simulations. Comparing the locations of max intracranial pressure, in the female skull the location shifted by 10.2 mm targeting the hand knob and 14.7 mm transmitting through the temporal window, while the points similarly translated by 6.00 mm and 4.33 mm in the male skull while targeting the hand knob and transmitting through the temporal window respectively. Additionally, the deviations of the points of peak pressure using 250 kHz and 500 kHz compared to their peak pressure location in free water is given in table 2. These results seem to indicate that while the transmission coefficient of lower frequencies through the skull may generally be higher, the skull geometry and thickness may influence the transmission of different frequencies to a varying degree, producing lower pressures and translating the point of peak pressure at 250 kHz compared to 500 kHz, as similarly observed in previous literature that took source location relative to the skull into account as well [38,39]. Regarding the point of peak pressure, table 2 shows that the focal point will deviate less from free field simulations at 250 kHz when propagating through the temporal window. This can be explained by the longer wavelength of 250 kHz in water (approximately 6 mm) being greater than the thickness of the temporal window, and therefore experiencing less aberration by the skull layer compared to 500 kHz (wavelength of approximately 3 mm in water). Also notable in the 250 kHz pressure maps of figure 10 is the greater prominence of standing waves throughout the intracranial space due to reflections compared to simulations at 500 kHz. To inspect these aspects further, the intracranial pressure maps were spatially filtered to allow features greater than half the wavelength of ultrasound in water at the frequency of interest. A qualitative comparison of figures 11 and 12 readily reveal that standing waves due to reflections within the intracranial space are more prevalent at 250 kHz compared to 500 kHz.

Variation between individuals
Due to differences in skull morphology and composition, the intracranial pressure maps from tFUS can vary between individuals. To investigate this, we simulated 500 kHz linear acoustic ultrasound targeting the hand knob and ultrasound transmitted through the temporal window in 10 additional individuals beyond the VHP female and male models. The intracranial max pressures and FWHM volumes were determined in each of these models and assembled for analysis ( figure 13). Overall, the mean and standard deviation of max pressures was 142.4 ± 37.1 kPa and of FWHM volumes was 1471.4 ± 944.7 mm 3 . A linear relationship was not found between skull thickness and max pressure or FWHM volume (p-values > 0.05 in figures 13(a) and (c)). To account for skull composition, max pressures and FWHM volumes were also analyzed with respect to the integral of density along the thickness of the skull through the central axis of the transducer, and a linear relationship was also not found using this metric (p-values > 0.05, data not shown). Figures 13(b) and (d) show the data paired by individual, where temporal skull thickness is thinner than the skull thickness when targeting the hand knob.

Discussion
Transcranial focused ultrasound for neuromodulation is appealing for its noninvasiveness, high spatial resolution, and ability to focus at shallow or deep focal lengths. However, the  current research landscape for ultrasonic neuromodulation uses diverse stimulation parameters with variable results. One source of variability is the uncertainty in knowing the exact location and dosing of ultrasound when used noninvasively.
To work towards addressing this issue, we have developed and assessed computational modeling methods based on medical images for transcranial focused ultrasound using k-Wave, a readily available MATLAB toolbox that has been extensively utilized and evaluated in previous literature. We assessed a number of acoustic models, and found linear acoustic simulations with acoustic properties linearly mapped based on medical image intensities to be suitable to solve for intracranial pressure maps. Elastic simulations may be desirable for higher accuracy as we found differences compared to linear acoustic simulations of 15% or less regarding the max pressure and FWHM volume, but the complexity of elastic simulations with their tighter sampling criteria and restrictions of simulation space to arrive at a solution preclude its utility. We also found that the skull base has effects on the intracranial pressure map, dependent on the overall size of the skull. Larger skulls can likely suffice without the skull base when solving for intracranial pressures, while smaller skulls will have more reflected interference with the main beam path and require the skull base in order to more accurately determine the pressures within the main beam path. Frequency of ultrasound will also play into consideration of the skull base, as the smaller female skull generally experienced greater changes at 250 kHz compared to 500 kHz, relative to the male skull. Geometry also seemed to play a factor in the exploration of homogenous skull models, as the FWHM volume and location of max pressure in the female skull model did not change overall from the baseline model, but did in the male skull model. Furthermore, we simulated the intracranial pressure maps at different ultrasound targeting locations in a collection of individuals to gain an appreciation of the variability between individuals and their skull geometries and thicknesses. By using computational models to investigate and examine the intracranial pressure maps of individuals, we can help explain and control variability to advance the utility of focused ultrasound methods in research for human neuromodulation. The values and functions used in this work to convert medical images into maps of acoustic properties were based on previous methods reported in literature, but variations in the assumed methods and values can have a large influence on the simulated intracranial pressure maps. Previous work has found that at 500 kHz the peak pressure, peak location, and FWHM volume can readily change when the speed of sound of acoustic property maps are varied, albeit less so when the density and attenuation maps are varied [25]. The inherent variations in maps of acoustic properties between individuals likely played a large role in our not finding any simple linear relations to predict max pressure and FWHM volume in our analysis of intracranial pressure maps between individuals. Indeed, just a change of 1 mm in skull thickness has been reported to change peak pressure by about 20%, peak location by 2 mm, and FWHM volume by about 5% [25]. However, due to the methods of assigning acoustic properties being constant between all models of the individuals used in this work, the generality of trends and observations made holds despite any inaccuracies in the methods at capturing the true values of acoustic properties.
One key factor that appeared to play a role in the intracranial pressure map between the two VHP skull models, and helps account for the variation in pressure maps in the lager dataset, is the size of the skull. Comparison of the VHP female and male skull models, and their baseless counterparts, suggest that smaller skulls interact with the main acoustic beam more, as further illustrated by our investigations with a primate skull as well. This is an important phenomena to consider when considering neuromodulatory ultrasound applied in small animal models, such as a rat model where the skull is substantially smaller than the ultrasound transducer. Previous numerical simulations of experimental setups with rats have found substantial intracranial pressure increases in rat skulls compared to the corresponding transducer's free-field pressures for an acoustic beam larger than the rat's skull [40]. Compared to rats, the skulls of humans and nonhuman primates are much larger, and these larger skull sizes likely factor into why we observed a decrease in intracranial acoustic pressure compared to the baseless skull models in this work. Furthermore, previous numerical simulations of experimental setups with primates have found differences in peak pressures and FWHM volumes between two differing monkey skull models, one female and the other male, with corresponding differences in weight and skull thickness [41]. This further supports our evidence of varying max intracranial pressures between individuals. Overall it is important to keep in mind that the greater ultrasound wave confinement in smaller skull models will result in stronger and more numerous interactions, and that wave interactions can be either constructive or destructive, dependent on the phase of intersecting waves. Knowing which type of interaction will happen is complicated to predict, especially for unideal geometries.
Due to the lower pressures produced by our models of neuromodulatory ultrasound, we observed virtually no difference between linear and nonlinear models, as expected. This held for maximal pressures, FWHM volumes, and the locations of intracranial max pressure. However, for models and applications utilizing higher pressures, or possibly small skull models where reflections interact to produce elevated pressures, there will be a point when linear acoustic theory no longer represents the wave propagation of the system accurately. In one computational study, it was found that the effect of nonlinear distortion on temperature rise in tissue ranged from negligible at 1 MPa source pressures to an 80% increase in temperature elevation at 10 MPa source pressures [42]. Another study using empirical temperature measurements in an agar phantom with a 1 MHz focused transducer similarly supports that the threshold source pressure for significant nonlinear heating effects is over 1 MPa [43]. These results reinforce that at source pressures around 1 MPa, wave propagation is well described by linear acoustic theory, and that applications utilizing higher pressures may need to consider nonlinear effects.
Regarding our model population's average age of 62 years, it is unlikely that the age of our model population could largely bias results due to changes of the skull with age. Previous studies have only shown a general and insignificant trend of increased skull thickness with age, although there is a statistically significant relationship between cortical thinning and age for the inner and outer tables of the skull in females [44]. Similarly, regarding skull density, male skull bone density remains constant with age while female skull bone density decays slowly starting at 20 years of age [45]. Given that the majority of our modeling population is male, we believe that the variations in intracranial pressure maps observed between individuals are largely due to the morphological differences in skulls between individuals.
A technology similar to tFUS that has noted variability but successful application is transcranial Doppler sonography. Transcranial Doppler sonography is used to evaluate cerebral blood flow velocities, and its applications include detection of cerebral emboli and prognosis in acute stroke. The successful use of transcranial Doppler sonography is dependent on placement of the typically 2 MHz probe over an adequate acoustic window; where the transmittance of 2 MHz ultrasound has been reported to be only 0.8% to 13.2% depending on temporal bone thickness [46]. The rate of inadequate temporal window is around 10% for transcranial Doppler sonography, with a higher prevalence in females than males and a confirmed influence of aging [47,48]. Since previously mentioned studies found no significant differences in skull thickness or density in relation to age in the male skull, the finding of inadequate temporal window in males in transcranial Doppler sonography, albeit at a low rate, indicates that there may be some factor beyond the composition and thickness of the skull affecting performance and leading to variability in successful application. Possible factors include the acoustic frequency of ultrasound, the placement of the transducer on the skull with respect to skull morphology, and the size of the ultrasound aperture.
There have been a number of additional reports in literature on the variability of the transitivity of the skull for transcranial ultrasound dependent on both the acoustic frequency as well as the placement of the ultrasound transducer relative to the skull. Sonothrombolysis is an actively studied method with recent clinical trials where irradiation with ultrasound is used to enhance the thrombolytic activity of recombinant tissue plasminogen activator specifically at the ultrasound irradiation site [49][50][51][52]. For this application, estimation of transmittance through the human skull is an important factor for system design and dosing calculation. If transmittance is unexpectedly high or low, the acoustic intensity in the brain will be correspondingly higher or lower than anticipated, and the risk of cerebral hemorrhage will be elevated or the effect of thrombolysis reduced. Nonetheless, the correct estimation of skull transmittance has proved difficult due to the complex properties of the human skull. Frequency dependent fluctuation of ultrasound transmittance in the skull has been found in calculations and measurements using large-aperture multiple-source transducer arrays [36]. Further complicating matters, transmittance has been reported to fluctuate with placement of the ultrasound source relative to the skull [38,39], and has even been noted to occur in multi-frequency transducer design [41]. Additionally, tests on the sensitivity of simplified transcranial models to the variation of acoustic parameters have been found to influence the pressure amplitude readily, while the site of peak pressure and the FWHM volume are more stable [25]. Despite the sensitivity of representation of the skull in computational models, as the ultrasound source and the equations used for definition of skull acoustic properties were held consistent across simulations in this study, the results obtained here easily allow the comparison of the effects of individual differences in skull morphology. These findings, as well as those of previous literature, underscore the importance of detailed models of the skull and accurate representation of the ultrasound source and its placement relative to the skull for precise estimation of ultrasound dosing.
An additional important observation is the prominence of standing waves found in simulations at 250 kHz. Using frequencies below 500 kHz can reduce the aberration and attenuation of incident ultrasound by the skull, but comes at the cost of a larger focal spot size and a greater influence of standing waves. These standing waves may become a competing focal spot of ultrasound effects in the context of neural modulation, or a possible site of damage. The TRUMBI (Transcranial Low-Frequency Ultrasound-Mediated Thrombolysis in Brain Ischemia) clinical trial evaluating sonothrombolysis using 300 kHz ultrasound was stopped prematurely due to a high ratio of cerebral hemorrhages [49]. The exact reason for elevated cerebral hemorrhages was not elucidated in the TRUMBI clinical trial, but it was speculated to be due to locally intense spots due to standing waves stemming from the use of lower acoustic frequencies of ultrasound [49,53]. Baron et al demonstrated numerically the likelihood of inducing standing waves with the low frequency TRUMPI setup [54], in agreement with our observations in simulations at 250 kHz.
We simulated transcranial focused ultrasound for noninvasive neuromodulation using detailed skull models derived from medical images to evaluate simulation methods and the variation of intracranial pressure maps between individuals. The exact location and dosing of transcranial neuromodulatory ultrasound was found to vary between individuals due to differences in skull morphology and composition, but can be readily estimated for an individual using computational methods paired with medical scans of the head. The methods presented here can help to ensure proper targeting and to more tightly control dosing, which will account for and reduce variability in research optimizing low intensity applications of transcranial focused ultrasound.