# A simulation study on the sensitivity of transcranial ray-tracing ultrasound modeling to skull properties

Robert Andrew Drainville, Sylvain Chatillon, David Moore, John Snell,

Frédéric Padilla, Cyril Lafon

► To cite this version:

Robert Andrew Drainville, Sylvain Chatillon, David Moore, John Snell, Frédéric Padilla, et al. A simulation study on the sensitivity of transcranial ray-tracing ultrasound modeling to skull properties. Journal of the Acoustical Society of America, 2023, 154 (2), pp.1211-1225. 10.1121/10.0020761. cea-04241156

HAL Id: cea-04241156

https://cea.hal.science/cea-04241156

Submitted on 13 Oct 2023

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
# A simulation study on the sensitivity of transcranial ray-tracing ultrasound modeling to skull properties

Robert Andrew Drainville ID; Sylvain Chatillon ID; David Moore; John Snell ID; Frederic Padilla ID;
Cyril Lafon ID

<figure><img src="image_2.png" /></figure>

Check for updates

J. Acoust. Soc. Am. 154, 1211–1225 (2023)

https://doi.org/10.1121/10.0020761

<figure><img src="image_3.png" /></figure>

View
Online

<figure><img src="image_4.png" /></figure>

Export
Citation

CrossMark

## Related Content

Simulation of ultrasound propagation through human skull: Experimental validation and application to treatment planning

J Acoust Soc Am (September 2018)

Performance of a simulation-based phase aberration correction technique in transcranial ultrasound modeling

J Acoust Soc Am (April 2016)

Acoustic and thermal numerical modelling of transcranial ultrasound thermal therapy

J Acoust Soc Am (October 2021)

## Take the Lead in Acoustics

The ability to account for coupled physics phenomena lets you predict, optimize, and virtually test a design under real-world conditions – even before a first prototype is built.

» Learn more about COMSOL Multiphysics®

<figure><img src="image_5.png" /></figure>

<figure><img src="image_6.png" /></figure>

COMSOL
# A simulation study on the sensitivity of transcranial ray-tracing ultrasound modeling to skull properties

Robert Andrew Drainville,<fnref n="1" />,ª ID Sylvain Chatillon,<fnref n="2" /> ID David Moore,<fnref n="3" /> John Snell,<fnref n="4" /> ID

Frederic Padilla,<fnref n="3" />,ᵇ ID and Cyril Lafon<fnref n="1" /> ID

## ABSTRACT:

In transcranial focused ultrasound therapies, such as treating essential tremor via thermal ablation in the thalamus, acoustic energy is focused through the skull using a phased-array transducer. Ray tracing is a computationally efficient method that can correct skull-induced phase aberrations *via* per-element phase delay calculations using patient-specific computed tomography (CT) data. However, recent studies show that variations in CT-derived Hounsfield unit may account for only 50% of the speed of sound variability in human skull specimens, potentially limiting clinical transcranial ultrasound applications. Therefore, understanding the sensitivity of treatment planning methods to material parameter variations is essential. The present work uses a ray-tracing simulation model to explore how imprecision in model inputs, arising from clinically significant uncertainties in skull properties or considerations of acoustic phenomena, affects acoustic focusing quality through the skull. We propose and validate new methods to optimize ray-tracing skull simulations for clinical treatment planning, relevant for predicting intracranial target's thermal rise, using experimental data from ex-vivo human skulls.

© 2023 Author(s). All article content, except where otherwise noted, is licensed under a Creative Commons Attribution (CC BY) license (http://creativecommons.org/licenses/by/4.0/). https://doi.org/10.1121/10.0020761

(Received 14 March 2023; revised 23 June 2023; accepted 8 August 2023; published online 23 August 2023)

[Editor: Bradley E. Treeby]

Pages: 1211–1225

## I. INTRODUCTION

Ultrasound has long been recognized as a therapeutic tool with the potential to treat disorders throughout the brain. However, the distortion and attenuation of acoustic waves caused by the presence of the skull have been significant challenges to the development and clinical application of transcranial ultrasound therapies.<fnref n="1" />,<fnref n="2" /> Acoustic waves experience significant distortion during propagation through the skull due to reflection, refraction, scattering, absorption, and mode conversion.<fnref n="3" /> Early attempts to use focused ultrasound in the brain required craniotomy for the successful delivery of energy to the brain tissues.<fnref n="4" />,<fnref n="5" /> The development of large phased-array transducers and high-performance computers, however, has enabled the calculation of phase corrections necessary to focus ultrasound waves through the skull at the intended target, making it possible to obtain quality acoustic focusing through the skull and enable the clinical application of transcranial focused ultrasound.⁶⁻⁹

Various methods exist to account for the distortion caused by the presence of the skull in transcranial focused ultrasound applications. The gold standard for transcranial

acoustic refocusing is the hydrophone-based time reversal approach, in which a hydrophone is implanted at the location of the desired focus to directly measure the phase shifts induced by the presence of the skull at the transducer location, which are then inverted and applied across all elements to ensure the emitted waves converge in phase at the focal spot and maximize the acoustic energy at the target.⁶ Time-reversal modeling based on the reciprocity principle assumes that forward and backward ultrasonic wave propagation has the same time-frequency response and that an acoustic signal that propagates from the desired focal point to the ultrasound transducer can be recorded, reversed, and re-emitted at the transducer to obtain optimal focusing at the target.¹⁰

The development of numerical model-based approaches that calculate the phase corrections necessary to compensate for skull-induced phase aberrations using computed tomography (CT) scans of individual patient skulls has enable completely non-invasive transcranial focusing. These models consider CT or magnetic resonance (MR) images of the skull from which the density and acoustic velocity in the skull are inferred and the position of the transducer can be registered.<fnref n="1" />,¹¹⁻¹³ These time-reversal numerical simulations may use full wave models, which employ a three-dimensional (3D) finite-difference time domain (FDTD)

<figure><img src="image_7.png" /></figure>

<footnote n="1">¹LabTAU, INSERM, Centre Léon Bérard, Université Lyon 1, Univ Lyon, F-69003, Lyon, France</footnote>

<footnote n="2">²Université Paris-Saclay, CEA, List, F-91120, Palaiseau, France</footnote>

<footnote n="3">³Focused Ultrasound Foundation, Charlottesville, Virginia 22903, USA</footnote>

<footnote n="4">⁴Histosonics, Ann Arbor, Michigan 48103, USA</footnote>

<footnote n="5">⁵University of Virginia, Charlottesville, Virginia 22904, USA</footnote>

<footnote>ªEmail: radrainv@lakeheadu.ca</footnote>

<footnote>ᵇAlso at: University of Virginia, Charlottesville, Virginia 22904, USA.</footnote>
code to solve the acoustic or viscoelastic wave equation,¹⁴⁻¹⁶ or a ray tracing model which calculates the geometric path and associated propagation time from source to target using Snell's law and acoustic properties derived from CT scans.⁶

The use of a homogeneous model of the skull has been validated for time-reversal ray-tracing modeling of transcranial propagation at 710 kHz.¹⁷ Ray tracing methods based on Snell-Descartes law consider the overall geometry of the skull, but typically not the internal fine structures,¹⁷ with the skull often being considered as a three-layer medium composed of the outer table bone, the middle cancellous layer, and the inner table bone.³,¹⁸ While full wave simulations typically provide superior results compared to ray tracing models, their implementation has been limited due to its significantly higher computational burden.¹⁹ Bancel et al.²⁰ compared the performance of phase corrections derived from the ray-tracing algorithm and those from full wave simulations to those derived from hydrophone-based correction, which found that full wave models and the ray-tracing algorithm performed similarly in restoring 86% and 84% of the maximum focal pressure relative to hydrophone-based correction, respectively.

Material properties that impact the simulations include bone geometry and the acoustic properties of density, longitudinal and shear wave speeds, and attenuation, and the precise level of fidelity required to obtain accurate simulation results remains unclear.²¹,²² A previous study using full-wave simulations of a single-element transducer and a heterogeneous skull model found the longitudinal speed of sound to be the most influential acoustic parameter for producing effective focusing, that should be predicted with an accuracy of 4% to obtain clinically-acceptable focalization quality in terms of peak pressure, focus position, and focal volume estimations.²¹ Multiple different conversion algorithms have been proposed to allow for the conversion from Hounsfield units (HU) obtained from CT images to acoustic velocity for phase correction through the skull.⁷,⁹,¹³,²¹,²³⁻²⁹ While there is some consensus between these models, significant disagreement remains on the value of acoustic velocity inferred based on the same HU values between the models. To further complicate the conversion of CT images to skull properties, there exists no consensus on the relationship between acoustic attenuation and HU at clinically-relevant ultrasound frequencies.²³

The estimation of the longitudinal speed of sound in the skull from HU is influenced by CT parameters, including the x-ray energy, reconstruction method, and reconstruction kernel.²⁴ An assessment of how these parameters affect the longitudinal speed of sound estimation in human skulls found that CT is able to account for 23%–53% of the variation. While accounting for these variations in skull property reconstruction algorithms can improve the accuracy of phase correction and improve the transmission of acoustic energy through the skull, the remaining uncertainties in material properties conversion may give rise to errors or selection of treatment parameters that are not

optimized, and affect clinical applications of transcranial focused ultrasound.

Changes to the intracranial ultrasound field due to the presence of shear waves within the skull are generally considered negligible when the ultrasound is applied at normal or near-normal incidences, as is the case with large hemispherical arrays.²¹,²³,³⁰ Previous examinations²¹ have considered the presence of shear waves within simulations only to verify the validity of neglecting them and found the contribution of shear waves to be negligible for frequencies of 500, 750, and 1000 kHz.

Previous work by Robertson et al.²¹ has presented an examination of the sensitivity of simulated transcranial ultrasound fields to uncertainties in the skull properties using the full wave open-source k-Wave model, which was substantiated using experimental comparison with ultrasound transmission through 3D printed and resin cast skull bone phantoms and a single element transducer. However, to date, no systematic parametric study has been presented that examines the sensitivity of transcranial ultrasound focusing on the acoustic properties of the skull using a ray-tracing model. In this study, we examine the use of the ray-tracing model as a tool for treatment planning to predict temperature rise and temperature distribution in the brain when considering a 1000-element clinical transducer system and complete human skull model, including how transcranial ultrasound focusing is influenced by imprecision in the model inputs, such as those that may arise from clinically significant uncertainties in skull properties or the consideration of acoustic phenomena such as shear mode conversion and internal reflection. The goal of this work is to examine the ray-tracing model as a method of clinical treatment planning for therapeutic ultrasound that can be used to predict the focal acoustic field and temperature rise in brain tissue while identifying the most influential parameters to provide a framework for enhancing simulation applicability in clinical scenarios. The simulation sensitivity analysis is then substantiated through a comparison of experimental sonication experiments with simulations for the prediction of both the acoustic field and thermal rise.

## II. METHODS

Simulations were conducted using the CIVA Healthcare platform,³¹ an expertise platform for conducting non-destructive testing simulations that incorporates unique tools specifically designed for the simulation of 3D ultrasonic wave propagation in geometrically-complex isotropic and anisotropic media generated by arbitrary transducers and the subsequent receiver response originating from their interactions with defects and boundaries. The field radiation model is based on the so-called “pencil model,”³²⁻³⁵ which consists in discretizing the radiating surface into a set of independent source points whose elementary contributions are summed at each calculation point to obtain the impulse response of the radiated field. Each contribution is evaluated along an ultrasonic ray-path between the source point and the
calculation point. This path is obtained by applying the Snell's laws at each interface, in transmission or reflection, with or without considering mode conversions. The amplitude of each elementary contribution is computed by a high frequency asymptotic formulation close to paraxial ray models.³⁶ The main approximation of this high frequency model consists in neglecting frequency dependent effects such as edge diffraction by sharp edges of interfaces, surface wave propagation along interfaces or head wave propagation beyond critical angles. In the present application, the description of the skull surface is regular and smooth and the angles of incidence are below critical angles. Therefore, the use of this type of paraxial ray model is appropriate, and very satisfactory results can be obtained. Recent work has presented benchmarking of different transcranial ultrasound simulation models³⁷ to examine whether different modeling techniques provide the same answer when the inputs to the model are well specified, and an extension of this work is underway to provide the same benchmark testing using the CIVA Healthcare ray-tracing algorithm.

Numerical simulations considered human skull samples, which were acquired in accordance with standard institutional policy for human tissue. The skull had previously been debrided of all soft tissue and degassed using an industrial vacuum chamber with degassed water. A 3D layered model was generated from CT scans as shown in Fig. 1, with a CT spatial resolution of 1 mm or less along each axis and a tube voltage of 120 kV. The skull was positioned relative to the transducer to be similar to transcranial clinical treatment of essential tremor,³⁸ with the acoustic focal point located at the geometric focus of the transducer. The anterior, superior, and left directions of the patient are aligned with the negative y-axis, positive z-axis, and positive x-axis, respectively.

The transducer source was modeled on a transcranial InSightec ExAblate Neuro (Insightec, Haifa, Israel) clinical system, which offers the ability to focus acoustic energy by adjusting the amplitude and phase of each of its 1024 transducer elements to compensate for the distortions introduced

<caption>FIG. 1. (Color online) 3D skull CAD model file used for simulating ultrasound field. Volume where field is calculated is shown as a cube at the geometric focus of the transducer. The anterior, superior, and left directions of the patient are aligned with the negative y-axis, positive z-axis, and positive x-axis, respectively.</caption>

by transmission through the skull while using the maximum available surface area and chilled degassed water is circulated around the head to cool the scalp and reduce skull heating.³⁹⁻⁴¹ The use of the ExAblate system has been approved for use in the United States by the Food and Drug Administration (FDA) for use in unilateral thalamotomy for the treatment of essential tremor.³⁹ The hemispherical transducer has a diameter of 30 cm, and simulations considered continuous wave sonications at 670 kHz. A model of the proprietary site-specific Exablate transducer was derived based on the location of the center of each element and their corresponding area.

Simulations were performed using phase delay profiles calculated by the built-in CIVA phase correction software, which is able to provide near-instantaneous phase delay correction profile for a given material configuration by calculating the geometric acoustic path between the intended focal point and each element. Phase delay calculations take into account energy transmission and the time of flight from each element to the focal point. Any elements for which a ray path cannot be found from the target location to the emitting surface are deactivated.

For each set of simulations, the focal quality was quantified using the value of the maximum amplitude ($A_{max}$), and the distance from the position at which the maximum amplitude occurs to the intended target at the geometric focus of the transducer, which is referred to as the focal offset ($r_{max} = |\vec{r}_{max}|$). The focal volume ($V_f$) was calculated based on the total volume of the largest continuous region with pressure greater than 50% of the maximum pressure. When calculated, the focal shift ($\Delta r$) was quantified as the change in the location of the maximum amplitude between different simulations induced by changes in simulation conditions, such as comparing how the inclusion of phenomena such as mode conversion or internal reflection affects focusing and was calculated according to

$$ \Delta r_i = |\vec{r}_{max,0} - \vec{r}_{max,i}|, \qquad (1) $$

where $\vec{r}_{max,i}$ is the position of the maximum acoustic amplitude for a given simulation, and $\vec{r}_{max,0}$ is the position of the maximum acoustic amplitude obtained in the reference simulation.

## A. Assessing focal quality with homogeneous skull model simulations

Simulations considering homogeneous models of two skull samples were used to examine the influence of variation in the overall skull properties. The homogeneous 3D models of the skull samples were obtained by segmenting CT data using the open source software package 3D Slicer (www.slicer.org).

Simulations were performed to calculate the acoustic field within a volume with dimensions 4 mm by 4 mm by 6 mm along x, y, and z, respectively, with an isotropic spatial resolution of 0.4 mm centered at the geometric focus of the transducer, as shown in Fig. 1. The calculated acoustic fields

<figure><img src="image_8.png" /></figure>
were then resampled with a spline interpolation of order three to increase the spatial resolution by a factor of 32 to precisely determine the maximum amplitude and its corresponding location.³⁷ All simulations used an acoustic power of 100 W.

# 1. Impact of longitudinal speed of sound conversion uncertainty

To examine uncertainty in the conversion of skull density to acoustic velocity, simulations were conducted at five values of density ranging from 1600 to 2400 kg/m³, consistent with the range of skull density values correlated with acoustic velocity by Leung et al.²³ For each skull density value, simulations were performed at three longitudinal wave velocity values, which were determined by examining the various density/velocity relations also provided by Leung et al.,²³ identifying the minimum ($c_{min}$) and maximum ($c_{max}$) predicted velocity values for each value of density across all conversion algorithms, and taking the average between the two extreme values ($c_{avg}$). The skull density values and their corresponding longitudinal wave velocities are presented in Table I. As has been done in previous simulation studies,²¹ the shear wave speed remained constant at 1500 m/s³⁰,⁴² for all simulations, regardless of skull density.

For each density value, the phase correction profile was calculated based on the average speed of sound, $c_{avg}(\rho)$, and simulations were performed at all three longitudinal speed of sound values. Thus, the simulations performed with a skull speed of sound set to $c_{avg}$, provide a case in which there is no mismatch between the skull speed of sound used for simulation and for calculations of the phase profile, against which the results from the other two simulations can be compared.

The experimentally measured longitudinal-wave attenuation coefficient for the whole skull, given by $\alpha_0$, has previously been estimated²¹,⁴³ to be given by

$$ \alpha_0 = 8.83 \frac{\text{dB}}{\text{cm} \cdot \text{MHz}^{1.43}}, \qquad (2) $$

which gives an attenuation coefficient of 0.498 dB/mm at the frequency of 670 kHz used in simulations. Shear wave mode conversion was included for all simulations, with longitudinal to shear wave conversion occurring at the fluid/bone interface, and shear to longitudinal conversion occurring at the bone/fluid interface. The shear wave attenuation values were set equal to twice the value used for

<table><caption>TABLE I. Range of skull density and corresponding longitudinal speed of sound values used for simulations.</caption><thead><tr><th>Density (kg/m³)</th><th>c<sub>min</sub>(m/s)</th><th>c<sub>avg</sub>(m/s)</th><th>c<sub>max</sub>(m/s)</th></tr></thead><tbody><tr><td>1600</td><td>2219</td><td>2260</td><td>2300</td></tr><tr><td>1800</td><td>2343</td><td>2431</td><td>2518</td></tr><tr><td>2000</td><td>2495</td><td>2634</td><td>2772</td></tr><tr><td>2200</td><td>2689</td><td>2858</td><td>3027</td></tr><tr><td>2400</td><td>2939</td><td>3238</td><td>3536</td></tr></tbody></table>

longitudinal wave attenuation, which is approximately what is observed in solids.²⁸

# 2. Independent longitudinal speed of sound variation in homogeneous skull models

To assess the impact of uncertainty in the acoustic velocity of the skull independent of skull density, simulations were performed in which the acoustic velocity was varied while other properties were held constant. The acoustic velocity varied from 2200 m/s to 3600 m/s - equivalent to the range of values shown in Table I. The phase delay profile was calculated with a skull density and longitudinal speed of sound value of 2000 kg/m³ and 2634 m/s, respectively, which are similar to average skull values used in previous simulations.²¹ Results from the simulation where the skull density and longitudinal speed of sound are not changed from these values can be considered as the reference, as there is no mismatch between the skull properties used for simulation and those used for the calculation of the phase profile. Simulations were performed with the longitudinal attenuation constant set to 0.498 dB/mm.

# 3. Independent attenuation variation in homogeneous skull models

To determine the effect of uncertainty in the acoustic attenuation in the skull independent of other parameters, simulations were performed in which the attenuation was varied with all other acoustic properties held constant. The phase delay profile was calculated with a skull density and acoustic velocity of 2000 kg/m³ and 2634 m/s, respectively. Attenuation coefficient values between 0 and 1.5 dB/mm were considered to account for the attenuation uncertainty in the literature.²³

# B. Layered skull model: Examining skull heterogeneity, internal reflection, and shear wave mode conversion

To examine the importance of considering skull bone micro-architecture, simulations were performed in which skull sample #1 was modeled as a homogeneous medium possessing properties derived from those of the whole skull, and compared to those in which the skull was modeled as a tri-layered material with properties determined by the characteristics of cortical and trabecular bone. The models of the cortical and trabecular layers were generated using the open source software platform 3D Slicer with CT data from the human skull sample. In order to generate the 3D models of the cortical bone layer, any voxels which exhibited a HU value greater than 250 were included as part of the 3D model of the cortical bone volume, and any cavities within the bone layer were filled to create a single solid layer. To generate the various trabecular bone layers, voxels below a specified HU threshold value that overlapped with the cortical model were used to create different 3D models of the trabecular bone, which could be composed of multiple separate islands and allow for internal cavities where the bone was
dense enough to be considered as cortical. The specific trabecular threshold value was varied to create multiple different models of trabecular bone. The various 3D models of the trabecular bone layer could then be overlaid on the single cortical layer model to simulate layered skull models with varying proportions of trabecular bone. The generation of skull models with different quantities of trabecular bone was quantified using the volume of trabecular bone relative to the total skull volume, which is referred to as the trabecular fraction.

A cross section of the CT scan of the skull and an example of the cortical and trabecular bone layer models are shown in Fig. 2. The volume of the trabecular bone relative to the total skull volume is shown in Fig. 2(c) as a function of the trabecular HU threshold value, which demonstrates how the volume of the trabecular varies with the maximum HU threshold value used to generate the 3D models. The trabecular volume relative to the whole skull ranged from a minimum value of 0% for a skull composed purely of cortical bone, to a maximum value of 55.7% at a HU threshold of 1300. The heterogeneous skull model with 0 trabecular fraction was assigned the properties of cortical bone, while the homogeneous skull model, included for comparison, used the same CAD file as the cortical bone but was assigned acoustic properties based on the average skull properties.

The CIVA Healthcare platform provides the user-determined option to consider multiple internal reflections in the calculation of the intracranial acoustic field, which can approximate the effects of acoustic resonance within different layers inside the skull. In order to isolate and examine the effects of internal reflection on transcranial ultrasound propagation in a heterogeneous skull model, results from simulations that included reflection from boundaries were compared to simulations in which they were neglected. The number of internal reflections considered within CIVA Healthcare ray-tracing simulations must be specified within the user-defined simulation parameters. The number of internal reflections was limited to two, as the exponential decrease in amplitude with further reflections was expected to make any additional contributions negligible. Due to the large increase in computational time required to consider both internal reflection and mode conversion, simulations considering internal reflection did not include the contribution arising from shear wave mode conversion within the skull.

To investigate the effects of mode conversion with a layered skull model, simulations that included the contribution of shear mode conversion were compared with simulations in which shear mode conversion was neglected. As with internal reflections, the number of mode conversions

<figure><img src="image_9.png" /><figcaption>FIG. 2. (Color online) (a) Cross section of the CT image of skull. (b) A corresponding cross section of the cortical (green) and trabecular (yellow) models derived from CT data. Illustrated trabecular model uses 1200HU as maximum HU threshold value. (c) Plot of the volume of cortical bone relative to the whole skull volume as a function of the maximum HU threshold for trabecular bone.</figcaption></figure>
considered within CIVA simulations must be specified by the user. Based on the increased degree of computational resources required to consider the possible combinations of mode conversions at different boundaries, including bone-water and cortical-trabecular, the number of mode conversions considered was limited to two, to allow for the possibility of longitudinal-transverse-longitudinal waves.

### C. Experimental verification of simulation outcomes

Measurements were performed through a human skull to mimic the clinical treatment of essential tremor.⁴⁴ Acquisition of cadaveric material was performed in accordance with the standard institutional policy for human tissue. Skulls were debrided of all soft tissue and degassed using an industrial vacuum chamber filled with degassed water.

### 1. Assessing acoustic field simulation accuracy

To provide a comparison between experimentally measured acoustic fields and the simulated acoustic fields obtained using the CIVA platform, simulations were performed to obtain the pressure field at the focus in the absence of a skull, as well as when the beam was steered along each axis, with results being compared to the experimental acoustic field as measured used a hydrophone under similar conditions. Simulations were performed with the transducer focused at the geometric focus as well as offset 10 mm along each axis, with experiments being performed under the same conditions. Simulations were then performed in which the human skull sample was placed at the geometric focus of the transducer, which were then compared to experimental measurements.

For acoustic field pressure measurements, the skulls were mounted inside the ExAblate transducer, in a similar positioning as would the head of a patient for a clinical treatment. The skull was maintained in place using a rigid acrylic frame, with the transducer positioned vertically facing in the upward direction and with a water retainer mounted on top of it to secure a water column. An MR scan was acquired, and together with pre-acquired CT scan of the skull, was loaded into the ExAblate software to compute phase aberration correction for each of the active elements, as previously described.⁴⁴ For pressure field characterization, a needle hydrophone (HNA-0400 needle hydrophone, Onda Corp., Sunnyvale, CA) was mounted on computer-controlled 3D positioning mechanical stages to scan a volume centered around the geometrical focus of the transducer. The acoustic signals were acquired with an oscilloscope (DSO7012B, Keysight Technologies) sampling at 100 MHz, and then converted in MPa, after deconvolution from the impulse response for the hydrophone. Acoustic power was set to 1.3 W to maintain pressure levels within the safe operating range of the needle hydrophone.

### 2. Evaluating temperature rise predictions in brain tissue phantom

The widely-used bio-heat transfer equation (BHTE) provides a method to calculate temperature increases within tissue due to ultrasound.⁴⁵,⁴⁶ The general form of the bioheat equation is given by⁴⁷

$$ \rho c \frac{\partial T}{\partial t} = \nabla(k\nabla T) + \omega_b c_b \rho_b (T_a - T) + q_m + Q, \quad (3) $$

where $ρ, c$, and $k$ are the density (kg/m³), the specific heat (J kg⁻¹ K⁻¹), and the tissue thermal conductivity. The second term on the right hand side of the equation represents changes due to blood perfusion, where $w_b$ is the mass flow rate of blood per unit volume of tissue (kg s⁻¹ m⁻³), $c_b$ is the blood specific heat, $q_m$ is the metabolic heat generation (J s⁻¹ m⁻³), and $(T_a - T)$ represents the difference between local tissue temperature and temperature of arterial blood. The final term $Q$ represents the heat source contribution, which in the current context is due to the absorption of acoustic energy. The heating source is related to the acoustic field according to⁴⁸

$$ Q = \frac{\alpha_l p^2}{\rho c_l}, \quad (4) $$

where $α_l$ is the acoustic absorption coefficient for the longitudinal wave (Np/m), $p$ is the amplitude of the acoustic pressure wave (Pa), and $c_l$ is the longitudinal speed of sound. Given the use of a homogeneous thermometric tissue phantom, attenuation was assumed to be entirely attributable to absorption.

A numerical BHTE model was coupled with the acoustic CIVA ray-tracing model to simulate the temperature rise produced at the acoustic focus. Acoustic simulations were performed with a spatial resolution of 0.5 mm along each axis, with domain dimensions of 16 × 16 × 21 mm along the x-, y-, and z-axes, respectively, and with a time step calculated to provide a Courant-Friedrieichs-Lewy number of 0.05. After the acoustic field had been calculated, the acoustic field was interpolated to increase the spatial resolution by a factor of 4. The temperature was fixed at the domain boundaries, and comparison with perfectly insulating boundaries did not show any meaningful differences in the resultant temperature distribution. The finite difference method employed an O(2,4) scheme for calculation of spatial and temporal derivatives, which is 2nd order in time and 4th order in space, as the use of higher-order schemes for the calculation of spatial derivatives has been found to produce lower numerical errors.⁴⁹,⁵⁰

For thermometry measurements, a slab of tissue mimicking material was cast inside a human skull, and positioned within the transducer so that the geometrical focus of the transducer was positioned 10 mm deep into the tissue mimicking material (TMM). The TMM was manufactured by CIRS (Norfolk, VA), with the following properties: density 1.04 g/cm³, speed of sound 1540 m/s, attenuation
coefficient 0.50 dB/cm/MHz. Sonications of 10-s durations were performed with acoustic power levels ranging from 50 to 325 W. Temperature maps were obtained using PRF-based MR thermometry.⁵¹

The properties of the TMM used in simulations were provided by the manufacturer (CIRS, Norfolk, VA), and are given in Table II. The skull properties were estimated using average whole-skull CT data and the conversion algorithm provided by McDannold *et al.*,⁵² and are also shown in Table II.

The experimentally measured temperature rise was then compared against the results obtained from numerical simulations, in terms of the maximum temperature achieved and the spatial temperature distribution. In order to account for the finite resolution provided by MR thermometry, a spatial averaging filter was applied to numerical results, which provided spatially averaged temperatures equivalent to a resolution of 1 × 1 × 3 mm along the x-, y-, and z-axes, respectively.

### III. RESULTS

#### A. Assessing focal quality with homogeneous skull model simulations

##### 1. Impact of longitudinal speed of sound conversion uncertainty

The results of the simulations quantify how existing uncertainty in the acoustic properties of the skull determined from CT scans may manifest in reduced focal quality for clinical applications. The disagreement among proposed relationships between density and acoustic velocity produces differences in the peak pressure of almost 30% when compared against simulations where material property uncertainty is not present, as shown in Fig. 3(a).

Although the uncertainty in the predicted speed of sound can vary by as much as ±20%, the resultant error in spatial positioning remains relatively low, as shown in Fig. 3(b). The focal volume plotted in Fig. 3(c) shows a strong increase for simulations that use high values of skull density, particularly for simulations using a higher speed of sound than was used for the calculation of the phase correction profile.

<table><caption>TABLE II. Material properties used in acoustic and thermal simulations for tissue mimicking material (TMM) and human skull.</caption><thead><tr><th>Property</th><th>Units</th><th>TMM</th><th>Skull</th></tr></thead><tbody><tr><td>Thermal Conductivity (k)</td><td>J kg<sup>-1</sup> K<sup>-1</sup></td><td>0.5</td><td>NA</td></tr><tr><td>Longitudinal speed of sound (c<sub>l</sub>)</td><td>m s<sup>-1</sup></td><td>1540</td><td>2403</td></tr><tr><td>Shear speed of sound (c<sub>s</sub>)</td><td>m s<sup>-1</sup></td><td>0</td><td>1500</td></tr><tr><td>Density (ρ)</td><td>kg m<sup>-3</sup></td><td>1040</td><td>2212</td></tr><tr><td>Perfusion coefficient (w<sub>b</sub>)</td><td>kg s<sup>-1</sup> m<sup>-3</sup></td><td>0</td><td>NA</td></tr><tr><td>Longitudinal attenuation coefficient (α<sub>l</sub>)</td><td>Np m<sup>-1</sup></td><td>0.05</td><td>1.52</td></tr><tr><td>Shear attenuation coefficient (α<sub>s</sub>)</td><td>Np m<sup>-1</sup></td><td>NA</td><td>1.00</td></tr></tbody></table>

<figure><img src="image_10.png" /></figure>

<figure><img src="image_11.png" /></figure>

<figure><img src="image_12.png" /><figcaption>FIG. 3. (Color online) (a) Maximum amplitude vs density for three different longitudinal speed of sound values. (b) Focal offset vs density. (c) Focal volume vs density. Results obtained using skull 1 are shown with a solid line, while results from skull 2 are shown as a dashed line.</figcaption></figure>

#### 2. Independent longitudinal speed of sound variation in homogeneous skull models

The effect of independently varying the longitudinal velocity of the skull on maximum amplitude is shown in Fig. 4(a). This figure illustrates that a shift from a typical average acoustic velocity of the skull of 2634 m/s to the
<figure><img src="image_13.png" /><figcaption>FIG. 4. (Color online) (a) Maximum amplitude vs longitudinal speed of sound in the skull. (b) Focal offset vs longitudinal speed of sound in the skull. (c) Focal volume vs density. Results obtained using skull 1 are shown in red, while results from skull 2 are shown in blue. Dashed line denotes the speed of sound value used for the calculation of phase correction.</figcaption></figure>

maximum acoustic velocity of 3600 m/s, which is within the range of values obtained for the outer layers of cortical bone,²¹ produces a 50% reduction in maximum amplitude with skull 1, and a 60% reduction with skull 2. The effect of

lowering the skull velocity showed an asymmetry in the effect on maximum amplitude, with skull 1 showing in increase in amplitude, while the maximum amplitude obtained with skull 2 showed a smaller decrease in maximum amplitude compared to when skull velocity was increased. Figure 4(c) shows the effect of variation of longitudinal wave speed of the skull on the focal volume, with little change in focal volume when speed of sound decreases, and a large increase in focal volume at higher speed of sound values in both skulls.

The effect of varying the acoustic velocity of the skull independent of other properties on the location at which the maximum acoustic amplitude occurs is shown in Fig. 4(b), which produced a maximum focal offset of 0.83 mm and 1.4 mm for skulls 1 and 2, respectively.

### 3. Independent longitudinal attenuation variation in homogeneous skull models

The effect of variations in skull attenuation on maximum amplitude, its corresponding distance from the geometric focus, and the focal volume are shown in Fig. 5. These results show an exponential dependence of the maximum amplitude on the attenuation coefficient of the skull, with a 50% variation in the skull attenuation from the average value of 0.498 dB/mm being found to produce an amplitude error of approximately 25%. Variation in attenuation did not produce a significant change in the focal spot location, leading to small values of focal offset and focal drift for both skull samples. The discrete changes evident in Fig. 5(b) are caused by the finite spatial resolution. Figure 5(c) shows little dependence of the focal volume on the longitudinal attenuation in the skull, with a 12% and 47% increase in skull 1 and 2, respectively.

### B. Layered skull model: Examining skull heterogeneity, internal reflection, and shear wave mode conversion

As shown in Fig. 6(a), when comparing simulations that consider mode conversion with those that neglect it in a heterogeneous model of the skull, the largest changes in maximum amplitude occurred for trabecular fraction values between 0.4 and 0.5, which demonstrates a reduction of almost 30% due to the inclusion of shear waves. Regardless of whether shear wave mode conversion or internal reflection are included in simulations, the maximum amplitude trends towards a local minimum between trabecular fraction values 0.3 and 0.4. Figure 6(c) shows the effect of an increased trabecular fraction on the focal volume, where similar to the effect on maximum amplitude, the variation causes a decrease in the spatial confinement of the focus initially, before decreasing and returning to near the values obtained for a homogeneous skull model.

Figure 6(b) shows the change in the position at which the maximum acoustic amplitude occurs, relative to the geometric center of the transducer. The peak position offset remains relatively low, with less than a 0.8 mm shift
<figure><img src="image_14.png" /></figure>

<figure><img src="image_15.png" /></figure>

<figure><img src="image_16.png" /><figcaption>FIG. 5. (Color online) (a) Maximum amplitude vs longitudinal attenuation coefficient. (b) Focal offset vs longitudinal attenuation coefficient. (c) Focal volume vs longitudinal attenuation coefficient. Results obtained using skull 1 are shown in red, while results from skull 2 are shown in blue.</figcaption></figure>

regardless of the amount of trabecular bone, though the positional change was greatest when shear waves are considered.

The focal shift of the acoustic focus, induced solely by the inclusion of internal reflection and calculated using

<figure><img src="image_17.png" /></figure>

<figure><img src="image_18.png" /></figure>

<figure><img src="image_19.png" /><figcaption>FIG. 6. (Color online) (a) Maximum acoustic amplitude at the focus as a function of skull trabecular bone fraction for simulations considering internal reflection (red), simulations considering shear wave conversion (blue), and simulations considering neither shear wave conversion nor internal reflection (green). Dashed lines show the results obtained using homogeneous skull profiles. (b) Distance from the location of maximum amplitude relative to the geometric center of the transducer, plotted as a function of the skull trabecular bone fraction.</figcaption></figure>

Eq. (1) using results obtained from simulations without internal reflection or shear wave mode conversion as the reference, is shown in Fig. 7. The additional contribution of internal reflection caused negligible shifts in the majority of
Focal Shift

<figure><img src="image_20.png" /><figcaption>FIG. 7. (Color online) Change in the position at which the maximum amplitude occurs induced by the simulation of internal reflection (red) and shear wave mode conversion (blue) at skull/water and cortical/trabecular boundaries, plotted as a function of the amount of trabecular bone relative to total skull volume. Dashed lines illustrate the value obtained with a homogeneous model of the skull.</figcaption></figure>

simulations, with only two results exhibiting a shift greater than 0.03 mm. The inclusion of internal reflections tended to cause a shift towards the geometric focus, which can be attributed to random variation and produced a smaller positional error compared to when the contributions due to internal reflection were neglected. The inclusion of shear waves tended to produce a slightly larger shift in the focal position, with a maximum value of approximately 0.1 mm, though the overall magnitude of the shift in all cases is considered negligible.

C. Experimental verification of simulation outcomes

1. Assessing acoustic field simulation accuracy

The simulated and experimentally measured acoustic field is presented in Fig. 8, and shows good agreement in both pressure magnitude and distribution. The maximum acoustic pressure was similar between simulation and experiment, with maximum amplitude values (MPa) in [simulation, experiment] being [1.84, 1.474], [1.535, 1.20], [1.54, 1.27], and [1.57, 1.25] when the focal spot was at the geometric center and displaced 10 mm along the x, y, and z axes, respectively. The simulated and experimentally measured acoustic field with transcranial focusing through skull sample #2 are shown in Fig. 9. The maximum acoustic pressure showed good agreement between simulation and experiment, while the spatial distribution of the simulated acoustic field showed better spatial confinement.

2. Evaluating temperature rise predictions in brain tissue phantom

The resultant acoustic field near the focus was found to have a focal volume of 6.96 mm³, and a focal offset of 0.269 mm. The maximum amplitude varied proportional to the square root of the total acoustic power, up to a maximum of 3.13 MPa for an acoustic power of 350 W.

<figure><img src="image_21.png" /><figcaption>FIG. 8. (Color online) Simulated (left) and experimentally measured (right) acoustic field when focusing in water without the presence of the skull at the geometric center of the transducer (first row), and displaced 10 cm from the geometric center 10 mm along x (second row), y (third row) and z (fourth row). Simulated fields are focused electronically using the CIVA-generated phase profile, while experiments employed ExAblate-derived phase profile. Simulated acoustic power was set to 1.3 W to be consistent with experimental conditions.</figcaption></figure>

A linear best fit model indicated that the experimental maximum temperature as a function of acoustic power, as shown in Fig. 10, rose at a rate of 0.0423 °C/W, while the simulated maximum temperature rose at a rate of 0.0599 °C/W. Once spatial averaging was applied to the simulated temperature data to provide the same spatial resolution as the experimental measurements, the maximum temperature rose
<figure><img src="image_22.png" /><figcaption>FIG. 9. (Color online) Experimental (left) and simulated (right) acoustic field at the focal plane while focusing through skull sample #2 with 19.5 W of acoustic power.</figcaption></figure>

at a rate of 0.0458 °C/W. The linear fit for the experimental temperature provided a temperature increase in 1.87 °C for zero acoustic power.

The temperature distribution profiles shown in Fig. 11 show good agreement between experimental and simulated results, though the simulated temperature rise is slightly lower than the experimentally measured values, and the application of spatial averaging to simulated temperature fields, further lowering the simulated temperature. The simulated fields show narrower temperature distributions, which may be attributed to the tighter focusing of the acoustic field, as demonstrated in Fig. 9.

# IV. DISCUSSION

## A. Assessing focal quality with homogeneous skull model simulations

### 1. Impact of longitudinal speed of sound conversion uncertainty

The change of amplitude with density and skull longitudinal speed of sound demonstrated a greater effect of change in skull density than the difference arising from variation in

the speed of sound conversion algorithms. The defocusing effect caused by a mismatch between the longitudinal speed of sound in the skull used in simulation versus the value used for the calculation of the phase profile has a smaller effect on the maximum amplitude than that caused by changes in the skull transmission due to changes in the skull acoustic impedance, leading the maximum amplitude to tend to occur in simulations with lower speed of sound at a given density. Simulations with both high skull density and longitudinal speed of sound produced the largest increase in focal volume, though this is also a reflection of the larger difference between $c_{min}$, $c_{avg}$, and $c_{max}$ at higher skull density values. The increase in focal volumes reinforces that the decrease in maximum amplitude is attributable to both decreased transmission through the skull due to increased acoustic impedance, defocusing effects due to velocity conversion uncertainty, and stronger distortion of the waves that impede tight acoustic focusing.

Figure 3(b) shows that the shift in focal position has little sensitivity to changes in density and longitudinal speed of sound variation, with all focal offset values being lower than 0.6 mm - approximately one quarter of the acoustic wavelength in water at 670 kHz. While the focal offset is low, it is interesting to note that the smallest focal offset did not always occur when the longitudinal speed of sound in the skull was set to $c_{avg}$, when there was no mismatch between the speed of sound used for the calculation of phase profile and the actual speed of sound used within the skull model. This may be attributed to the fact that the calculation of the phase profile maximizes the field amplitude at the defined focal point, rather than ensuring the position of maximum amplitude occurs at the intended focal point, though these two cases would ideally coincide. Though maximizing the amplitude at the focal point, the complex geometry of the skull and its effects on the acoustic field may inadvertently optimize the acoustic amplitude at the focus while still producing an acoustic field with a maximum that deviates from the precise intended focal point.

### 2. Independent longitudinal speed of sound variation in homogeneous skull models

The results shown in Sec. III A 2 quantify the effects of uncertainty of acoustic velocity in homogeneous models of the human skull on the quality of focusing. A previous examination²¹ of the sensitivity of transcranial ultrasound fields to variations in acoustic property maps observed that increasing acoustic velocity leads to lower peak pressure as well as some oscillatory behavior. The current results demonstrate that peak pressure exhibited some of the same oscillatory behavior in the peak pressure and decrease in focal pressure with increasing velocity, though the amount of oscillation was reduced, and a decrease in acoustic velocity eventually led to a reduction in the maximum focal amplitude. This may be attributed primarily to the fixed phase correction profile, the complex geometry of the human skull models used in simulations, and the calculation of a fixed phase profile based on the acoustic properties of the skull to

<figure><img src="image_23.png" /><figcaption>FIG. 10. (Color online) Experimental and simulated maximum temperature as a function of acoustic power.</figcaption></figure>
<figure><img src="image_24.png" /><figcaption>FIG. 11. (Color online) (a) Temperature as a function of lateral position (x/y) obtained experimentally and from simulations with 100 W of acoustic power. (b) Temperature as a function of lateral position (x/y) obtained experimentally and from simulations with 325 W of acoustic power.</figcaption></figure>

adapt focusing to the skull, whereas previous examinations used a simpler bone geometry and a single element transducer that did not adapt to bone geometry.

Existing literature suggests that the effects of refraction and phase delay may exert contradictory effects on the accuracy of the simulated pressure field within homogeneous models of the skull. Using k-Wave based simulations of transcranial ultrasound propagation with a planar phased array, Jiang *et al.*¹⁰ found that simplifying the skull bone as a homogeneous solid is more accurate for transcranial ultrasound path estimation when the properties of the skull bone are idealized as that of the cortical bone to accurately account for wave diffraction, as opposed to the properties averaged over the entire skull bone. These findings may suggest the usefulness of a weakly heterogeneous model of the skull that uses cortical bone properties at the water/bone interface to correctly model refractive effects and an internal layer with properties more similar to the average skull properties that properly account for time delay and attenuation. Further work examining the effect of considering bone micro-architecture on transcranial propagation and comparison with experimental results may provide further insight into how models of transcranial ultrasound propagation can be optimized for clinical applications.

### 3. Independent longitudinal attenuation variation in homogeneous skull models

The results shown in Sec. III A 2 quantify the effects of uncertainty of acoustic attenuation in homogeneous models of the human skull on the quality of focalization. The relationship shown in Fig. 5 indicates that variation of the attenuation coefficient of the human skull causes a corresponding variation in the maximum amplitude at the focus, with very limited effect on the focal position or focal volume.

Spatial variation in the acoustic attenuation has previously been found to have almost no effect on the intracranial

field,²¹ which is consistent with what is observed in the simulation results. This fact, combined with the smaller effects that uncertainties in attenuative properties have on focal quality, indicates that there should be a greater focus on properly modeling the spatial variations of acoustic properties of density and velocity and that precisely determining the relationship between density attenuation is of less importance to obtain high quality focusing than accurately determining the average attenuation over the whole skull.

### B. Layered skull model: Examining skull heterogeneity, internal reflection, and shear wave mode conversion

The results shown in Sec. III B demonstrate the effects that the degree of heterogeneity within the skull can have on simulated transcranial ultrasound focusing. Both the maximum amplitude and the spatial confinement exhibited an initial reduction as the trabecular fraction was increased from zero, with this effect largely reverting at higher values. This may be explained due to the presence of a more irregular and complex geometry of the trabecular bone at lower HU values that gives rise to more diffuse scattering.

The role of mode conversion was shown to be significant even with the homogeneous models of the skull, though the importance of considering mode conversion becomes more significant with greater degrees of heterogeneity.

When considering the overall effects of skull heterogeneity on transcranial focusing, mode conversion appears to have a much more significant effect on both amplitude and focal position than internal reflection between skull layers. Neglecting mode conversion in heterogeneous models caused an average amplitude change of 6.1%, while neglecting internal reflection caused an average amplitude change of only 2.3%. Similarly, the positional shift associated with neglecting internal reflection was on average 0.049 mm, which was significantly less than the 0.20 mm average
positional shift associated with mode conversion. The relatively large decrease in focal amplitude that occurred when shear waves were included with the homogeneous model of the skull can be attributed to destructive interference between the contributions from longitudinal waves and shear waves in the skull.

The role of internal reflection between different skull layers appeared to have little effect on the maximum focal acoustic amplitude or the spatial positioning of the focus. While some simulations indicated a tendency of the position at which the maximum amplitude occurs towards the geometric center of the transducer, this effect is largely negligible and may be attributed to stochastic distortion of the acoustic field arising from the complexity of the skull. The shift in spatial position of the focus due to internal reflection was smaller for the heterogeneous skull model than for the homogeneous model in most cases, indicating that internal reflection between skull layers plays a relatively small role in the focal quality.

The significance of internal reflection may have been minimized by the values of longitudinal speed of sound assumed for cortical and trabecular bone of 2900 and 2500 m/s, respectively. These numbers were obtained from published work by Clement et al.,⁵³ which found them to be optimal when using a tri-layered model of the skull to experimentally correct for phase aberration, but exhibit lower speed of sound contrast than the corresponding values provided by other sources, such as the ITIS database⁵⁴ (3515 m/s for cortical bone, 2118 m/s for trabecular). It is likely that increased reflection between bone layers provided by increased acoustic impedance contrast would increase the contribution of internally-reflected rays.

To test the hypothesis that low longitudinal speed of sound contrast between the skull layers could underestimate the role of internal reflection, simulations were performed in which the higher value for the longitudinal speed of sound in cortical bone of 3515 m/s was used in cortical bone, while the longitudinal speed of sound in the trabecular bone was varied between 2118 and 2500 m/s. The maximum amplitude and the

shift in maximum acoustic amplitude position are shown in Fig. 12. As expected, the difference in maximum amplitude caused by the inclusion of internally reflected waves is greater when the acoustic velocity contrast between cortical and trabecular bone is increased, but the overall difference between the two remains small. The decrease in amplitude at the focus compared to previous results may be due in part to the decreased longitudinal transmission at the cortical/water and cortical/trabecular interfaces caused by greater acoustic impedance contrast between these layers. The focal shift plotted as a function of trabecular fraction shown in Fig. 12(b) shows little focal position dependence regardless of the acoustic impedance between cortical and trabecular layers.

There may also be an effect of ray tracing-based simulation of a tri-layered skull model to overestimate the amount of reflection, as it is computed with the assumption of a discrete boundary existing between cortical and trabecular regions of bone, while real bone may not exhibit such clearly demarcated boundaries. The addition of multiple layers may also increase the overall attenuation of the skull in heterogeneous models compared to homogeneous models, even when the average absorption coefficient remains the same.

Previous examinations of the influence of skull heterogeneity on transcranial ultrasonic focusing have found that multi-layer reflection and resonance can be neglected, since the energy of reflected waves is negligible compared with that of the primary wavefront as a result of attenuation in bone and energy loss due to reflection at the tissue-bone boundary.¹⁰ This finding is supported by the investigation of the effects of internal reflection presented in this study.

# C. Experimental comparison

## 1. Acoustics

Comparison of experimentally measured and simulated acoustic field presented in Sec. III C 1 showed good agreement in both the magnitude and spatial distribution of the pressure field. When focusing through the skull, the

<figure><img src="image_25.png" /></figure>

<figure><img src="image_26.png" /><figcaption>FIG. 12. (Color online) (a) Maximum acoustic amplitude at the focus as a function of the longitudinal speed of sound in the trabecular bone layer. (b) Change in the position at which the maximum amplitude occurs induced by the inclusion of internal reflection at skull/water and cortical/trabecular boundaries, plotted as a function of the longitudinal speed of sound in the trabecular bone layer.</figcaption></figure>
simulated pressure field showed considerably better spatial confinement than the experimentally measured acoustic field—an expected outcome due to the simplicity of the homogeneous skull model representation.

## 2. Temperature modeling

Good agreement is shown between experimental and simulated temperature increase presented in Sec. III C 2. Consistent with the comparison of the experimental and simulated acoustic field, the simulated temperature distribution showed better confinement of the temperature increase compared to experimental measurements, even after spatial averaging of the simulated temperature field was applied.

The rate of change of the maximum temperature differed by 1.68% between experiments and simulations after spatial averaging was applied to provide comparable spatial resolutions for temperature measurements. However, the experimental data showed an estimated temperature rise of 1.865 °C for zero acoustic power. This may be partially explained by noise during MR thermometry, which would also account for spatial fluctuations when the temperature is plotted as a function of lateral position for 100 W acoustic power in Fig. 11(a).

## V. CONCLUSION

The result of the homogeneous parametric study quantifies the sensitivity of focal intensity and position to uncertainty in acoustic properties of the human skull used for numerical simulations for transcranial applications. These findings illustrate the importance of accurately determining the relationship between density and acoustic velocity in ray-tracing models of transcranial ultrasound.

We demonstrated that using a ray tracing method, the predictions of peak pressure and the location of the focal spot are not highly sensitive to uncertainties and disagreement in the precise relationship between density, longitudinal speed of sound, and attenuation. Prediction of amplitude, and by extension bioeffects such as temperature rise, may require accurate characterization of tissue properties.

Overall, the amount of trabecular bone in the skull relative to the overall skull volume has been shown to decrease the maximum acoustic amplitude produced at the focus, which is related to the amount of heterogeneity within the skull. The findings of this work illustrate how increased skull heterogeneity can have a deleterious effect on the ability to effectively transmit energy through the skull. It may be possible that a homogeneous model of the skull may be sufficient for use within simulations, provided the material properties of the model are generated by accounting for the overall properties of the physical skull.

Internal reflection between skull layers was shown to have a relatively small impact on transcranial focal quality and shows a smaller significance on the acoustic field than the effects of mode conversion. Results indicate that accurate estimation of acoustic parameters using a ray-tracing model may require consideration of the heterogeneity of the skull.

Good agreement was observed between experimental and simulated temperature rise when a numerical temperature model was combined with the ray-tracing acoustic simulations.

## ACKNOWLEDGMENTS

This work was supported by the FUS Foundation and benefited from the technical support of the company InSightec.

<footnote n="1">¹K. Hynynen and G. Clement, “Clinical applications of focused ultrasound—The brain,” *Int. J. Hyperthermia* **23**(2), 193–202 (2007).</footnote>

<footnote n="2">²J. G. Lynn, R. L. Zwemer, A. J. Chick, and A. E. Miller, “A new method for the generation and use of focused ultrasound in experimental biology,” *J. General Physiol.* **26**(2), 179–193 (1942).</footnote>

<footnote n="3">³W. Qiu, A. Bouakaz, E. E. Konofagou, and H. Zheng, “Ultrasound for the Brain: A Review of Physical and Engineering Principles, and Clinical Applications,” *IEEE Trans. Ultrason. Ferroelectr. Freq. Control* **68**(1), 6–20 (2021).</footnote>

<footnote n="4">⁴D. Jeanmonod, B. Werner, A. Morel, L. Michels, E. Zadicario, G. Schiff, and E. Martin, “Transcranial magnetic resonance imaging-guided focused ultrasound: Noninvasive central lateral thalamotomy for chronic neuropathic pain,” *Neurosurg. Focus* **32**(1), E1–E11 (2012).</footnote>

<footnote n="5">⁵W. J. Elias, D. Huss, T. Voss, J. Loomba, M. Khaled, E. Zadicario, R. C. Frysinger, S. A. Sperling, S. Wylie, S. J. Monteith, J. Druzgal, B. B. Shah, M. Harrison, and M. Wintermark, “A pilot study of focused ultrasound thalamotomy for essential tremor,” *New England J. Med.* **369**(7), 640–648 (2013).</footnote>

<footnote n="6">⁶A. Kyriakou, E. Neufeld, B. Werner, M. M. Paulides, G. Szekely, and N. Kuster, “A review of numerical and experimental compensation techniques for skull-induced phase aberrations in transcranial focused ultrasound,” *Int. J. Hyperthermia* **30**(1), 36–46 (2014).</footnote>

<footnote n="7">⁷A. Kyriakou, E. Neufeld, B. Werner, G. Székely, and N. Kuster, “Full-wave acoustic and thermal modeling of transcranial ultrasound propagation and investigation of skull-induced aberration correction techniques: A feasibility study,” *J. Therapeutic Ultrasound* **3**(1), 11 (2015).</footnote>

<footnote n="8">⁸G. T. Clement, P. J. White, R. L. King, N. McDannold, and K. Hynynen, “A magnetic resonance imaging-compatible, large-scale array for transskull ultrasound surgery and therapy,” *J. Ultrasound Med.* **24**(8), 1117–1125 (2005).</footnote>

<footnote n="9">⁹J.-F. Aubry, M. Tanter, M. Pernot, J.-L. Thomas, and M. Fink, “Experimental demonstration of noninvasive transskull adaptive focusing based on prior computed tomography scans,” *J. Acoust. Soc. Am.* **113**(1), 84–93 (2003).</footnote>

<footnote n="10">¹⁰C. Jiang, D. Li, F. Xu, Y. Li, C. Liu, and D. Ta, “Numerical evaluation of the influence of skull heterogeneity on transcranial ultrasonic focusing,” *Front. Neurosci.* **14**, 1–12 (2020).</footnote>

<footnote n="11">¹¹G. T. Clement and K. Hynynen, “A non-invasive method for focusing ultrasound through the human skull,” *Phys. Med. Biol.* **47**(8), 1219–1236 (2002).</footnote>

<footnote n="12">¹²G. T. Clement and K. Hynynen, “Correlation of ultrasound phase with physical skull properties,” *Ultrasound Med. Biol.* **28**(5), 617–624 (2002).</footnote>

<footnote n="13">¹³C. W. Connor, G. T. Clement, and K. Hynynen, “A unified model for the speed of sound in cranial bone based on genetic algorithm optimization,” *Phys. Med. Biol.* **47**(22), 3925–3944 (2002).</footnote>

<footnote n="14">¹⁴S. Pichardo, C. Moreno-Hernández, R. A. Drainville, V. Sin, L. Curiel, and K. Hynynen, “A viscoelastic model for the prediction of transcranial ultrasound propagation: Application for the estimation of shear acoustic properties in the human skull,” *Phys. Med. Biol.* **62**(17), 6938–6962 (2017).</footnote>

<footnote n="15">¹⁵R. A. Drainville, L. Curiel, and S. Pichardo, “Superposition method for modelling boundaries between media in viscoelastic finite difference time domain simulations,” *J. Acoust. Soc. Am.* **146**(6), 4382–4401 (2019).</footnote>

<footnote n="16">¹⁶B. E. Treeby, J. Jaros, D. Rohrbach, and B. T. Cox, “Modelling elastic wave propagation using the k-Wave MATLAB Toolbox,” in *Proceedings of the 2014 IEEE International Ultrasonics Symposium*, Chicago, IL (September 3–6, 2014), pp. 146–149.</footnote>

<footnote n="17">¹⁷G. W. Miller, M. Eames, J. Snell, and J.-F. Aubry, “Ultrashort echo-time MRI versus CT for skull aberration correction in MR-guided transcranial</footnote>
focused ultrasound: In vitro comparison on human calvaria," *Medical Phys.* **42**(5), 2223-2233 (2015).

¹⁸F. J. Fry and J. E. Barger, "Acoustical properties of the human skull," *J. Acoust. Soc. Am.* **63**(5), 1576-1590 (1978).

¹⁹M. Pérez-Liva, J. L. Herraiz, J. M. Udías, E. Miller, B. T. Cox, and B. E. Treeby, "Time domain reconstruction of sound speed and attenuation in ultrasound computed tomography using full wave inversion," *J. Acoust. Soc. Am.* **141**(3), 1595-1604 (2017).

²⁰T. Bancel, A. Houdouin, P. Annic, I. Rachmilevitch, Y. Shapira, M. Tanter, and J.-F. Aubry, "Comparison between ray-tracing and full-wave simulation for transcranial ultrasound focusing on a clinical system using the transfer matrix formalism," *IEEE Trans. Ultrason. Ferroelectr. Freq. Control* **68**(7), 2554-2565 (2021).

²¹J. Robertson, E. Martin, B. Cox, and B. E. Treeby, "Sensitivity of simulated transcranial ultrasound fields to acoustic medium property maps," *Phys. Med. Biol.* **62**(7), 2559-2580 (2017).

²²Q. Grimal, D. Rohrbach, J. Grondin, R. Barkmann, C. C. Glüer, K. Raum, and P. Laugier, "Modeling of femoral neck cortical bone for the numerical simulation of ultrasound propagation," *Ultrasound Med. Biol.* **40**(5), 1015-1026 (2014).

²³S. A. Leung, T. D. Webb, R. R. Bitton, P. Ghanouni, and K. Butts Pauly, "A rapid beam simulation framework for transcranial focused ultrasound," *Sci. Rep.* **9**(1), 7965 (2019).

²⁴T. D. Webb, S. A. Leung, J. Rosenberg, P. Ghanouni, J. J. Dahl, N. J. Pelc, and K. B. Pauly, "Measurements of the Relationship Between CT Hounsfield Units and Acoustic Velocity and How It Changes With Photon Energy and Reconstruction Method," *IEEE Trans. Ultrasonics, Ferroelectrics, Frequency Control* **65**(7), 1111-1124 (2018).

²⁵U. Vyas and D. Christensen, "Ultrasound beam simulations in inhomogeneous tissue geometries using the hybrid angular spectrum method," *IEEE Trans. Ultrasonics Ferroelectr. Freq. Control* **59**(6), 1093-1100 (2012).

²⁶S. Pichardo, V. W. Sin, and K. Hynynen, "Multi-frequency characterization of the speed of sound and attenuation coefficient for longitudinal transmission of freshly excised human skulls," *Phys. Med. Biol.* **56**(1), 219-250 (2011).

²⁷A. Pulkkinen, B. Werner, E. Martin, and K. Hynynen, "Numerical simulations of clinical focused ultrasound functional neurosurgery," *Phys. Med. Biol.* **59**(7), 1679-1700 (2014).

²⁸G. Pinton, J.-F. Aubry, E. Bossy, M. Muller, M. Pernot, and M. Tanter, "Attenuation, scattering, and absorption of ultrasound in the skull bone," *Med. Phys.* **39**(1), 299-307 (2011).

²⁹T. D. Mast, "Empirical relationships between acoustic parameters in human soft tissues," *Acoust. Res. Lett. Online* **1**(2), 37-42 (2000).

³⁰P. J. White, G. T. Clement, and K. Hynynen, "Longitudinal and shear mode ultrasound propagation in human skull bone," *Ultrasound Med. Biol.* **32**(7), 1085-1096 (2006).

³¹S. Chatillon, R. Loyet, L. Brunel, F. Chavrier, N. Guillen, and S. Le Berre, "Applications of intensive HIFU simulation based on surrogate models using the CIVA HealthCare platform," *J. Phys. Conf. Ser.* **1761**(1), 012007 (2021).

³²N. Gengembre and A. Lhémery, "Pencil method in elastodynamics: Application to ultrasonic field computation," *Ultrasonics* **38**(1), 495-499 (2000).

³³A. Lhémery, P. Calmon, I. Lecoeur-Taïbi, R. Raillon, and L. Paradis, "Modeling tools for ultrasonic inspection of welds," *NDT E Int.* **33**(7), 499-513 (2000).

³⁴A. Gardahaut, K. Jezzine, D. Cassereau, and N. Leymarie, "Simulation of ultrasonic wave propagation in welds using ray-based methods," *J. Phys. Conf. Ser.* **498**(1), 012008 (2014).

³⁵J.-F. de Belleval, N. Gengembre, N. Lhémery, and P. Calmon, "Modeling of ultrasonic beams," in *Materials and Acoustics Handbook*, edited by M. Bruneau and C. Potel (ISTE & Wiley, London, UK, 2009), pp. 359-398.

³⁶V. Cerveny, *Seismic Ray Theory*, 1st ed. (Cambridge University Press, Cambridge, UK, 2001).

³⁷J.-F. Aubry, O. Bates, C. Boehm, K. Butts Pauly, D. Christensen, C. Cueto, P. Gélat, L. Guasch, J. Jaros, Y. Jing, R. Jones, N. Li, P. Marty, H. Montanaro, E. Neufeld, S. Pichardo, G. Pinton, A. Pulkkinen, A. Stanziola, A. Thielscher, B. Treeby, and E. van't Wout, "Benchmark problems for transcranial ultrasound simulation: Intercomparison of compressional wave models," *J. Acoust. Soc. Am.* **152**(2), 1003-1019 (2022).

³⁸W. J. Elias, N. Lipsman, W. G. Ondo, P. Ghanouni, Y. G. Kim, W. Lee, M. Schwartz, K. Hynynen, A. M. Lozano, B. B. Shah, D. Huss, R. F. Dallapiazza, R. Gwinn, J. Witt, S. Ro, H. M. Eisenberg, P. S. Fishman, D. Gandhi, C. H. Halpern, R. Chuang, K. Butts Pauly, T. S. Tierney, M. T. Hayes, G. R. Cosgrove, T. Yamaguchi, K. Abe, T. Taira, and J. W. Chang, "A randomized trial of focused ultrasound thalamotomy for essential tremor," *New England J. Med.* **375**(8), 730-739 (2016).

³⁹A. Jameel, P. Bain, D. Nandi, B. Jones, and W. Gedroyc, "Device profile of exAblate Neuro 4000, the leading system for brain magnetic resonance guided focused ultrasound technology: An overview of its safety and efficacy in the treatment of medically refractory essential tremor," *Expert Rev. Med. Dev.* **18**(5), 429-437 (2021).

⁴⁰A. Franzini, S. Moosa, F. Prada, and W. J. Elias, "Ultrasound ablation in neurosurgery: Current clinical applications and future perspectives," *Neurosurgery* **87**(1), 1-10 (2020).

⁴¹K. W. Chang, Y. S. Park, and J. W. Chang, "Skull factors affecting outcomes of magnetic resonance-guided focused ultrasound for patients with essential tremor," *Yonsei Med. J.* **60**(8), 768-773 (2019).

⁴²P. B. Rosnitskiy, P. V. Yuldashev, O. A. Sapozhnikov, L. R. Gavrilov, and V. A. Khokhlova, "Simulation of nonlinear trans-skull focusing and formation of shocks in brain using a fully populated ultrasound array with aberration correction," *J. Acoust. Soc. Am.* **146**(3), 1786-1798 (2019).

⁴³B. E. Treeby and B. T. Cox, "Modeling power law absorption and dispersion for acoustic propagation using the fractional Laplacian," *J. Acoust. Soc. Am.* **127**(5), 2741-2748 (2010).

⁴⁴C. Lafon, D. Moore, M. D. C. Eames, J. Snell, R. A. Drainville, F. Padilla, A. Drainville, and F. Padilla, "Evaluation of pseudorandom sonifications for reducing cavitation with a clinical neurosurgery HIFU device," *IEEE Trans. Ultrason. Ferroelectr. Freq. Control* **68**(4), 1224-1233 (2021).

⁴⁵H. H. Pennes, "Analysis of tissue and arterial blood temperatures in the resting human forearm," *J. Appl. Physiol.* **1**(2), 93-122 (1948).

⁴⁶C. P. Phenix, M. Togtema, S. Pichardo, I. Zehbe, and L. Curiel, "High intensity focused ultrasound technology, its scope and applications in therapy and drug delivery," *J. Pharmacy Pharm. Sci.* **17**(1), 136-153 (2014).

⁴⁷A. Lakhssassi, E. Kengne, and H. Semmaoui, "Modified pennes' equation modelling bio-heat transfer in living tissues: Analytical and numerical analysis," *Nat. Sci.* **02**(12), 1375-1385 (2010).

⁴⁸I. M. Hallaj, R. O. Cleveland, and K. Hynynen, "Simulations of the thermo-acoustic lens effect during focused ultrasound surgery," *J. Acoust. Soc. Am.* **109**(5), 2245-2253 (2001).

⁴⁹R. van Vossen, J. O. a Robertsson, and C. H. Chapman, "Finite-difference modeling of wave propagation in a fluid-solid configuration," *Geophysics* **67**(2), 618-624 (2002).

⁵⁰R. A. Drainville, "Principles of superposition in viscoelastic Finite-Difference modelling," Ph.D. thesis, Lakehead University, Thunder Bay, Canada (2019).

⁵¹V. Rieke and K. Butts Pauly, "MR thermometry," *J. Magn. Resonance Imag.* **27**(2), 376-390 (2008).

⁵²N. McDannold, P. J. White, and R. Cosgrove, "Elementwise approach for simulating transcranial MRI-guided focused ultrasound thermal ablation," *Phys. Rev. Res.* **1**(3), 033205 (2019).

⁵³G. Clement, J. Sun, and K. Hynynen, "The role of internal reflection in transskull phase distortion," *Ultrasonics* **39**(2), 109-113 (2001).

⁵⁴P. Hasgall, E. Neufeld, M. C. Gosselin, A. Klingenböck, and N. Kuster, "IT'IS Database for thermal and electromagnetic parameters of biological tissues, Version 4.1," *itis.swiss/database* (Last viewed March 14, 2023).