# Transcranial ultrasound simulations: A review

Célestine Angla, Benoit Larrat, Jean-luc Gennisson, Sylvain Chatillon

► To cite this version:

Célestine Angla, Benoit Larrat, Jean-luc Gennisson, Sylvain Chatillon. Transcranial ultrasound sim-
ulations: A review. Medical Physics, 2022, 1, pp.22. 10.1002/mp.15955 . hal-03857081

HAL Id: hal-03857081

https://hal.science/hal-03857081

Submitted on 17 Nov 2022

**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>
TRACK & TREND
YOUR QA DATA WITH

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

RITtrend™

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

RITtrend™ is the all-in-one SQL statistical database solution that tracks and trends data for all of your department's QA measurements.

RITtrend allows you to view your current and historical data. Then, based on any identified trends, you can establish well-defined, well-reasoned tolerance specifications that are particular to your machine. Within the RITtrend interface, you can set minimum and maximum values, view moving means, identify outliers, and monitor trends or drift.

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

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

Set baseline measurements for statistical process control to track historical data over time.

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

Easily import your analysis results into the RITtrend database with a single click.

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

Set fully-customizable parameters to track and trend any pass/fail data in RIT's Tolerance Manager.

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

Enable email alerts to be instantly alerted of any data discrepancies in real-time.

RITtrend is a built-in feature with the following RIT Family of Products software packages:

RIT Complete

RITG142

Radia

RIT Classic

RITG148†

RITG135

CLICK TO VISIT RADIMAGE.COM TODAY
TO DEMO THIS SOFTWARE FEATURE
FROM RADIOLOGICAL IMAGING
TECHNOLOGY, INC.

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

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

+1 (719) 590-1077, OPT. 4

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

SALES@RADIMAGE.COM

Connect with us on social media

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

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

<figure><img src="image_14.png" /></figure>
REVIEW ARTICLE

MEDICAL PHYSICS

# Transcranial ultrasound simulations: A review

Célestine Angla<fnref n="1" />

Benoit Larrat<fnref n="2" />

Jean-Luc Gennisson<fnref n="3" />

Sylvain Chatillon<fnref n="1" />

**Correspondence**

Célestine Angla, Université Paris-Saclay, CEA, List, F-911120 Palaiseau, France.

Email: celestine.angla@cea.fr

**Funding information**

Commissariat à l'Énergie Atomique et aux Énergies Alternatives (CEA)

**Abstract**

Transcranial ultrasound is more and more used for therapy and imaging of the brain. However, the skull is a highly attenuating and aberrating medium, with different structures and acoustic properties among samples and even within a sample. Thus, case-specific simulations are needed to perform transcranial focused ultrasound interventions safely. In this article, we provide a review of the different methods used to model the skull and to simulate ultrasound propagation through it.

**KEYWORDS**

CT-scan acoustical properties, skull acoustical properties, skull modeling, transcranial acoustic simulations, transcranial ultrasound focusing, transcranial ultrasound simulations

## 1 | INTRODUCTION

Transcranial ultrasound is more and more used for therapy and imaging of the brain. Therapeutic applications include thermal ablation,<fnref n="1" /> blood-brain-barrier (BBB) opening,<fnref n="2" /> and neuro-modulation.<fnref n="3" /> Transcranial ultrasound imaging has been used for a long time for vascular and flow imaging through the temporal windows, but its generalization is still impeded by the strong image distortion and sensitivity loss at typical imaging frequencies.⁴

Making ultrasound cross the skull is of high interest, but it is also a big challenge as the skull is highly attenuating and aberrating. In addition, the pores located in the trabecular bone act as scatterers for ultrasound. As a result, ultrasound beams are both spatially shifted and attenuated. As the structure and acoustic properties of the skull differ between samples and even within a sample, case-specific simulations are needed to perform transcranial focused ultrasound interventions safely. The first issue when performing transcranial acoustic simulations is skull modeling. Due to the big variations within a skull and among specimens, most studies use images of the skull to extract its geometric description and derive acoustic properties using empirical relationships. However, these relationships are only hypothetical and the imaging resolution is often not precise enough to fully describe the porous structure of the skull. Once the

skull is modeled, several simulation methods can be used, which can be divided into two main categories: numerical methods and semi-analytical ones. The type of method used usually depends on the application and the aim is to find the best balance between accuracy and computation performances. Likewise, approximations, such as neglecting shear waves and nonlinear effects, can be made to speed up the simulations but the loss in accuracy has to be evaluated. After a brief overview of a few applications of transcranial acoustic simulations in the first part, parts one and two will focus on the different skull modeling methods that were previously introduced, and on the simulation methods used for transcranial ultrasound. Finally, the methods used to compare simulations and experiments will be described in the last part. Indeed, another main issue with transcranial acoustic simulations is to evaluate their accuracy. The only way to do so is to compare the results of the simulations with the results of experiments in a setup as close as possible to the simulated one.

## 2 | APPLICATIONS OF TRANSCRANIAL SIMULATIONS

Transcranial acoustic simulations have many applications. Their main goal is to optimize the ultrasound treatment, by predicting the acoustic dose, the

This is an open access article under the terms of the Creative Commons Attribution License, which permits use, distribution and reproduction in any medium, provided the original work is properly cited.

© 2022 The Authors. Medical Physics published by Wiley Periodicals LLC on behalf of American Association of Physicists in Medicine.

<footnote n="1">¹Université Paris-Saclay, CEA, List, Palaiseau, France</footnote>

<footnote n="2">²Université Paris Saclay, CNRS, CEA, DRF/JOLIOT/NEUROSPIN/BAOBAB, Gif-sur-Yvette, France</footnote>

<footnote n="3">³Université Paris Saclay, CNRS, Inserm, BioMaps, Orsay, France</footnote>
focalization characteristics, and for operation planning. When combined with a multielement transducer or an acoustic lens,⁵ simulations can be used to correct the phase and amplitude (especially at high frequencies where the distortion induced by the skull is higher) of the emitted signal in order to achieve good focusing through the skull or for imaging applications.⁶⁻⁸ For thermal therapies, such as subthalamic nuclei ablation for the treatment of essential tremor,⁹ heat simulations are also often performed in order to predict whether heating occurs at the planned location and to avoid burning surrounding tissues. Acoustic simulations are thus needed in this case, as the pressure field obtained is the input of the thermal simulations. Thermal simulations are also used for nonthermal applications, in order to predict that no unwanted heating happens. Generally speaking, in addition to optimizing the treatment, acoustic simulations are often performed to ensure the safety of ultrasound therapies. Apart from heating prediction, they are used for standing waves and cavitation prediction.¹⁰⁻¹⁵ Finally, acoustic simulations can be used to optimize the positioning of the transducer.¹⁶

## 2.1 | Phase and amplitude correction

Transcranial simulations can be used to correct the phase difference induced by the skull in order to obtain nearly optimal focusing. To do that, several strategies can be used: classical time reversal method and cross-correlation.

Time reversal consists in simulating wave propagation from a virtual source point located inside the brain at the intended focusing location and recording the signals on each element of the transducer. Then, if only time derivatives of even order are considered in the simulation equation, one can achieve focusing on the virtual source point by emitting the reversed signals received by each element of the transducer.

The cross-correlation method aims to facilitate the experimental backpropagation, by using only the phase delays (and not the whole reversed signals). It is based on the following equation¹⁷:

$$p_{\text{reference}} \cdot p_n(\tau) = \int p_{\text{reference}}(t) p_n(t + \tau) dt \quad (1)$$

where $\tau$ is the phase delay.

As the skull is a strongly absorbing medium with varying thickness across its surface, one can find a way to distribute the emitted energy across the surface, so the transducer elements near the thickest parts, and thus the most absorbing parts, of the skull, emit less power than the elements in front of the thinnest locations. Doing so will ensure a better transmission of the waves while avoiding unnecessary heating of the skull in the most absorbing locations. This distribution of energy can be achieved with amplitude correction.¹⁵,¹⁷⁻²⁰

Acoustic simulations can also be combined with the emission of a cavitation bubble that reflects the ultrasound waves.²¹ This method uses acoustic simulations as a first step to focus near the intended location that enables the creation of a cavitation bubble, whose signal is recorded by the transducer array. The bubble is then used an invasive source inside the brain to obtain accurate focusing. For safety reasons, the relevance of this approach for clinical applications is debatable.

## 2.2 | Heating prediction

The pressure field computed by the transcranial acoustic simulation can be entered as input to a heat transfer equation to simulate the heating of the skull and brain, and avoid burning tissues.¹⁹,²⁰,²²⁻²⁸

This heat transfer equation is called the Pennes bioheat equation and is given by

$$\rho C \frac{\partial T}{\partial t} = \nabla \cdot (k \nabla T) - \rho_b w_b C_b (T - T_b) + \frac{\alpha p^2}{\rho c} \quad (2)$$

where $\rho$ is the density of the tissue, $C$ its specific heat capacity, $k$ its thermal conductivity, $\alpha$ its absorption coefficient, and $c$ its speed of sound. $\rho_b$ is the blood density, $C_b$ its specific heat capacity, $T_b$ its temperature, and $\omega_b$ the perfusion rate. $p$ is ultrasound pressure derived from the acoustic model.

The first term of the right side models the thermal diffusion, the second one models the effects of perfusion, and the last one models the acoustical power (found with the acoustic simulation). In Kyriakou et al.²⁰ and Guo et al.,²⁶ add another term to the equation to model the body heat production. The bioheat is often solved using 3D finite differences.

The criteria for a safe transcranial treatment are that the temperature at the focal region is higher than 60°C and that there are no regions under the skull whose thermal dose is higher than 90 equivalent minutes at 43°C.¹⁹

Using the temperature found after solving the bioheat equation, the thermal dose can be computed as

$$TD = \int_{0}^{t} R(T(t'))^{T_{\text{ref}}-T(t')} dt' \quad (3)$$

With $T_{\text{ref}} = 43^{\circ}\text{C}$ and $R(T) = \begin{cases} 0.50 & \text{for } \tau \ge T_{\text{ref}} \\ 0.25 & T < T_{\text{ref}} \end{cases}$.

The thermal dose is defined as cumulative equivalent minutes at 43°C. For thermal ablation, the thermal dose can be used to compute the power required to achieve a thermal dose high enough to cause irreversible damage in the brain tissues. For instance, Pulkkinen et al.²⁷ performed the sonifications with the power required to achieve an equivalent thermal dose
of 25 min at the focus. In addition to ensuring that the target location is heated enough to perform ablation, the thermal dose can also be used to check whether heating occurs at undesired locations. For example, in Pulkkinen et al.,²⁷ target positions inducing a thermal dose of more than 5 min, at a location adjacent to the skull, were considered to be untreatable.

For nonthermal applications such as BBB opening, heat simulations can be performed to ensure that no undesired heating occurs during the treatment. For instance, Marquet et al.²² estimated the maximum heat dissipation increase in the skull and find that after 60 s of sonication, the temperature increase is lower than 0.03°C, which is negligible. In fact, this is due to the very short duty cycles, which let the skull bone cool enough between each pulse.

During most transcranial treatments, active cooling of the skull is performed in order to limit skull heating. Pulkkinen et al.²⁷ simulated the circulating water at 15°C for 15 min prior to sonifications using the Dirichlet boundary condition with constant temperature of 15°C around the skull. The simulated field is then used as the input temperature field for the sonifications simulations. Kyriakou et al.²⁰ used a similar method to simulate the cooling procedure. Pulkkinen et al.²⁷ also simulated the active cooling of the skull base. This cooling method is based on circulating cooled water prior to sonifications in a large nasal cavity at the center of the skull base. It is very important to apply such a method when targets are located near the skull base, as measurements have shown that the temperature in the soft tissue adjacent to the bone can exceed that of the focus.²⁷ In this study,²⁷ the nasal cavity cooling was found to increase the treatment envelope for the uncorrected sonifications, whereas for the phase corrected sonifications, the nasal cavity cooling was found to sometimes have negative effect on the sonifications if the focus was to close the cooled area (in those cases, the temperature also decreased at the focus).

Most thermal simulations are performed assuming that thermal tissue properties are not temperature dependent. However, the high temperature used for thermal ablation can cause vascular shutdown, which prevents the perfusion in those locations and cause the temperature to increase faster. Kyriakou et al.²⁰ simulated this effect by assuming that perfusion decreases linearly when the tissue temperature is above 50°C and stops above 51°C. The results show that after 20 s of sonication, only 1°C of the total additional temperature increase was observed in tissues where vascular shutdown occurred, suggesting that this effect can be neglected.

Heat simulations can also be used to predict thermal lesions. McDannold et al.²⁴ gathered the posttreatment images from 40 clinical TcMRgFUS treatments, including 16 for which bone marrow lesions were observed. They predicted the presence of lesions using a thresh-

old for the acoustic energy of 18.1–21.1 kJ (maximum acoustic energy used) and 97.0–112.0 kJ (total acoustic energy applied over the whole treatment). The results show that the size of the lesions was not always predicted by the acoustic energy. However, the locations, sizes, and shapes of the heated regions estimated by the thermal simulations were qualitatively similar to those of the lesions and the lesions generally appeared in areas with high predicted temperatures.

## 2.3 | Transcranial FUS safety assessment

Ultrasound can produce cavitation phenomenon and thus create bubbles that will oscillate with the varying pressure field. If the bubble oscillates too much and collapses, the pressure and temperature can increase very fast, and blood vessels can be damaged, causing hemorrhages. Simulations can be used to predict cavitation phenomenon¹³⁻¹⁵ and avoid hemorrhages in the brain. Several studies have used simulation to try to understand phenomenon that occurred during clinical trials.¹³,²⁹ It can help avoid these effects in the future. However, simulation can also be used to prepare case-specific treatments and prevent those undesirable effects.

The standard indicator to evaluate the likelihood of cavitation risks is the mechanical index (MI), defined by

$$MI = \frac{P_{-}}{\sqrt{f}} \qquad (4)$$

where $P_{-}$ the peak rarefactional pressure and $f$ the frequency. This is only an indicator as cavitation depends on other factors that are currently unknown. To avoid adverse biological effects related to acoustic cavitation, the FDA (Food & Drug Administration—510K Norm 1992) imposes the diagnostic devices to ensure an MI less than 1.9.¹³

In particular, studies²²,⁵⁵ have focused on the TRUMBI (transcranial low-frequency ultrasound-mediated thrombolysis in brain ischemia) clinical trial that has been stopped prematurely because hemorrhages were observed on several patients. Baron et al.¹³ simulated the pressure fields used in the trial, and they computed the associated MI. They found that the peak negative pressure is higher than the cavitation threshold in large areas of the brain, which is mainly due to the presence of standing waves. Using acoustic simulations, Pinton et al.²⁹ found that a volume of 2.7 cm³ was above the MI threshold in the TRUMBI trial. They also investigated the risk of cavitation at 220 kHz and 1 MHz. They found that for an equivalent energy deposition rate and the same geometry, the brain volume above the MI ($0.4 \le MI \le 3$) is three to four orders of magnitude larger at 220 kHz than at
1 MHz. This high correlation between frequency and the volume with a high probability of cavitation is due to three reasons. First, the focal volume decreases with the frequency. Second, the heat deposition increases when frequency increases (as absorption increases). Third, the pressure needed to obtain a given MI decreases as the square root of the frequency. Both studies claim that unfocused transducers (as the one used in the TRUMBI trial) should not be used, in order to avoid creating hot spots at unwanted locations due to the focusing effect of the skull.

Top et al.¹⁵ used simulation to investigate BBB disruption that was observed in the pre-focal region during previous experiments at 220 kHz. Their results are in agreement with the disruption in the pre-focal region. However, they observed side lobes in the post-focal region in the simulations that did not seem to have had any effects in the experiments. They suggested this difference might be due to shielding of the ultrasound field due to microbubble activity in the focal region. They also simulated the field produced by the passive cavitation detector at its resonant frequency (610 kHz) and at the subharmonic of the transducer (110 kHz), as the sensitivity pattern is proportional to the transmitted field.

Transcranial acoustic simulations can also help predicting standing waves, which are contributors to undesired cavitation effects. Deffieux and Konofagou¹⁰ and Mueller et al.¹¹ estimated standing waves using a filter that detects their characteristic pattern. First, the fast oscillations of the field are extracted along each dimension using a high-pass spatial filter. Second, a Hilbert transform is applied to obtain the envelope of this signal. Third, only the maximum of the three dimensions is kept so as to have a unique map. The resulting field estimates the maximum of the spatial modulation pattern that is often associated with the presence of standing waves. Zhang et al.¹² used what they call the standing-wave ratio to predict the intensity of standing waves. It is defined as

$$R_a = \frac{P_{max} - P_{min}}{P_{avg}} \qquad (5)$$

Methods for reducing standing waves have also been developed and verified thanks to simulations. Baron et al.¹³ suggested using higher frequencies (more absorption so less interferences), reducing the duration of the pulses (to avoid interferences), and applying modulation. Deffieux and Konofagou¹⁰ proposed to use fast periodic linear chirps to reduce standing waves. Chirps are signals with a varying frequency over time. It induces a time-dependent phase difference between the incident and reflected waves which ensures that, even if interferences occur, their position will change over time and thus, after some time, the effects of constructive and destructive interference will even out. When using

a chirp with a period of 23 µs and a frequency oscillating between 450 and 550 kHz, they manage to reduce the peak amplitude of standing waves (as a percentage of the peak pressure) from 19% to 12% and to reduce the standing wave volume (that is to say the volume with standing wave amplitude higher than 5% of the peak pressure) from 0.87% to 0.17%. Although the standing wave volume reduction is quite satisfying (5 times reduction), the standing wave amplitude reduction is not as significant (1.5 times reduction). The authors suggested it may be due to interferences occurring close to the skull surface, because the time delay between the incident and reflected waves is too small compared to the frequency shift. Zhang et al.¹² used the Randi function to generate a random phase between 0 and π, which is added to the original signal at regular intervals to break the standing wave formation condition.

## 3 | SKULL MODELING

Skull modeling is an important factor for accurate simulations. As the shape and acoustical properties of the skull vary between individuals, and even within one individual (across the years and among the whole surface of the skull), it is very tricky to create good skull models. Most skull models are based on skull imaging of the considered individual, as, for example, X-rays.

### 3.1 | Shape

When finite differences or other voxel-based simulation algorithms are used, the skull shape does not need to be extracted. However, for a geometry-based simulation algorithm, one needs to model the skull shape.

Hayner and Hynynen³⁰ modeled the skull as a medium contained between two flat, but not necessarily parallel, interfaces, as in the frequency range considered, the minimum radius of curvature of a skull is greater than the wavelength.

However, in most studies, the skull is segmented from either computed tomography (CT) scans or magnetic resonance (MR) images, in order to have a more accurate modeling of the geometry. To obtain the shape of the skull from a tensor of intensity voxels, two steps are usually needed: a segmentation step and a mesh-making step.

Segmentation of the bone voxels is often performed using a threshold method. Clement and Hynynen,³¹ Clement et al.,³² and Yin and Hynynen³³ identified the inner and outer surfaces of the skull by searching for the innermost and outermost densities greater than 1400 kg m⁻³ along each line of a CT slice. Pichardo and Hynynen³⁴ performed the segmentation with the FMRIB Software Library.³⁵
Mesh-making is often performed using an iso-surface algorithm.34,36–39 In some cases, a condition is imposed so that no surface area element is greater than $(\lambda/2)^2$, with $\lambda$, the acoustic wavelength in water.37,38

For skull segmentation from magnetic resonance imaging (MRI), two UTE (ultrashort echo time) images are generally combined (using simple arithmetic operations) to separate skull voxels from other tissues using various thresholds.26,40 In Miller et al.,40 after thresholding, a spatial connectivity requirement is applied to the pixels in the bone class, to eliminate isolated pixels or groups of pixels that are clearly not part of the skull.

## 3.2 | Acoustical properties

Several studies have measured skull acoustical properties (density, speed of sound, and attenuation).30,41,42 Even though properties of pure bone could be tabulated, most of the bone volume is heterogeneous and it is hard to predict the effective properties due to partial volume effect of skull images. However, most studies try to deduce the acoustical properties of a given skull from images, acquired with either CT or MR.

### 3.2.1 | Density

#### Linear relationship between density and Hounsfield units (HU)

In Rho et al.,43 ultrasonic measurements show a linear relationship between density and HU. HU are arbitrary units that depend on the CT scanner used. It is thus necessary to have a specific conversion from HU to any physical measure for each sample.

Connor et al.44 introduced a method to compute the skull density from HU. Assuming that the relation between HU and density is linear, and ensuring that a sample of both water and air (which have known densities) appears in the CT scan, one has

$$ \rho = \kappa_1 H + \kappa_0 \quad (6) $$

With $\kappa_1 = \frac{1}{H_{water} - H_{air}}$ and $\kappa_0 = \frac{-H_{air}}{H_{water} - H_{air}}$.

This method is used in Refs. [6, 7, 28, 36–38, 45–47].

To confirm the linear relationship obtained, they make the hypothesis that density and HU are related by a second-order polynomial. They use the relationship relating the total mass of the skull M with the volume v of a voxel and the sum of all densities across the skull: $M = v \sum_{skull} \rho$. Using this additional relationship, they are able to fit a second-order polynomial for density, but they find that the quadratic term is negligible in front of the others (by five orders of magnitude). Thus, their linear approximation seems correct.

#### Porosity-based and equivalent relationships

Aubry et al.18 proposed a different method based on the bone porosity. HU are defined by

$$ H = 1000 \frac{\mu_x - \mu_{water}}{\mu_{bone} - \mu_{water}} \quad (7) $$

where $\mu_x$, $\mu_{water}$, and $\mu_{bone}$ denote the photoelectric linear attenuation of the explored tissue, water, and bone. They proposed a linear relationship:

$$ \mu_x = \Psi \mu_{water} + (1 - \Psi) \mu_{bone} \quad (8) $$

with $\Psi$ the bone porosity. Thus they have

$$ \psi = 1 - \frac{H}{1000} \quad (9) $$

The density can be computed from the porosity with

$$ \rho = \psi \rho_{water} + (1 - \psi) \rho_{bone} \quad (10) $$

This method is used in Ref. [11–13, 15, 19–21, 48–51].

This relationship implies a linear relationship between the bone density and the HU, as it was assumed in Equation (6). However, instead of building the linear fit from $H_{water}$ and $H_{air}$, they build it from $\rho_{water}$ and $\rho_{bone}$. In the first case, $H_{water}$ and $H_{air}$ are derived from direct measurements on the CT scans. Thus, they are more reliable than the values of $\rho_{water}$ and $\rho_{bone}$, which are taken from measurements from the literature that were probably not conducted in the exact same conditions as the experiments conducted in the study. On the other hand, if the density is assumed to be linearly related to HU in the range corresponding to bone, it is probably not the case for all values of HU. Thus, building the linear fit from values in water and cortical bone (which are the lower and upper bounds for HU in the skull) seems more accurate than using values from air and water.

Guo et al.26 used the same equations as before, but they changed the definition of porosity:

$$ \psi = 1 - \frac{H}{\max(H)} \quad (11) $$

Indeed, with Aubry's definition of Hounsfield units, one has $0 \le H \le 1000$. In this study, they generalize Aubry's formulas for raw CT data, where the maximum value of H is not necessarily 1000.

Deffieux and Konofagou10 normalized the apparent CT density to have $0 \le \rho_{CT} \le 1$, and then the density is given by

$$ \rho = \rho_{min} + (\rho_{max} - \rho_{min}) \rho_{CT} \quad (12) $$
This is in fact equivalent to the previous equation based on porosity as $ρ_{CT} = \frac{H}{\max(H)} = 1 - \Psi$.

Marsac et al.⁵² introduced a new relationship for density:

$$ρ = ρ_{min} + (ρ_{max} - ρ_{min}) \frac{H - H_{min}}{H_{max} - H_{min}} \quad (13)$$

This equation is similar to the porosity-based Equation (10) and is equivalent to it when $H_{min} = 0$, Indeed, if $H_{min} = 0$, one has $1 - \Psi = \frac{H - H_{min}}{H_{max} - H_{min}}$. Thus this equation is the more general form of the porosity based equation. This equation is used in Refs. [5, 53].

### Bone fraction-based relationship

Vyas et al.⁵⁴ proposed a method based on the fraction of bone in the voxel, which is defined by

$$f = \frac{\frac{H}{1000} \times \frac{\mu}{\rho_{water\_eff}}}{\rho_{bone} \times \frac{\mu}{\rho_{bone\_eff}} - \frac{\mu}{\rho_{water\_eff}}} \quad (14)$$

With $\frac{\mu}{\rho_{eff}} = \frac{\Sigma S(E) \times \frac{\mu}{\rho}(E)_{water}}{\Sigma S(E)}$, where $S(E)$ is the beam spectrum as a function of energy. Then the density is defined by

$$ρ = f \times ρ_{bone} + (1 - f) ρ_{water} \quad (15)$$

This method is similar to the porosity-based Equation (10). However, bone fraction is computed in a more complicated way than porosity, and Equation (14) seems not to be homogeneous.

### 3.2.2 | Speed of sound

#### Linear relationships between speed of sound and Hounsfield units

Fry and Barger⁵⁵ and in Vyas et al.⁵⁴ used an empirical linear relationship linking directly the speed of sound to Hounsfield units:

$$c = 1460 + 0.7096 H \quad (16)$$

Clement and Hynynen³¹ used a relationship between the average speed of sound in a skull and the average density of that same skull, obtained experimentally from 1000 measurements with 10 skulls at 0.51 MHz:

$$c = 2.06 \rho - 1540 \quad (17)$$

where c in m s⁻¹, for densities ρ between 1820 and 2450 kg m⁻³. Thus, c ∈ [2209; 3507] m s⁻¹. This

roughly corresponds to the range of values found in the literature.

As a linear relationship is assumed between density and HU, this relationship is similar to Equation (16) as it implies a linear relationship between speed of sound and HU.

#### Porosity-based and equivalent relationships

Similarly to the porosity-based relationship for density, Aubry et al.¹⁸ proposed a sound speed relationship based on porosity:

$$c = c_{min} + (1 - \psi) (c_{max} - c_{min}) \quad (18)$$

The linear relationship between velocity and porosity is justified by Carter and Hayes who showed that the elastic modulus of bone is proportional to the apparent density cubed. This method is used in Refs. [11–13, 15, 19–21, 26, 48–51]. This relationship is similar to Equations (16) and (17) as it implies a linear fit between speed of sound and HU. However, Equation (18) is more general as porosity does not depend on the CT scanner, whereas HU do.

In the same way as for the density equations, other studies have improved the porosity-based equation for the speed of sound so as to make it more general. In order to account for CT scans where HU are not normalized, Deffieux and Konofagou¹⁰ proposed

$$c = c_{min} + (c_{max} - c_{min}) ρ_{CT} \quad (19)$$

which is equivalent to the equation introduced by Marsac et al.⁵²:

$$c = c_{min} + (c_{max} - c_{min}) \frac{\rho - \rho_{min}}{\rho_{max} - \rho_{min}} \quad (20)$$

These equations are also used in Refs. [5, 53].

Marsac et al.⁵² searched for the best value of $c_{max}$ such that the simulations fit the experiments.

#### Relationships obtained with genetic algorithms

Connor et al.⁴⁴ used a genetic algorithm to find the optimal relationship between density and speed of sound at 0.74 MHz, using a success function that compares the phase difference between simulation and experiment.

They compared their result function with what could be expected for a cellular solid, as they approximated the skull bone as one. The speed of sound in a solid is given by $c = \sqrt{\frac{E}{\rho}}$ where E is the Young modulus of the material and ρ its density. For open-pore cellular materials, it can be shown that $E \propto \rho^2$, whereas for closed-pore materials, $E \propto \rho^3$. The curve of speed of sound against bone density obtained in this study agrees with these models. Indeed, it initially has a form similar to that of a square root function (as it would be for
an open-celled porous solid), then the model transitions into a linear function (as it would be for a closed-celled porous solid). In addition, the transitional density region is located between the density of trabecular bone and that of cortical bone.

Pichardo et al.³⁶ used a similar method, but they investigated the sound speed relationships at various frequencies as they had found from measurements that the skull is a dispersive medium. They performed a two-step optimization (to speed up the process while maintaining a good simulation accuracy) at several frequencies: 0.27, 0.836, 1.402, 1.965, and 2.525 MHz. The sound speed functions found in this study were used (either directly or interpolated to fit the needed frequency) in Refs. [6, 7, 28, 37, 38, 45, 46].

### Polynomial relationship

McDannold et al.²³,²⁴ assumed that the relationship between the inverse of the speed of sound and the skull density can be approximated by a series of polynomials:

$$ \frac{1}{c} = \sum_{m} B_{m} \rho^{m} \qquad (21) $$

The advantage of this formulation is that they can estimate the phase shifts resulting from a change in speed of sound with

$$ \begin{aligned} \Delta\phi_k(x_i, y_i) &\approx 2\pi f \sum_{j=0}^{i} \frac{\Delta z}{c(x_i, y_i, z_j)} \\ &= 2\pi f \Delta z \sum_{m} B_{m} \sum_{j=0}^{i} \rho^{m}(x_i, y_i, z_j) \quad (22) \end{aligned} $$

Indeed, when thermal therapy is performed with ultrasound, they observed changes in the speed of sound inside the skull due to temperature changes; thus this formula allows adapting the phase shift during therapy. However, this polynomial relationship seems to go against all previous models that proposed a relationship proportional to density.

### 3.2.3 | Attenuation

#### Constant attenuation

Deffieux and Konofagou¹⁰ considered a constant attenuation across the whole skull as they claimed that more complex models for attenuation can be inconsistent.

#### Porosity-based and normalized density-based relationships

Similarly to their density and sound speed relationships, Aubry et al.¹⁸ proposed a bone porosity-based

relationship for attenuation:

$$ \alpha = \alpha_{min} + (\alpha_{max} - \alpha_{min}) \psi^{\beta} \qquad (23) $$

The constants were adjusted by comparing simulations with measurements. This method is used in Refs. [11–13, 15, 19–21, 48–51].

Yoon et al.⁵⁶ used a method based on the normalized density, as was done by Deffieux and Konofagou¹⁰ for density and speed of sound:

$$ \alpha = \alpha_{min} + (\alpha_{max} - \alpha_{min}) \rho_{CT} \qquad (24) $$

Unlike the relationships for density and speed of sound proposed with the normalized apparent density, this formula is not equivalent to the porosity-based attenuation formula. Indeed, this is a linear relationship and not a power law, and $\rho_{CT} = \frac{H}{\max(H)} = 1 - \Psi \neq \Psi$.

#### Bone fraction-based relationship

Vyas et al.⁵⁴ proposed a method based on the fraction of bone in the voxel, which is defined by

$$ f = \frac{\frac{H}{1000} \times \frac{\mu}{\rho_{water\_eff}}}{\rho_{bone} \times \frac{\mu}{\rho_{bone\_eff}} - \frac{\mu}{\rho_{water\_eff}}} \qquad (25) $$

With $\frac{\mu}{\rho_{eff}} = \frac{\Sigma S(E) \times \frac{\mu}{\rho}(E)_{water}}{\Sigma S(E)}$, where $S(E)$ is the beam spectrum as a function of energy. Then the absorption component of attenuation is defined by

$$ \alpha_a = f \times \alpha_{bone} + (1 - f) \alpha_{soft tissue} \qquad (26) $$

For the scattering component of the attenuation, they used data found experimentally by Tavakoli.⁵⁷ This data consists in experimental values of attenuation for various porosity values. They fit the data to linear curves and obtained four linear equations (for four ranges of porosity) relating the scattering component of attenuation with porosity. Then they used those equations to compute the scattering component of attenuation, and they used the voxel bone fraction to compute porosity:

$$ \Psi = \frac{f_{maxsubject} - f}{f_{maxsubject}} + 0.05 \qquad (27) $$

#### Relationships obtained with genetic algorithms

Pichardo et al.³⁶ used a genetic algorithm to find a relationship between attenuation and density at several frequencies: 0.27, 0.836, 1.402, 1.965, and 2.525 MHz. Similarly to what they did for the speed of so, they found the relationship that minimizes the phase difference between simulation and experiment. The attenuation
functions found in this study were used (either directly or interpolated to fit the needed frequency) in Refs. [6, 7, 28, 37, 38, 45, 46].

### Polynomial relationship

McDannold et al.²³,²⁴ assumed that the relationship between the attenuation and the skull density can be approximated by a series of polynomials:

$$ \alpha = \sum_{m} A_{m} \rho^{m} \qquad (28) $$

The advantage of such a formulation is that, assuming that attenuation mainly occurs along the propagation axis z, one can compute the pressure field without attenuation and deduce an approximation of the pressure field with a given attenuation model without having to run again the whole simulation. Indeed, if attenuation mainly occurs along the z direction, the pressure distribution for element k is approximately

$$ P_k(x_i, y_i, z_i) \approx P_{k0}(x_i, y_i, z_i) \times \exp\left(-\sum_{j=0}^{i} \alpha(x_i, y_i, z_j) \Delta z\right) \qquad (29) $$

where $P_{k0}$ is the pressure computed without attenuation for the transducer element number k. This formula shows that the pressure field without attenuation can be computed separately from the attenuation contribution with this method, allowing to compare several attenuation models without any loss of computation time.

### Frequency dependency of attenuation

In Pichardo et al.,³⁶ the relationships between density and attenuation found at several frequencies suggest that attenuation generally increases with frequency. Several studies have focused on the relationship between attenuation and frequency.

### Linear relationships.
Connor et al.⁴⁴ used linear increase of attenuation with frequency with various coefficients, one for cortical bone:

$$ \alpha = 167 \times f \times 10^{-6} \text{ Np m}^{-1} \qquad (30) $$

And one for trabecular bone:

$$ \alpha = 300 \times f \times 10^{-6} \text{ Np m}^{-1} \qquad (31) $$

where trabecular bone voxels are identified as those containing less than 70% of bone.

Bossy et al.⁵⁸ performed acoustic simulations through trabecular bone at 0.4–1.2 MHz. They found that the

attenuation linearly increases with frequency and can thus be modeled by

$$ \alpha = \alpha_0 + nBUA \times f \qquad (32) $$

They also found that nBUA values are strongly and positively correlated with the bone volume fraction (BV/TV) and that speed of sound is positively correlated with the bone volume fraction.

Haïat et al.⁵⁹ confirmed that both broadband ultrasonic attenuation (BUA) and speed of sound exhibit a strong positive correlation with the bone volume fraction. In addition, they claimed that BUA and speed of sound vary quasi-linearly with the bone volume fraction.

### Power laws.
Attenuation is a combination of absorption and scattering. According to Pinton et al.,¹⁴ 86% of attenuation is due to scattering (and thus only 14% is due to absorption) in skull bone at 1 MHz. Thus, a more efficient model of attenuation could be defined, by separating the effects of absorption and scattering. That is what is done in Yousefian et al.,⁶⁰ where they numerically studied attenuation in porous media (mimicking cortical bone). They claimed that total attenuation is not described by a linear combination of scattering and absorption anymore in the presence of multiple scattering. Thus, they proposed a nonlinear formula for total attenuation:

$$ \alpha_{tot} = \alpha_{scat} + c f^{\beta_{app}} \qquad (33) $$

where $\beta_{app}$ and c are to be determined and depended on pore diameter and pore density.

Webb et al.⁶¹ compared CT and MR images as predictors of attenuation. They measured the acoustic attenuation at 0.5, 1, and 2.25 MHz in 89 samples taken from two ex vivo human skulls and found the best parameters to fit the equation:

$$ \alpha = \alpha_0 f^{\beta} e^{cp} \qquad (34) $$

where f is the frequency, p the imaging parameter (HU value, UTE and ZTE magnitude, or $T_2^*$ value), and $\alpha_0$, $\beta$, c the model parameters. This equation assumes the traditional power law between frequency and attenuation, and that the imaging parameter provides a rough estimate of the pore structure in each sample. HU provide the best prediction of attenuation with a minimum standard error of 1.7 Np cm⁻¹. The ZTE, UTE, and $T_2^*$ values had standard errors of 2.0, 2.0, and 2.1 respectively.

### 3.2.4 | Deriving acoustical properties from MRI

Although acoustical properties are most often derived from CT-scans, several studies have investigated the
possibility of using MR images (as MR is often used to guide the focused ultrasound therapy), in order to reduce the whole procedure and to avoid the patient exposure to X-rays. These methods consist in building a virtual CT data from MRI.

Wintermark et al.⁶² compared three MRI pulse sequences: T₁-weighted 3D volumetric interpolated breath-hold examination (VIBE), proton density- weighted 3D sampling perfection with application-optimized contrasts using different flip-angle evolution (SPACE), and 3D true fast imaging with steady-state precession T₂-weighted imaging. The MR modality giv- ing a total thickness, the closest to the CT-based total thickness, was identified. Random coefficient regression was used to predict CT total skull thickness based on the optimal MRI sequence. The same method was used for the thickness of each of the three layers (inner table, diploe, and outer table). Similarly, a regression model was used to predict reference standard CT average den- sity based on the optimal MRI sequence average inten- sity. Virtual CT datasets derived from the MRI datasets were built using the models described before. The T₁ 3D VIBE sequence was the MRI sequence coming closest to the reference standard CT in terms of measuring total skull thickness and was thus selected for the rest. In experiments made with a human skull, the mean absolute difference between the phase shifts calculated with standard CT and virtual CT was 0.8 ± 0.6 rad.

Miller et al.⁴⁰ investigated the feasibility of using UTE MRI instead of CT to calculate and apply aberration corrections on a clinical TcMRgFUS system. A 3D map of the skull is created from the MR images by apply- ing a segmentation algorithm and is used to construct virtual CT data, by assigning 1000 HU to bone pixels and −1000 HU to all other pixels. The virtual CT data was entered in the planning software of an MRgFUS system to calculate aberration correction phases. There was no significant difference between the sonication results achieved using CT- and MR-based aberration correction.

### 3.2.5 | Shear waves parameters

As the acoustical parameters of shear waves in bone have not been studied a lot, most studies that account for the presence of shear waves assume that the shear wave's speed is half the one of the longitudinal waves, as it is a typical assumption for solids. However, a few studies have looked more in detail into shear wave parameters.

White et al.⁴¹ measured an average speed of sound for a 1.0-MHz longitudinal wave of 2820 ± 40 and 1500 ± 140 m s⁻¹ for shear waves. They found that the shear attenuation coefficient in skull bone is on aver- age higher by 115 Np m⁻¹ than the longitudinal one for the frequency range studied. So, although the speed of

sound of shear waves in the skull is closer to the one of the surrounding media (brain and tissues), which allows a more efficient transmission through the boundaries, the high shear attenuation coefficient tends to restrict this beneficial effect.

Pichardo et al.⁴⁷ used a genetic algorithm method (as they did before³⁶ for compressional waves parameters) to establish the relationship between apparent density calculated from CT scans and shear speed of sound and attenuation at 270 and 836 kHz. They assume a linear relationship for speed of sound and a constant attenuation.

### 3.2.6 | Comparison of the methods and discussion

*Comparison of a few methods for acoustical parameters modeling*

McDannold et al.²³ compared Pichardo relationships³⁶ for the speed of sound and attenuation to their poly- nomial formulas at 660 kHz. Although there was no significant difference between their relationships and Pichardo's in the resulting dimensions or obliquity of the simulated focal region, the temperatures simulated using Pichardo relationships were less accurate (less close to the experiments) than those simulated using their relationships.

Bancel et al.⁵³ used Marsac density relationship⁵² and compared the speed of sound relationships given by Pichardo,³⁶ Marsac⁵² and McDannold.²³ Marsac rela- tionships seem to perform a bit better for aberration correction simulations, with a restored pressure, com- pared to the hydrophone-based correction, of 85% (compared to 81% with Pichardo relationship and 82.5% with McDannold's relationship) and a peak to side lobe ratio of 45.5% (compared to 51.5% with Pichardo or McDannold relationship).

*Impact of the error in skull modeling*

In Robertson et al. paper,⁴⁸ the impact of changes in bone layer geometry and the speed of sound, density, and acoustic absorption values is quantified through a numerical sensitivity analysis. The errors in the field resulting from noisy variations (around a mean value) in medium properties are smaller than those from linear changes (variations of the mean value). Noisy variations for absorption have almost no effect: Homogeneous maps of absorption are good. Linear variations in HU values resulted in errors lower than the errors for indi- vidual changes in assigned medium properties, thus suggesting that the primary concern should be the robustness of any method for the conversion of CT images to acoustic properties. Speed of sound is shown to be the most influential acoustic property and must be defined with less than 4% error. Changes in the skull thickness of as little as 0.1 mm can cause an error in
peak pressure of greater than 5%, while smoothing with a 1-mm³ kernel (skull maps obtained from low resolution images) causes an increase of over 50% in peak pressure.

### Considering X-ray energy

Although many studies have investigated the relationship between acoustical properties and Hounsfield units, they have ignored the impact of X-ray energy (HU are a function of the linear attenuation coefficient that itself depends on the X-ray energy), reconstruction method, and reconstruction kernel on the measured HU. Webb et al.⁶³ imaged 91 human skull fragments by 80 CT scans with a variety of energies and reconstruction methods. The average HU from each fragment is found for each scan and correlated with the speed of sound measured. The results show that both the energy and the reconstruction method have a significant influence on the relationship between velocity and HU. The main issue is that it is difficult to have an estimate of the real X-ray energy except in dual-energy CT. Aubry and Marsac relationships are based on porosity, which depends on the linear attenuation and thus accounts for X-ray energy. So these formulas work well when the energy is known. Pichardo and Connor relationships do not account for X-ray energy. Considering CT parameters can be obtained with dual-energy CT, by obtaining calibration measurements or by working with the vendor to obtain accurate estimates of the X-ray spectrum of a given scan. The results provide estimates of the relationship between HU and velocity for a variety of different CT parameters. They also show that the relationship between HU and velocity is patient-specific. The measurements show that CT is able to capture only about one-half of the variation in acoustic velocity within the skull. Some of the remaining variation is likely due to errors in the measured velocity, but it is also likely an indication that some of the variation in velocity is not well modeled by HUs. This could be because the variation in velocity results from changes in the chemical composition of the skull, a variation not necessarily captured by the measured HU.

## 3.3 | Level of heterogeneity

The skull is made of two types of bone: cortical bone and trabecular bone. The inner and outer tables are made of cortical bone, which is a dense and nearly homogeneous medium, whereas the middle layer is made of trabecular bone, which is a complex porous structure, highly scattering. Thus, the skull is very heterogeneous. However the imaging modalities, be it CT or MR, cannot capture heterogeneities below their resolution. Thus, when modeling the skull, one can use levels of heterogeneity (Figure 1) going from assign-

ing different acoustical properties to each voxel, to assigning the same acoustical properties to the whole skull.

### 3.3.1 | Fully heterogeneous

In this model, each voxel has a given density, sound speed, and attenuation. It is used in Refs. [5, 6, 10, 13, 15, 18–21, 23–26, 28, 36, 44–46, 48–52, 56, 64–66].

### 3.3.2 | Heterogeneous but binary

The skull is assumed to be a complex porous structure made of bone and with pores filled with marrow (whose acoustical properties are similar to those of water). Simple image thresholding allows labeling voxels as either being bone or marrow and constant acoustical properties are assigned to bone or marrow. It is used in Refs. [14, 17, 29, 58, 59, 67]. However, this method requires high resolution scans, as the structure of the diploe is not visible at the millimeter scale.

### 3.3.3 | N layers

In this model,³⁶ the skull is divided into *N* homogeneous layers. Density along the acoustic axis (only normal incidence is considered) was computed by averaging the density from the voxel on the acoustic axis, and the one of its four neighbors in the plane perpendicular to the acoustic axis. In a given layer, the acoustical properties are averaged from the CT scans. Each layer thickness is set to one fourth of the voxel spacing to allow a smooth variation.

### 3.3.4 | Three layers

The skull is often described as a three-layered medium representing the inner and outer tables made of cortical bone and the middle layer made of trabecular bone. In each layer, the acoustical properties are averaged from the CT scans.³⁰,⁶²

### 3.3.5 | One layer

The simplest assumption that can be made is considering the skull as a homogenous and isotropic medium. The acoustical properties are either averaged from CT scans or taken from measurements from the literature. This model is used in Refs. [6, 20, 31–34, 40, 48, 51].

In Jiang et al.,⁵¹ the velocity on the skull surface is taken to be that in cortical bone (which influences refraction computations) and the velocity inside the skull is
<figure><img src="image_15.png" /><figcaption>**FIGURE 1** The different levels of homogenization: (a) original skull, (b) heterogeneous model, (c) binary model, (d) N-layered model, (e) three-layered model, and (f) one layered model</figcaption></figure>

taken as being the averaged value on the ray paths (which influences the time-of-flight computations). The density is taken as being the maximum density as it influences the refraction computations.

The main issue with those methods is that if the transducer has a large aperture, the acoustical properties may vary quite a lot along the surface intersected by the ultrasound beam. In Jones et al.³⁸ and Pajek and Hynynen,³⁷ a locally averaging method was proposed. It consists in spatially averaging the acoustical properties independently for each ray path going through the skull.

### 3.3.6 | Comparison of the levels of heterogeneity

Jones and Hynynen⁶ experimentally compared CT-based aberration corrections with two skull models: one homogeneous layer and a fully heterogeneous model. For both models, they compared the corrections computed using a full-wave method and using a ray-tracing method. By averaging the metrics on the two methods, the homogeneous model increases the peak pressure shift error (compared to the hydrophone-based shift) by 4%, the -3-dB volume error by 27.5% and the peak-to-side-lobe ratio error by 22%, compared to the heterogeneous model. Thus, the heterogeneous model

performs better, as the homogeneous one underestimates attenuation and the focal spot spreading.

Kyriakou et al.²⁰ found similar results when comparing heterogeneous and homogeneous skull modeling with a full-wave simulation. Indeed, the prediction of the focusing (shift and volume) seems unaffected by skull inhomogeneity, whereas the peak pressure is reduced by 25% with the heterogeneous model.

Jiang et al.⁵¹ compared the performances of two homogenous models: one with the average properties of the skull, and the other with the properties of the cortical bone on the surface (as refraction computations are influenced only by the properties on the surface on the skull) and with the average properties of the skull on the inside (for the time-of-flight computations). Shifts of 9.5 mm in the axial direction and 1.5 mm in the focal plane compared with the virtual source were observed with the first one, compared to 0.5 mm in the axial direction and 0.5 mm in the focal plane with the optimized model.

Robertson et al.⁶⁷ studied the influence of homogenization of CT maps. Progressive homogenization of acoustic property maps leads to an overestimation in the amplitude of transmitted US, an underestimation of time-of-flight, and a loss of fine spatial details in the intracranial field. Inflating the simulated attenuation coefficient of the skull layer reduces the error in transmitted pressure amplitude to around 40%; however, this
<figure><img src="image_16.png" /><figcaption>FIGURE 2 Three-layered model of wave of transmission through the skull</figcaption></figure>

is unable to correct fully for errors in time-of-flight and the pressure distribution of the transmitted field.

One must also take into account that the frequency will influence the level on heterogeneity needed. Indeed, at lower frequencies, heterogeneities much smaller than the wavelength do not need to be taken into account and parameters at those scales can be homogenized.

# 4 SIMULATION METHODS

## 4.1 Semi-analytical methods

### 4.1.1 Ray tracing

A fast common method to simulate wave propagation is ray tracing.32–34,37,38,41,51 It considers reflection, refraction, and mode conversion at each interface. The ray paths are obtained using Snell-Descartes’ laws of a three-layered model, composed of water, skull, and brain. This method is illustrated by Figure 2.

In general, possible reflections from inside the skull are not considered due to the high attenuation of the skull.34 Moreover, as the skull thickness varies across its surface, it acts as a random phase aberration to the transmitted waves, thus making constructive interferences less probable.

Clement et al.32 computed the amplitudes of the reflected wave, the transmitted longitudinal wave, and the transmitted shear wave, with respect to the incident

wave amplitude using the fact that the normal component of the particle displacement, the normal stress, and the shear stress must be continuous at the boundary.

In White et al.,41 the skull interfaces are supposed to be parallel and the pressure (longitudinal and shear) transmission coefficients are computed with

$$ \left\{ \begin{aligned} T_L &= \frac{\rho_f}{\rho_S} \frac{2Z_L \cos(2\theta_S)}{Z_L \cos^2(2\theta_S) + Z_S \sin^2(2\theta_S) + Z_f} \\ T_S &= -\frac{\rho_f}{\rho_S} \frac{2Z_S \cos(2\theta_S)}{Z_L \cos^2(2\theta_S) + Z_S \sin^2(2\theta_S) + Z_f} \end{aligned} \right. \quad (35) $$

where $\rho_f$ and $\rho_S$ are the densities of the fluid and the solid, $Z_L$ and $Z_S$ are the longitudinal and shear impedances of the solid, $Z_f$ is the fluid impedance. The ultrasound beam is approximated as a one-dimensional ray. The values obtained after the first interface are then used as the input for the second interface.

Yin and Hynnenen33 divided each skull interface into small rectangular planar patches. The plane-wave theory is then applied on each patch and the contribution of each patch is then summed using the Rayleigh-Sommerfeld integral. The computations are made for each interface successively: once the computations are obtained for an interface, this interface is considered the new source for the next interface. As the patch surfaces are small (about a quarter wavelength), the acoustic pressure in tissue can be computed as the product of the particle normal velocity and the specific acoustic impedance of the tissue.

Pichardo et al.34 also divided the interfaces into small elements, but the amplitude computations are different. The contribution coefficients due to the longitudinal propagation from the transducer surface to the outer face of the skull $U_{LL,1 \to 2}$, from the outer face to the inner face $U_{LL,2 \to 3}$ and from the inner face to a point $q$ inside the brain $U_{LL,3}$ are computed with

$$ \left\{ \begin{aligned} U_{LL,m \to n} &= \frac{j k_{lm}^* \rho_m c_{lm} e^{-j k_{lm}^* R_{lm} \to q}}{2\pi R_{lm} \to q} \left( 1 - \frac{j}{k_{lm}^* R_{lm} \to q} \right) \\ T_{lm \to in} \cos(\theta_{lm \to in}) ds_{lm} \\ U_{LL,3} &= \frac{j k_{l3}^* \rho_3 c_{l3} e^{-j k_{l3}^* R_{l3} \to q}}{2\pi R_{l3} \to q} ds_{l3} \end{aligned} \right. \quad (36) $$

where $c_{lm}$, and $k_{lm}^* = k_{lm} + j \alpha_{lm}$ are the longitudinal speed of sound and wave number ($\alpha_{lm}$ being the attenuation coefficient), $\rho_m$ is the density, $R_{lm \to in}$ is the distance between the surface elements $i_m$ and $i_n$, $\theta_{lm \to in}$ is the incident angle, $ds_{lm}$ is the surface area of element $i_m$, and $T_{lm \to in}$ is the transmission coefficient. Similarly, the contribution coefficients due to shear-mode conversion are
computed. The pressure amplitude of the transmitted wave is then obtained using the Rayleigh–Sommerfeld integral for a multilayered case:

$$p(q) = p_{LL}(q) + p_{SL}(q) \quad (37)$$

With $p_{LL}(q) = \sum_{i_1=1}^{N_1} u_{i_1} U_{LL,1 \to 2} U_{LL,2 \to 3} U_{LL,3}$ and $p_{SL}(q) = \sum_{i_1=1}^{N_1} u_{i_1} U_{LS,1 \to 2} U_{LS,2 \to 3} U_{LL,3}$, where $N_1, N_2, N_3$ correspond to the number of element of each interface (transducer, outer and inner surfaces of the skull) and $u_{i_1}$ is the amplitude of the particle displacement of the element $i_1$.

Pajek and Hynynen³⁷ used the same method, but shear waves are not taken into account due to their high attenuation at the high frequencies used (higher than 1 MHz). Jones et al.³⁸ also used this method, but with spatially averaged acoustical properties for each ray path going through the skull.

In Jiang et al.,⁵¹ the frequency deviation caused by the acoustic attenuation is calculated in the frequency domain and the ultrasound amplitude after refraction is computed using the transmission coefficient at each interface. After the first interface, the amplitude of the velocity potential is given by

$$Q_2 = Q_s T_r \frac{l_2''}{l_2'' + l_{2r}} \quad (38)$$

where $Q_s$ is the velocity potential of the source, $T_r$ is the plane-wave transmission coefficient, $l_1$ is the length of the incident ray, $l_{2r}$ is the length of the refracted ray, and $l_2'' = \frac{l_1 \sin\alpha \cos^2\beta}{\sin\beta \cos^2\alpha}$ represents the divergence of the beam. Then after propagation through the medium layer, the velocity potential at the second interface is considered a new source to compute the velocity potential after crossing the second interface.

### 4.1.2 | Transfer function methods

#### Transfer matrix method for transmission coefficient computation

In Pichardo et al.,³⁶ as the skull is divided into N layers and is surrounded by water, the transmission coefficient $T$, at normal incidence, is computed using

$$T = \frac{2Z_w}{(m_{2,2} + m_{3,3} + Z_w m_{2,3}) Z_w + m_{3,2}} \quad (39)$$

where $Z_w = \rho_w c_w$ is the acoustic impedance of water, $m_{j,k} = A_{j,k} - \frac{A_{j,1}A_{4,k}}{A_{4,1}}$. ($A_{j,k}$ is the coefficient of the transfer matrix $A$ giving the stress and particle velocity at the last interface from the stress and particle velocity at the first interface. $A$ is obtained by the product of the

matrices of the sound transmission coefficients of the N skull layers.)

#### Projection algorithm in the wavevector-frequency domain

Clement and Hynynen³¹ proposed a projection algorithm in the wavevector-frequency domain that works for normal incidence or any angle of incidence. By taking the Fourier transform of the linearized Stoke's equation, they obtained a simple equation whose solution is given by

$$\tilde{p}(\tilde{k}, z) = \tilde{p}(\tilde{k}, z_0) e^{i\tilde{k} \cdot (z-z_0)} \quad (40)$$

Thus, the field recorded in a plane $z = z_0$ is related to the field at any other plane by a transfer function. The pressure field is then obtained by taking the inverse Fourier transform. For non-normal incidence, a rotation of the planar field at $z = z_0$ is performed in the wavevector-frequency domain.

#### Hybrid angular spectrum (HAS) method

The traditional angular spectrum method assumes that the medium is homogeneous. The pressure in the initial plane is encoded into a spectrum of traveling plane waves in the spatial-frequency domain using fast Fourier transform. Then, the propagation of the waves to the next plane is calculated in the spatial-frequency domain by multiplying the initial spectrum by a propagation transfer function. Finally, IFFT is performed to obtain the pressure in the final plane.

Vyas and Christensen⁶⁸ extended this method to calculate linear ultrasound wave propagation in inhomogeneous tissue geometries. In this new method, wave propagation through each plane is calculated in two steps: one in the space domain and the other in the spatial-frequency domain.

The transmission from one plane to another is computed using

$$p_n(x, y) = p_{n-1}(x, y) \exp(jb_n(x, y) r' - a_n(x, y) r) \quad (41)$$

where $b_n(x, y) = \frac{2\pi f}{c_n(x, y)}$ is the propagation constant at a given voxel $(x, y)$ with a speed of sound $c_n(x, y)$, and $a_n(x, y)$ is the pressure attenuation coefficient of a given voxel. $r'$ is the effective path length between two wavefronts and $r$ is the entire propagation path.

By dividing the phase change into an average phase shift $b'_n r'$ inside a plane and the difference $\Delta b_n(x, y) r'$ from this average phase shift for each voxel of the plane, the transmission equation can be rewritten as

$$p_n(x, y) = p_{n-1}(x, y) \exp(jb'_n r' + j\Delta b_n(x, y)r' - a_n(x, y)r) \quad (42)$$
<table><caption>**TABLE 1** Effects taken into account in the equations used in finite differences studies</caption><thead><tr><th>Reference</th><th>Nonlinearity</th><th>Shear waves</th><th>Absorption</th></tr></thead><tbody><tr><td>19, 20, 26, 29, 36</td><td>X</td><td></td><td>X</td></tr><tr><td>6, 10, 13, 18, 20, 28, 44, 45, 52, 56, 69</td><td></td><td></td><td>X</td></tr><tr><td>21</td><td></td><td></td><td></td></tr><tr><td>58</td><td></td><td>X</td><td></td></tr><tr><td>14</td><td></td><td>X</td><td>X</td></tr><tr><td>17</td><td>X</td><td>X</td><td>X</td></tr></tbody></table>

Propagation from one plane to another then is achieved in the space domain using

$$p'_{n}(x, y) = p_{n-1}(x, y) \exp(j\Delta b_{n}(x, y) r' - a_{n}(x, y) r) \quad (43)$$

And propagation across a given plane is performed in the spatial-frequency domain using:

$$p_n(x, y) = F^{-1} \left[ F[p'_{n}(x, y)] \exp \left( j b'_{n} \sqrt{1 - \alpha^2 - \beta^2 \Delta z} \right) \right] \quad (44)$$

With $\alpha = \lambda f_x$ and $\beta = \lambda f_y$.

This method is used in Refs. [25, 46, 54].

## 4.2 | Numerical methods

### 4.2.1 | Finite differences

Many studies<sup>6,10,13,14,17–21,26,28,29,36,44,45,49,52,56,58,64,65,69,70</sup> used finite difference methods to simulate wave propagation. This allows modeling the skull as a heterogeneous medium but it is computationally intensive. Unlike semi-analytical methods, finite differences methods need a spatial step fine enough in order to converge, as the pressure field at a given position depends on the pressure field at previous positions.

Several wave equations can be considered when using finite differences, depending on which effects one wants to take into account or neglect, such as linearity, shear mode conversion, absorption, and others. The effects taken into account by several studies are listed in Table 1. All those equations are derivatives of the Westervelt equation, which is given by

$$\frac{\partial^2 p}{\partial z^2} - \frac{1}{c^2} \frac{\partial^2 p}{\partial t^2} + \frac{\mu}{c^4} \frac{\partial^3 p}{\partial t^3} + \frac{2\beta}{\rho c^4} \left[ p \frac{\partial^2 p}{\partial t^2} + \left( \frac{\partial p}{\partial t} \right)^2 \right] \\ - \frac{\partial p}{\partial z} \cdot p \frac{\partial (\ln p)}{\partial z} = 0 \quad (45)$$

where $p$ is the wave pressure, $c$ is the local speed of sound, $\rho$ is the local density, $\mu = 2 \alpha c^3 (2\pi f)^{-2}$ is the local diffusivity parameter ($\alpha$ is the local attenuation coefficient) and $\beta$ is the nonlinearity parameter.

### 4.2.2 | k-Space

In 2012, Jing et al.<sup>50</sup> proposes to use a k-space method to simulate wave propagation inside the skull. Like finite differences, it allows a heterogeneous description of the skull, but k-space is faster as it needs less spatial resolution to achieve convergence with a reasonable accuracy.

The k-space methods are part of the bigger family of pseudo-spectral finite difference time domain methods (PSTD). The PSTD methods transform the fields to the spectral domain at each time step in order to compute the spatial derivatives more easily. Those methods are called global methods because the entire simulation domain is used in approximating the derivative of a single point. Those methods are therefore more accurate than standard finite difference methods and theoretically allow discretizing the domain as low as two points per wavelength. However, the standard PSTD methods only improve the efficiency in the spatial domain, as finite differences schemes are still needed in the time domain, and thus fine time steps are needed to avoid dispersion. The k-space method aims to allow larger time steps without lowering the accuracy too much. It consists in multiplying the time step by the k-space operator $\kappa = \text{sinc}(c_{\text{ref}} k \frac{\Delta t}{2})$, which is derived from the exact solution of the homogeneous and lossless wave equation.

The k-space method has been used in other studies,<sup>5,11,15,23,24,48,53,65–67</sup> thanks to the development of the Matlab toolbox k-wave, even if some studies have developed their own k-space algorithm.<sup>51</sup>

## 4.3 | Hybrid methods

A good balance is to combine numerical methods and analytical ones. For instance, one can use an analytical method to compute the wave propagation before and after the skull (as those media can reasonably be assumed homogeneous) and then use a numerical method only within the skull.

When simulations are used for phase correction computations, the waves are propagating from a virtual source inside the brain toward the transducer outside the skull. In a few papers,<sup>52,53,69,71,72</sup> the finite difference or k-space algorithms are only performed from an arbitrary plane or spherical surface as close as possible to the inner surface of the skull, and propagation in brain tissues is modeled by a ray-tracing code to save time. In Refs. [49, 53, 72], the numerical simulations are only performed to an arbitrary plane or a receiving
surface just in front of the skull, and the propagation from there to the transducer is realized by a ray tracing code.

In direct simulations (from the transducer to the brain), hybrid methods are also used. In Wu et al.,⁶⁶ they developed a method that combines the *k*-space-corrected PSTD method with an acoustic holography approach based on the Rayleigh integral. The *k*-wave stage is used to calculate the sound field in the skull, which is divided into multiple sound paths (one per array element). The propagation in each sound path is simulated with *k*-wave from the element to a virtual sensor located near the skull–brain interface. The simulations in all sound paths are run in parallel. The relative amplitude and phase of each discrete point in the sensors are extracted by taking an FFT of an integer number of cycles of the signal after it reaches steady state. Propagation from the virtual sensors is then performed with the Rayleigh integral based method. In Pulkkinen et al.,²⁷ the propagation from the phased-array element through the water to a parallel plane near the skull is simulated with the Rayleigh integral. Then the propagation through the skull is performed with a finite element code and propagation in the brain is simulated using the angular spectrum method.

## 4.4 | Comparison of the methods

All these methods simulate wave propagation through the skull with a variable degree of realism, depending on the acoustic effects taken into account or neglected. Most of the time analytical methods are faster but less accurate than numerical ones. A few papers have compared some of these methods.

In Jing et al.,⁵⁰ they compared the *k*-space method and the finite differences time domain (FDTD) method. It was found that for very fine spatial resolution (more than 10 grid points per wavelength), these two methods match very well. However, at a low spatial resolution, the *k*-space method was observed to produce considerably less numerical error. The computation time at a fine spatial resolution (7.68 grid points per wavelength) can be reduced by a factor of around 80 using 2.56 grid points per wavelength, without altering too much the resulting field.

Jones et al.⁶ compared full-wave finite difference phase correction with analytical phase correction. The results show that the finite difference method outperforms the analytical one for every metrics (shift decrease of 25%, –3-dB volume decrease of 20%, peak-to-side- lobe ratio increase of 26%, signal-to-noise ratio increase of 37% with the full-wave method compared to the ana- lytical one). However, the FDTD method took 7.8 h, whereas the analytical phase correction took only 27 s.

Kyriakou et al.²⁰ compared distance-based phase cor- rection, ray-tracing based phase correction, and finite

difference-based phase correction. The ray-tracing based phase correction only performs a bit better than the distance based one and gives much worse results than the finite difference-based phase correction (shift ×6, peak pressure /2, focal volume ×5).

Robertson et al.⁶⁵ compared FDTD and *k*-space methods for phase aberration correction. Numerical dis- persion has a serious effect on the accuracy of FDTD scheme, resulting in high temporal sampling require- ments to reduce positional error, whereas for the *k*-space scheme, only three pixel per wavelength (PPW) will serve to limit dispersion sufficiently for transcranial transmission for any stable Courant–Friedrichs–Lewy value. To reduce error in the intensity below 10% follow- ing transcranial transmission, *k*-space scheme requires 4.3 PPW, whereas FDTD requires 5.9 PPW.

Jiang et al.⁵¹ compared FDTD, *k*-space and ray- tracing methods, for phase and amplitude correction. The *k*-space method (4 PPW) has a 0.7% phase error compared with the FDTD method (16 PPW). The ray- tracing method has an amplitude error of 5.35% and a phase error of 1.2% compared to the FDTD method. The *k*-space method took 23h35min, whereas the ray-tracing method took 37 min.

Bancel et al.⁵³ compared *k*-space and ray-tracing methods. The restored pressure compared to the hydrophone-based correction is 85% for the *k*-space method and 83.5% for the ray-tracing method. Similarly, the other metrics show results a bit better with the *k*- space method, but the differences are not as big as other studies have shown before.

In terms of computation speed, it is quite hard to compare the methods for various papers as the condi- tions (CPUs, GPUs, parallel computing, etc.) differ from one paper to another and relevant information (such as domain size, time step) are not always available. Table 2 is an attempt to compare the computation times in several studies.

# 5 | EXPERIMENTAL VALIDATION OF THE SIMULATIONS

Experiments are the only way to evaluate the accuracy of simulations. They are most of the time considered gold standards or references even though they can also be enticed with errors as discussed later. In this section, the classical setups for experimental valida- tion of transcranial acoustic simulations are described and the experiments versus simulations results are compared.

## 5.1 | Experimental setups

The typical setup for transcranial acoustic measure- ments is composed of a transducer and an hydrophone
<table><caption>TABLE 2 Comparison of the computation time of several simulation methods</caption><thead><tr><th>References</th><th>Simulation method</th><th>Skull modeling</th><th>Computer</th><th>Domain size (mm³)</th><th>Number of points</th><th>Time</th><th>Time for 1e9 points</th></tr></thead><tbody><tr><td>31</td><td>Layered wavenumber</td><td>Homogeneous</td><td>1 GHz</td><td></td><td></td><td>5 h</td><td></td></tr><tr><td>18</td><td>Finite differences</td><td>Heterogeneous</td><td>500 MHz</td><td>70 x 10 x 30</td><td>2.1e7</td><td>20 h</td><td>950 h</td></tr><tr><td>32</td><td>Ray-tracing</td><td>Locally homogeneous</td><td>1 GHz</td><td></td><td>1.6e4</td><td>30 s</td><td>521 h</td></tr><tr><td>34</td><td>Ray-tracing</td><td>Homogeneous</td><td></td><td>160 x 160 x 160</td><td>3.2e7</td><td>106 min</td><td>55 h</td></tr><tr><td>21</td><td>Finite differences</td><td>Heterogeneous</td><td></td><td>180 x 180 x 150</td><td>1.5e9</td><td>90 min</td><td>60 min</td></tr><tr><td>64</td><td>Finite differences</td><td>Heterogeneous</td><td></td><td></td><td>2.2e9</td><td>120 min</td><td>54 min</td></tr><tr><td>29</td><td>Finite differences</td><td>Binary</td><td>2.8 GHz</td><td></td><td>1.0e8</td><td>12 h</td><td>120 h</td></tr><tr><td>70</td><td>Finite differences</td><td>Heterogeneous</td><td></td><td></td><td>1.3e9</td><td>3 h</td><td>2h18 min</td></tr><tr><td>37</td><td>Ray-tracing</td><td>Homogeneous per ray</td><td></td><td></td><td></td><td>24 h</td><td></td></tr><tr><td>50</td><td>k-Space</td><td>Heterogeneous</td><td>2.6 GHz</td><td>47 x 100 x 31</td><td>7.5e5</td><td>17 s</td><td>6h20 min</td></tr><tr><td>6</td><td>Finite differences</td><td>Heterogeneous/Homo</td><td></td><td>100 x 130 x 40</td><td>2.1e7</td><td>7h48 min</td><td>371 h</td></tr><tr><td>54</td><td>HAS</td><td>Heterogeneous</td><td>2.7 GHz</td><td></td><td></td><td>20 min</td><td></td></tr><tr><td>15</td><td>k-Space</td><td>Heterogeneous</td><td>2.5 GHz</td><td></td><td>1.9e6</td><td>15 min</td><td>131h30 min</td></tr><tr><td>46</td><td>HAS</td><td>Heterogeneous</td><td></td><td></td><td>4.6e7</td><td>45 min</td><td>16 h</td></tr><tr><td>11</td><td>k-Space</td><td>Heterogeneous</td><td></td><td></td><td>1.2e8</td><td>30 h</td><td>250 h</td></tr><tr><td>67</td><td>k-Space</td><td>Heterogeneous</td><td>3.3 GHz</td><td></td><td>2.7e8</td><td>19h53 min</td><td>74 h</td></tr><tr><td>67</td><td>k-Space</td><td>Binary</td><td>3.3 GHz</td><td></td><td>4.2e6</td><td>9 h</td><td>2143 h</td></tr><tr><td>56</td><td>Finite differences (GPU)</td><td>Heterogeneous</td><td>3.6 GHz</td><td>180 x 200 x 175</td><td>6.3e7</td><td>150 s</td><td>40 min</td></tr><tr><td>25</td><td>HAS (7 CPU cores)</td><td>Heterogeneous</td><td>1.4 GHz</td><td></td><td>1.4e9</td><td>180 s</td><td>128 s</td></tr><tr><td>51</td><td>k-Space</td><td>Heterogeneous</td><td></td><td></td><td>1.3e8</td><td>23 h</td><td>177 h</td></tr><tr><td>51</td><td>Ray-tracing</td><td>Homogeneous</td><td></td><td></td><td>1.3e8</td><td>37 min</td><td>4h40 min</td></tr></tbody></table>

<figure><img src="image_17.png" /><figcaption>FIGURE 3 Typical setup for transcranial acoustic measurements</figcaption></figure>

mounted on a three axis positioning system, as illustrated by Figure 3. A first scan, without skull is often performed as a reference result. Then, a skull or a phantom is placed between the transducer and the hydrophone. Prior to the experiments, the skulls are put in water and degassed (typically for 24 h) in order to remove all the air trapped inside and to ensure the skull pores are filled with water to mimic in vivo conditions. When the skull modeling is deduced from CT scan, CT-scans of the skulls should be performed under water as

well. Likewise, the experiments are made in a tank filled with degassed water.

Most of the time, accurate positioning of the skull is performed using a stereotactic frame, or a home-made frame, which contains fiducial markers. The CT scans of the skulls are acquired with the positioning frame and the position of the fiducial markers is determined from the CT images. Likewise, the fiducial marker positions are recorded with the three axis positioning system,¹⁸,⁴⁵ or with an optical-based position tracking system.⁵⁶ From this, a rigid transformation matrix representing the rotation and translation between the experimental and simulation frames can be computed.³¹ Otherwise, the transducer can be attached to the skull positioning frame⁵,³¹,³⁶,⁴⁹,⁶⁶,⁶⁹,⁷² in order to obtain accurate relative positioning more easily. On the other hand, this method does not allow for different positions of the transducer. When experiments are performed under MR guidance, MR tracking coils can be used to find the transducer location within the MR frame.²⁴ MR images can also be used to find the relative position of the transducer with respect to the skull frame.⁴⁶

A few studies have estimated the positioning error between their experiments and their simulations. Those errors are displayed in Table 3.
<table><caption>TABLE 3 Positioning method and error for measured made in several studies</caption><thead><tr><td>References</td><td>Positioning method</td><td>Positioning error</td></tr></thead><tbody><tr><td>31</td><td>Skull attached to the transducer</td><td>0.5 mm</td></tr><tr><td>18</td><td>Position of the skull markers recorded by the three axis system</td><td>1 mm</td></tr><tr><td>49</td><td>Skull attached to the transducer</td><td>A few millimeters</td></tr><tr><td>69</td><td>Skull attached to the transducer</td><td>0.7 mm</td></tr><tr><td>46</td><td>Relative position between transducer and skull measured with MRI</td><td>0.25 mm</td></tr><tr><td>56</td><td>Position of the skull markers recorded by an optical based position tracking</td><td>0.41 mm</td></tr><tr><td>66</td><td>Skull attached to the transducer</td><td>0.02 mm</td></tr></tbody></table>

O'Reilly et al.⁷³ used high frequency ultrasound measurements to localize the skull surface and register CT data to the ultrasound treatment space. The results show on an average submillimeter (0.9 ± 0.2 mm) displacement and sub-degree (0.8° ± 0.4°) rotation registration errors.

## 5.2 | Performance versus simulation

Though many studies have checked their simulation results through experiments, only a few of them have directly compared the pressure field obtained with experiments with the simulated one. Indeed, most studies have compared simulations and experiments through their application. For instance, many studies experimentally compared the field produced by a multielement transducer with a simulated phase correction with the field produced by a multielement transducer with a hydrophone-based correction (which is considered the optimal phase correction).⁵,¹⁷,¹⁸,³¹,⁴⁶,⁵²,⁵³,⁶⁶,⁶⁹,⁷² Another way to assess the accuracy of transcranial acoustic simulations is to compare the results of heat simulations with MR thermometry data²⁵ or with lesion locations.²⁴

In Bouchoux et al.,⁴⁵ Yoon et al.,⁵⁶ and Bancel et al.,⁵³ they directly compared the acoustic field obtained with experiments with the simulated one. The results are summed up in Table 4.

The phase correction results obtained in several studies are displayed in Table 5.

Almquist et al.⁴⁶ performed experiments with phantoms and human skulls in order to separate the errors coming from the simulation algorithm from the errors coming from the skull modeling (as the acoustical properties of the phantoms are known). The results are presented in Table 6.

Simulation-based correction is able to reach 94% of the optimal correction (hydrophone-based correc-

tion) when the acoustical properties are known (phantom experiments), whereas only 70% when they are unknown and deduced from CT-scans with empirical relationships (skull experiments).

In a few studies,²³,²⁵⁻²⁷ MR thermometry data was compared with the results of thermal simulations. Qualitatively speaking, the predicted and measured temperature fields agree well in all those studies. In Pulkkinen et al.,²⁷ for the seven sonications performed on a single skull, the error between the simulated and measured temperatures at the two locations of interest (focus and tissue near the skull base) is on average 1°C and never exceeds 3°C. The other studies collected data from several patients (at least 9) who had undergone several sonications and were thus able to make statistical analysis of the results, as shown in Table 7.

McDannold et al.²³ also compared the dimensions and obliquity of the heating. They found that the predicted and measured dimensions ($R^2 = 0.62$) and obliquity ($R^2 = 0.74$) of the heating are correlated.

## 5.3 | Reasons for error between experiment and simulation

On the simulation side, as shown by Almquist et al.,⁴⁶ imprecise skull modeling based on CT scans is an important source of mismatch between simulations and experiments. In addition, many approximations, such as neglecting shear waves (especially if the incidence angles are greater than 30°³²), attenuation, and nonlinearity, degrade the accuracy of the results. For instance, in the Ding et al. study,¹⁹ the thermal dose at the focus computed with nonlinear wave propagation was almost twice the linear thermal dose, suggesting that nonlinearity effects have to be taken into account. Jiang et al.⁵¹ compared simulations with and without shear waves for incident angles smaller than 20°. They found a shift of 0.5 mm of the focus when shear waves are neglected, whereas perfect focusing is reached when shear waves are included. The maximum pressure reached with shear waves is approximately 3.65 Pa, whereas that when shear waves are not considered is approximately 4.85 Pa.

On the experimental side, very close positioning of the experiments relatively to the simulations is crucial in order to be sure to be able to compare them. One also needs to keep in mind that hydrophones are generally not very precise (up to 10% errors) and need to be well calibrated.

# 6 | DISCUSSION

Almquist et al.⁴⁶ showed that transcranial acoustic simulations accuracy is highly dependent on skull modeling, which can be divided into two steps: the geometric description and the acoustic properties. The shape is
<table><caption>**TABLE 4** Comparison between experimental and simulated acoustic field in different studies</caption><thead><tr><th>References</th><th>Focus shift</th><th>Pressure metric</th><th>Pressure value</th></tr></thead><tbody><tr><td>45</td><td>Below the scanning step (2.5 mm)</td><td>$$ \sum_{x,y,z} \frac{\max_t |p_{simu} - \max_t |p_{meas}|}{\max_t |p_{meas}|} $$</td><td>-9.17 ± 4.55%, averaged on four skulls (120 kHz)</td></tr><tr><td>56</td><td>1.43 ± 0.8 mm</td><td>$$ \frac{R^{XY} - \overline{R^{XY}}}{R^{XY}} + \frac{R^{XZ} - \overline{R^{XZ}}}{R^{XZ}} + \frac{R^{ZY} - \overline{R^{ZY}}}{R^{ZY}} $$<br>With R the ratio of the peak pressure with and without the skull</td><td>5.8 ± 6.4%, averaged on three sheep skulls (250 kHz)</td></tr><tr><td>53</td><td>0.3 ± 0.1 mm</td><td>$$ \frac{1}{N_{elements}} \sum_{i=1}^{N_{elements}} p_{i}^{simu} - p_{i}^{meas} $$</td><td>-0.9 ± 5.4%, averaged on five skulls (800 kHz)</td></tr></tbody></table>

<table><caption>**TABLE 5** Comparison of simulation-based and invasive aberration correction in several studies</caption><thead><tr><th rowspan="2">References</th><th colspan="3">Shift</th><th colspan="2">Percentage of peak pressure compared to the uncorrected case</th></tr><tr><th>Without correction (mm)</th><th>With simulation-based correction (mm)</th><th>With hydrophone-based correction (mm)</th><th>With simulation-based correction (%)</th><th>With hydrophone-based correction (%)</th></tr></thead><tbody><tr><td>31</td><td>1.1</td><td>0.48</td><td></td><td>135</td><td>294</td></tr><tr><td>17</td><td></td><td>0</td><td></td><td></td><td></td></tr><tr><td>46</td><td>1.77</td><td>0.71</td><td>0.25</td><td>151</td><td>217</td></tr><tr><td>52</td><td>1.05</td><td>0.54</td><td></td><td>190</td><td>220</td></tr><tr><td>66</td><td>0.5</td><td>0</td><td></td><td>137</td><td>147</td></tr><tr><td>53</td><td>0.9</td><td>0.2</td><td>0</td><td>130</td><td>153</td></tr><tr><td>72</td><td></td><td>0.63</td><td></td><td>149</td><td>165</td></tr></tbody></table>

<table><caption>**TABLE 6** Comparison of simulated-based and invasive aberration correction for a phantom and for a human skull</caption><thead><tr><th rowspan="2">Metric</th><th colspan="3">Phantom</th><th colspan="3">Skull</th></tr><tr><th>Without correction</th><th>With simulation-based correction</th><th>With hydrophone-based correction</th><th>Without correction</th><th>With simulation-based correction</th><th>With hydrophone-based correction</th></tr></thead><tbody><tr><td>Shift (mm)</td><td>1.9</td><td>0.56</td><td>0.25</td><td>1.77</td><td>0.71</td><td>0.25</td></tr><tr><td>Percentage of peak pressure compared to the uncorrected case (%)</td><td>100</td><td>141</td><td>150</td><td>100</td><td>151</td><td>217</td></tr></tbody></table>

<table><caption>**TABLE 7** Correlation between simulated and measured temperature in several studies</caption><thead><tr><th>References</th><th>Metric</th><th>Correlation between simulated and measured</th></tr></thead><tbody><tr><td>23</td><td>Peak temperature rise</td><td>r<sup>2</sup> = 0.71</td></tr><tr><td>26</td><td>Peak temperature</td><td>r<sup>2</sup> = 0.75</td></tr><tr><td>25</td><td>Focal spot temperature rise</td><td>0.613 ≤ r<sup>2</sup> ≤ 0.947</td></tr></tbody></table>

not needed for numerical methods and is usually segmented from CT or MR images for semi-analytical methods. Over the years, image-derived acoustic properties have replaced measured ones, because of their case-specific adaptability. Most studies agree that both density and speed of sound depend linearly on HU values, but various empirical relationships have been

proposed and improved. Marsac et al.⁵² seem to have reached the more general porosity-based relationships and their accuracy is demonstrated by Bancel et al.⁵³ The main issue with these relationships remains to choose the constants ($ρ_{max}$, $ρ_{min}$, $c_{max}$, $c_{min}$). Although $ρ_{min}$ and $c_{min}$ are often taken as the values of water, $ρ_{max}$ and $c_{max}$ are sample-dependent and thus hard to choose. Marsac et al.⁵² tried to find the best $c_{max}$ such that simulations fit the experiments, but this kind of approach is hard to generalize as it depends mostly on the samples used, the frequency and the experimental conditions. Indeed, Pichardo et al.³⁶ showed that the speed-of-sound relationships are frequency dependent. A polynomial relationship has been proposed by McDannold et al.²³ for the speed of sound, with constants found by comparing simulations with experiments. Once again, the problem of such strategies is that the relationships are case specific to the experimental
conditions of the study and are not always well generalizable. To avoid any bias, the relationships need to be compared in a different study, such as the one from Bancel et al.,⁵³ who showed that Marsac relationships were more accurate than Pichardo's and McDannold's. But more comparison studies are needed. Finally, attenuation is the parameter that is the hardest to determine as it is not well understood. Various empirical relationships have been proposed, relating attenuation to HU values with a linear or a power law. Robertson et al.⁴⁸ showed that attenuation can be defined by a constant across the whole skull, but of course this constant is skull specific. Linear and power laws relating attenuation with frequency have also been investigated. In particular, in trabecular bone, where most of the scattering occurs, linear relationships have been established with a slope that varies quasi-linearly with the bone fraction volume. However, Yousefian et al.⁶⁰ claimed that total attenuation is not described by a linear combination of scattering and absorption anymore in the presence of multiple scattering. They proposed a relationship for attenuation, which is the sum of a constant absorption and a frequency dependent power law for scattering. However, to our knowledge, no comparison studies have shown that any attenuation relationship is better than the others. Although most studies derived the acoustic properties of the skull from CT scans, Webb et al.⁶³ warned that X-ray energy, reconstruction method, and reconstruction kernel should also be taken into account. After computing the acoustic properties of the skull, one can consider several levels of heterogeneity, going from global properties for the whole skull to individual properties for each voxel. Comparisons of these levels of heterogeneity are complicated as it is closely linked to the kind of simulation method used. Indeed, numerical methods tend to use fully heterogeneous skull model, whereas semi-analytical methods (in particular ray-tracing) often use fully homogeneous ones. Jones and Hynynen⁶¹ and Kyriakou et al.²⁰ both compared homogeneous and heterogeneous models with a same simulation method for correction aberration and found that, as can be expected, the heterogeneous performs better, especially in terms of predicted pressure amplitude, whereas both models seem nearly equivalent in terms of predicted shift. This is in agreement with the Robertson et al. study,⁶⁷ who showed that the homogenization of acoustic property maps leads to an overestimation in the amplitude of ultrasound and an underestimation of time-of-flight. However, smarter homogeneous models can be used to improve computation time without losing too much accuracy. For example, Jiang et al.⁷⁴ showed that taking different averaged properties for the refraction computations and for the time-of-flight computations significantly improves the accuracy without increasing the computation time.

In terms of simulation methods, the most used ones are ray-tracing, hybrid angular spectrum, FDTD and

$k$-space. The $k$-space method seems to have replaced FDTD, as it produces less numerical error at low spatial resolution.⁵⁰ Thus, the $k$-space method can be used with quite coarse meshes without losing to much accuracy, which results in a computation time decrease (by a factor of around 80⁵⁰). Several studies compare finite difference based⁶,²⁰ and $k$-space based⁵¹,⁵³ phase correction with ray-tracing-based phase correction. The results show that both FDTD and $k$-space-based phase correction outperforms ray-tracing-based phase correction, but the gap between numerical methods and ray-tracing varies between the studies (going from a 2% difference⁵³ to a 50% difference²⁰ in restored pressure at focus). In terms of computation time, FDTD is slower than ray-tracing by three orders of magnitude (Jones et al.⁶), whereas the $k$-space method is slower than ray-tracing by one order of magnitude.⁵¹

## 7 | CONCLUSION

Transcranial focused ultrasound is a promising method for several therapeutic applications. Simulations are needed to optimize such treatments and to assess their safety. The main challenge of such simulations is to determine the skull acoustic properties. Then, depending on level of precision of the skull modeling, different types of simulation methods can be used. Another big challenge is to make measurements precise enough to be able to compare them with simulation results and thereby to assess the accuracy of the simulations.

## CONFLICT OF INTEREST

The authors have no conflict of interest to disclose.

## REFERENCES

1. White E. Focused ultrasound foundation state of the field report 2021. 2021.

2. Reinhard M, Hetzel A, Krüger S, et al. Blood-brain barrier disruption by low-frequency ultrasound. *Stroke*. 2006;37(6):1546-1548. https://doi.org/10.1161/01.STR.0000221813.27519.0B

3. Wagner T, Valero-Cabre A, Pascual-Leone A. Noninvasive human brain stimulation. *Annu Rev Biomed Eng*. 2007;9:527-565. https://doi.org/10.1146/ANNUREV.BIOENG.9.061206.133100

4. Yasuda J, Yoshikawa H, Tanaka H. Phase aberration correction for focused ultrasound transmission by refraction compensation. *Jpn J Appl Phys*. 2019;58(SG):SGGE22. https://doi.org/10.7567/1347-4065/ab19aa

5. Maimbourg G, Houdouin A, Deffieux T, Tanter M, Aubry JF. 3D-printed adaptive acoustic lens as a disruptive technology for transcranial ultrasound therapy using single-element transducers. *Phys Med Biol*. 2018;63(2):025026. https://doi.org/10.1088/1361-6560/aaa037

6. Jones RM, Hynynen K. Comparison of analytical and numerical approaches for CT-based aberration correction in transcranial passive acoustic imaging. *Phys Med Biol*. 2015;61(1):23-36. https://doi.org/10.1088/0031-9155/61/1/23
7. Jones RM, O'Reilly MA, Hynynen K. Experimental demonstration of passive acoustic imaging in the human skull cavity using CT-based aberration corrections. *Med Phys*. 2015;42(7):4385-4400. https://doi.org/10.1118/1.4922677

8. Gennisson J-L, Chatillon S. ECHOBRAIN: aberration correction in transcranial ultrasonic imaging using CT data and simulation-based focusing algorithms. *J Acoust Soc Am*. 2020;148:2485.

9. Elias WJ, Huss D, Voss T, et al. A pilot study of focused ultrasound thalamotomy for essential tremor. *N Engl J Med*. 2013;369(7):640-648. https://doi.org/10.1056/NEJMOA1300962

10. Deffieux T, Konofagou E. Numerical study of a simple transcranial focused ultrasound system applied to blood-brain barrier opening. *IEEE Trans Ultrason Ferroelectr Freq Control*. 2010;57(12):2637-2653. https://doi.org/10.1109/TUFFC. 2010.1738

11. Mueller JK, Ai L, Bansal P, Legon W. Numerical evaluation of the skull for human neuromodulation with transcranial focused ultrasound. *J Neural Eng*. 2017;14(6):066012. https://doi.org/10.1088/ 1741-2552/aa843e

12. Zhang Y, Zhang H, Pan T, Jian X. Simulation of standing wave reduction in HIFU transcranial tumors therapy. IEEE Int Ultrason Symp IUS. Octob 2019:552-555. https://doi.org/10.1109/ ULTSYM.2019.8926033

13. Baron C, Aubry JF, Tanter M, Meairs S, Fink M. Simulation of intracranial acoustic fields in clinical trials of sonothrombolysis. *Ultrasound Med Biol*. 2009;35(7):1148-1158. https://doi.org/10. 1016/j.ultrasmedbio.2008.11.014

14. Pinton G, Aubry JF, Bossy E, Muller M, Pernot M, Tanter M. Attenuation, scattering, and absorption of ultrasound in the skull bone. *Med Phys*. 2012;39(1):299-307. https://doi.org/10.1118/1. 3668316

15. Top CB, White PJ, McDannold NJ. Nonthermal ablation of deep brain targets: a simulation study on a large animal model. *Med Phys*. 2016;43(2):870-882. https://doi.org/10.1118/1.4939809

16. Park TY, Pahk KJ, Kim H. Method to optimize the placement of a single-element transducer for transcranial focused ultrasound. *Comput Methods Programs Biomed*. 2019;179:104982. https://doi.org/10.1016/j.cmpb.2019.104982

17. Narumi R, Matsuki K, Mitarai S, et al. Focus control aided by numerical simulation in heterogeneous media for high-intensity focused ultrasound treatment. *Jpn J Appl Phys*. 2013;52(7 pt 2):07HF01. https://doi.org/10.7567/JJAP52.07HF01

18. Aubry J-F, Tanter M, Pernot M, Thomas J-L, Fink M. Experimental demonstration of noninvasive transskull adaptive focusing based on prior computed tomography scans. *J Acoust Soc Am*. 2003;113(1):84-93. https://doi.org/10.1121/1.1529663

19. Ding X, Wang Y, Zhang Q, et al. Modulation of transcranial focusing thermal deposition in nonlinear HIFU brain surgery by numerical simulation. *Phys Med Biol*. 2015;60(10):3975-3998. https://doi.org/10.1088/0031-9155/60/10/3975

20. Kyriakou A, Neufeld E, Werner B, Székely G, Kuster N. Full-wave acoustic and thermal modeling of transcranial ultrasound propagation and investigation of skull-induced aberration correction techniques: a feasibility study. *J Ther ultrasound*. 2015;3(1):11. https://doi.org/10.1186/s40349-015-0032-9

21. Gâteau J, Marsac L, Pernot M, Aubry JF, Tanter M, Fink M. Transcranial ultrasonic therapy based on time reversal of acoustically induced cavitation bubble signature. *IEEE Trans Biomed Eng*. 2010;57(1):134-144. https://doi.org/10.1109/TBME. 2009.2031816

22. Marquet F, Tung YS, Konofagou EE. Feasibility study of a clinical blood-brain barrier opening ultrasound system. *Nano Life*. 2010;1:309. https://doi.org/10.1142/S1793984410000286. FEASIBILITY

23. McDannold N, Jason White P, Cosgrove R. An element-wise approach for simulating transcranial MRI-guided focused ultrasound thermal ablation. *arXiv*. Published online September 5, 2019. https://doi.org/10.1103/physrevresearch.1.033205

24. McDannold N, Jason White P, Cosgrove R. Predicting skull lesions after clinical transcranial MRI-guided focused ultrasound with acoustic and thermal simulations. *arXiv*. 2019.

25. Leung SA, Webb TD, Bitton RR, Ghanouni P, Butts Pauly K. A rapid beam simulation framework for transcranial focused ultrasound. *Sci Rep*. 2019;9(1):7965. https://doi.org/10.1038/s41598- 019-43775-6

26. Guo S, Zhuo J, Li G, et al. Feasibility of ultrashort echo time images using full-wave acoustic and thermal modeling for transcranial MRI-guided focused ultrasound (tCMRGUS) planning. *Phys Med Biol*. 2019;64(9):095008. https://doi.org/10.1088/1361- 6560/ab12f7

27. Pulkkinen A, Huang Y, Song J, Hynynen K. Simulations and measurements of transcranial low-frequency ultrasound therapy: skull-base heating and effective area of treatment. *Phys Med Biol*. 2011;56(15):4661-4683. https://doi.org/10.1088/0031-9155/ 56/15/003

28. Bouchoux G, Shivashankar R, Abruzzo TA, Holland CK. In silico study of low-frequency transcranial ultrasound fields in acute ischemic stroke patients. *Ultrasound Med Biol*. 2014;40(6):1154-1166. https://doi.org/10.1016/j.ultrasmedbio.2013.12.025

29. Pinton G, Aubry JF, Fink M, Tanter M. Numerical prediction of frequency dependent 3D maps of mechanical index thresholds in ultrasonic brain therapy. *Med Phys*. 2012;39(1):455-467. https:// doi.org/10.1118/1.3670376

30. Hayner M, Hynynen K. Numerical analysis of ultrasonic transmission and absorption of oblique plane waves through the human skull. *J Acoust Soc Am*. 2001;110(6):3319-3330. https://doi.org/ 10.1121/1.1410964

31. Clement GT, Hynynen K. A non-invasive method for focusing ultrasound through the human skull. *Phys Med Biol*. 2002;47(8):1219-1236. https://doi.org/10.1088/0031-9155/47/8/ 301

32. Clement GT, White PJ, Hynynen K. Enhanced ultrasound transmission through the human skull using shear mode conversion. *J Acoust Soc Am*. 2004;115(3):1356-1364. https://doi.org/10.1121/ 1.1645610

33. Yin X, Hynynen K. A numerical study of transcranial focused ultrasound beam propagation at low frequency. *Phys Med Biol*. 2005;50(8):1821-1836. https://doi.org/10.1088/0031-9155/ 50/8/013

34. Pichardo S, Hynynen K. Treatment of near-skull brain tissue with a focused device using shear-mode conversion: a numerical study. *Phys Med Biol*. 2007;52(24):7313-7332. https://doi.org/ 10.1088/0031-9155/52/24/008

35. Smith SM, Jenkinson M, Woolrich MW, et al. Advances in functional and structural MR image analysis and implementation as FSL. *Neuroimage*. 2004;23(suppl 1):S208-S219. https://doi.org/ 10.1016/J.NEUROIMAGE.2004.07.051

36. Pichardo S, Sin VW, Hynynen K. Multi-frequency characterization of the speed of sound and attenuation coefficient for longitudinal transmission of freshly excised human skulls. *Phys Med Biol*. 2010;56(1):219-250. https://doi.org/10.1088/0031-9155/56/1/014

37. Pajek D, Hynynen K. The design of a focused ultrasound transducer array for the treatment of stroke: a simulation study. *Phys Med Biol*. 2012;57(15):4951-4968. https://doi.org/10.1088/0031- 9155/57/15/4951

38. Jones RM, O'Reilly MA, Hynynen K. Transcranial passive acoustic mapping with hemispherical sparse arrays using CT-based skull-specific aberration corrections: a simulation study. *Phys Med Biol*. 2013;58(14):4981-5005. https://doi.org/10.1088/0031- 9155/58/14/4981

39. Park TY, Pahk KJ, Kim H. A novel numerical approach to stimulation of a specific brain region using transcranial focused ultrasound. Proc Annu Int Conf IEEE Eng Med Biol Soc EMBS. July 2018:3697-3700. https://doi.org/10.1109/EMBC.2018.8513331

40. Miller GW, Eames M, Snell J, Aubry JF. Ultrashort echo-time MRI versus CT for skull aberration correction in MR-guided
transcranial focused ultrasound: in vitro comparison on human calvaria. *Med Phys.* 2015;42(5):2223-2233. https://doi.org/10.1118/1.4916656

41. White PJ, Clement GT, Hynynen K. Longitudinal and shear mode ultrasound propagation in human skull bone. *Ultrasound Med Biol.* 2006;32(7):1085-1096. https://doi.org/10.1016/j.ultrasmedbio.2006.03.015

42. Duck FA. Physical properties of tissue: a comprehensive reference book. Published online 1990:346.

43. Rho JY, Hobatho MC, Ashman RB. Relations of mechanical properties to density and CT numbers in human bone. *Med Eng Phys.* 1995;17(5):347-355. https://doi.org/10.1016/1350-4533(95)97314-F

44. Connor CW, Clement GT, Hynynen K. A unified model for the speed of sound in cranial bone based on genetic algorithm optimization. *Phys Med Biol.* 2002;47(22):3925-3944. https://doi.org/10.1088/0031-9155/47/22/302

45. Bouchoux G, Bader KB, Korfhagen JJ, et al. Experimental validation of a finite-difference model for the prediction of transcranial ultrasound fields based on CT images. *Phys Med Biol.* 2012;57(23):8005-8022. https://doi.org/10.1088/0031-9155/57/23/8005

46. Almquist S, Parker DL, Christensen DA. Rapid full-wave phase aberration correction method for transcranial high-intensity focused ultrasound therapies. *J Ther Ultrasound.* 2016;4(1):30. https://doi.org/10.1186/s40349-016-0074-7

47. Pichardo S, Moreno-Hernández C, Andrew Drainville R, Sin V, Curiel L, Hynynen K. 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.* 2017;62(17):6938-6962. https://doi.org/10.1088/1361-6560/aa7ccc

48. Robertson J, Martin E, Cox B, Treeby BE. Sensitivity of simulated transcranial ultrasound fields to acoustic medium property maps. *Phys Med Biol.* 2017;62(7):2559-2580. https://doi.org/10.1088/1361-6560/aa5e98

49. Marquet F, Pernot M, Aubry JF, Montaldo G, Tanter M, Fink M. Non-invasive transcranial ultrasound therapy guided by CT-scans. Annu Int Conf IEEE Eng Med Biol - Proc. Published online 2006:683-687. https://doi.org/10.1109/IEMBS.2006.259962

50. Jing Y, Meral FC, Clement GT. Time-reversal transcranial ultrasound beam focusing using a k-space method. *Phys Med Biol.* 2012;57(4):901-917. https://doi.org/10.1088/0031-9155/57/4/901

51. Jiang C, Li D, Xu F, Li Y, Liu C, Ta D. Numerical evaluation of the influence of skull heterogeneity on transcranial ultrasonic focusing. *Front Neurosci.* 2020;14:317. https://doi.org/10.3389/fnins.2020.00317

52. Marsac L, Chauvet D, La Greca R, et al. Ex vivo optimisation of a heterogeneous speed of sound model of the human skull for non-invasive transcranial focused ultrasound at 1 MHz. *Int J Hyperthermia.* 2017;33(6):635-645. https://doi.org/10.1080/02656736.2017.1295322

53. Bancel T, Houdouin A, Annic P, et al. 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.* 2021;68:2554-2565. Published online 2021. https://doi.org/10.1109/TUFFC.2021.3063055

54. Vyas U, Ghanouni P, Halpern CH, Elias J, Pauly KB. Predicting variation in subject thermal response during transcranial magnetic resonance guided focused ultrasound surgery: comparison in seventeen subject datasets. *Med Phys.* 2016;43(9):5170-5180. https://doi.org/10.1118/1.4955436

55. Fry FJ, Barger JE. Acoustical properties of the human skull. *J Acoust Soc Am.* 1978;63(5):1576-1590. https://doi.org/10.1121/1.381852

56. Yoon K, Lee W, Croce P, Cammalleri A, Yoo SS. Multi-resolution simulation of focused ultrasound propagation through ovine skull from a single-element transducer. *Phys Med Biol.* 2018;63(10):105001. https://doi.org/10.1088/1361-6560/aabe37

57. Tavakoli MB. The effect of bone structure on ultrasonic attenuation and velocity. *Ultrasonics.* 1992;30(6):389-395.

58. Bossy E, Padilla F, Peyrin F, Laugier P. Three-dimensional simulation of ultrasound propagation through trabecular bone structures measured by synchrotron microtomography. *Phys Med Biol.* 2005;50(23):5545-5556. https://doi.org/10.1088/0031-9155/50/23/009

59. Haïat G, Padilla F, Peyrin F, Laugier P. Variation of ultrasonic parameters with microstructure and material properties of trabecular bone: a 3D model simulation. *J Bone Miner Res.* 2007;22(5):665-674. https://doi.org/10.1359/jbm.070209

60. Yousefian O, Karbalaeisadegh Y, Muller M. Frequency-dependent analysis of ultrasound apparent absorption coefficient in multiple scattering porous media: application to cortical bone. *Phys Med Biol.* 2021;66(3):35026. https://doi.org/10.1088/1361-6560/abb934

61. Webb TD, Leung SA, Ghanouni P, Dahl JJ, Pelc NJ, Pauly KB. Acoustic attenuation: multifrequency measurement and relationship To CT and MR imaging. *IEEE Trans Ultrason Ferroelectr Freq Control.* 2021;68(5):1532-1545. https://doi.org/10.1109/TUFFC.2020.3039743

62. Wintermark M, Tustison NJ, Elias WJ, et al. T1-weighted MRI as a substitute to CT for refocusing planning in MR-guided focused ultrasound. *Phys Med Biol.* 2014;59(13):3599-3614. https://doi.org/10.1088/0031-9155/59/13/3599

63. Webb TD, Leung SA, Rosenberg J, et al. Measurements of the relationship between CT Hounsfield units and acoustic velocity and how it changes with photon energy and reconstruction method. *IEEE Trans Ultrason Ferroelectr Freq Control.* 2018;65(7):1111-1124. https://doi.org/10.1109/TUFFC.2018.2827899

64. Aubry JF, Marsac L, Pernot M, et al. High intensity focused ultrasound for transcranial therapy of brain lesions and disorders. *IRBM.* 2010;31(2):87-91. https://doi.org/10.1016/j.irbm.2010.02.013

65. Robertson JLB, Cox BT, Jaros J, Treeby BE. Accurate simulation of transcranial ultrasound propagation for ultrasonic neuromodulation and stimulation. *J Acoust Soc Am.* 2017;141(3):1726-1738. https://doi.org/10.1121/1.4976339

66. Wu N, Shen G, Qu X, et al. An efficient and accurate parallel hybrid acoustic signal correction method for transcranial ultrasound. *Phys Med Biol.* 2020;65(21):215019. https://doi.org/10.1088/1361-6560/abaa25

67. Robertson J, Urban J, Stitzel J, Treeby BE. The effects of image homogenisation on simulated transcranial ultrasound propagation. *Phys Med Biol.* 2018;63:145014 Published online 2018.

68. Vyas U, Christensen D. Ultrasound beam simulations in inhomogeneous tissue geometries using the hybrid angular spectrum method. *IEEE Trans Ultrason Ferroelectr Freq Control.* 2012;59(6):1093-1100. https://doi.org/10.1109/TUFFC.2012.2300

69. Marquet F, Pernot M, Aubry JF, et al. Non-invasive transcranial ultrasound therapy based on a 3D CT scan: protocol validation and in vitro results. *Phys Med Biol.* 2009;54(9):2597-2613. https://doi.org/10.1088/0031-9155/54/9/001

70. Marsac L, Chauvet D, Larrat B, et al. MR-guided adaptive focusing of therapeutic ultrasound beams in the human head. *Med Phys.* 2012;39(2):1141-1149. https://doi.org/10.1118/1.3678988

71. Maimbourg G. Thérapie du cerveau par ultrasons focalisés transcraniens Rationalisation du nombre de transducteurs p. Published online 2018.

72. Maimbourg G, Guilbert J, Bancel T, et al. Computationally efficient transcranial ultrasonic focusing: taking advantage of the
high correlation length of the human skull. *IEEE Trans Ultrason Ferroelectr Freq Control.* 2020;67(10):1993-2002. https://doi.org/10.1109/TUFFC.2020.2993718

73. O'Reilly MA, Jones RM, Birman G, Hynynen K. Registration of human skull computed tomography data to an ultrasound treatment space using a sparse high frequency ultrasound hemispherical array. *Med Phys.* 2016;43(9):5063-5071. https://doi.org/10.1118/1.4960362

74. Jiang C, Li Y, Li B, et al. Ray theory-based transcranial phase correction for intracranial imaging: A phantom study.

IEEE Access. 2019;7:163013-163021. https://doi.org/10.1109/ACCESS.2019.2951152

How to cite this article: Angla C, Larrat B, Gennisson J-L, Chatillon S. Transcranial ultrasound simulations: A review. *Med Phys.* 2022;1-22. https://doi.org/10.1002/mp.15955