Skip to main content

Full-wave acoustic and thermal modeling of transcranial ultrasound propagation and investigation of skull-induced aberration correction techniques: a feasibility study



Transcranial focused ultrasound (tcFUS) is an attractive noninvasive modality for neurosurgical interventions. The presence of the skull, however, compromises the efficiency of tcFUS therapy, as its heterogeneous nature and acoustic characteristics induce significant distortion of the acoustic energy deposition, focal shifts, and thermal gain decrease. Phased-array transducers allow for partial compensation of skull-induced aberrations by application of precalculated phase and amplitude corrections.


An integrated numerical framework allowing for 3D full-wave, nonlinear acoustic and thermal simulations has been developed and applied to tcFUS. Simulations were performed to investigate the impact of skull aberrations, the possibility of extending the treatment envelope, and adverse secondary effects. The simulated setup comprised an idealized model of the ExAblate Neuro and a detailed MR-based anatomical head model. Four different approaches were employed to calculate aberration corrections (analytical calculation of the aberration corrections disregarding tissue heterogeneities; a semi-analytical ray-tracing approach compensating for the presence of the skull; two simulation-based time-reversal approaches with and without pressure amplitude corrections which account for the entire anatomy). These impact of these approaches on the pressure and temperature distributions were evaluated for 22 brain-targets


While (semi-)analytical approaches failed to induced high pressure or ablative temperatures in any but the targets in the close vicinity of the geometric focus, simulation-based approaches indicate the possibility of considerably extending the treatment envelope (including targets below the transducer level and locations several centimeters off the geometric focus), generation of sharper foci, and increased targeting accuracy. While the prediction of achievable aberration correction appears to be unaffected by the detailed bone-structure, proper consideration of inhomogeneity is required to predict the pressure distribution for given steering parameters


Simulation-based approaches to calculate aberration corrections may aid in the extension of the tcFUS treatment envelope as well as predict and avoid secondary effects (standing waves, skull heating). Due to their superior performance, simulationbased techniques may prove invaluable in the amelioration of skull-induced aberration effects in tcFUS therapy. The next steps are to investigate shear-wave-induced effects in order to reliably exclude secondary hot-spots, and to develop comprehensive uncertainty assessment and validation procedures.


Transcranial focused ultrasound (tcFUS) under magnetic resonance imaging guidance (tcMRgFUS) has attracted the interest of the scientific and clinical community in recent years as a noninvasive and promising treatment modality. Initial clinical trials have indicated successful treatment of patients with brain tumors [1, 2], neuropathic pain [3, 4], and essential tremor [57]. In addition, clinical trials for the treatment of movement disorders, gliomas, obsessive compulsive disorder (OCD), and Parkinson’s disease are being set up in medical centers around the world [8]. Apart from the neurosurgical applications of FUS, this technology has been extensively evaluated for applications such as thrombolysis [913], blood-brain barrier (BBB) disruption for increased drug delivery [1416], and even neuro-modulation [1721].

Despite the substantial benefits of tcFUS when employed in the clinical setting, complications have been reported in human trials, typically in the form of unforeseen brain hemorrhaging [4, 22, 23]. The skull poses the largest barrier to the use of tcFUS; the combination of its complex heterogeneous nature and its significantly different acoustic properties compared to the soft tissue, i.e., twice the speed of sound and density and at least an order of magnitude higher absorption [24], can cause several undesirable effects. These include distortion of the acoustic energy deposition, shifting of the focus, and a significant decrease of the treatment’s focal gain, i.e., the ratio between the energy deposition at the focus and the energy deposition on the scalp and skull bone [25, 26].

The problems related to skull-induced aberrations have been partially resolved by employing large aperture, hemispherical transducer arrays that feature hundreds to more than 1000 elements, where each element is driven with an individual phase and amplitude [22, 27]. The large number of elements permits the acoustic energy to be distributed on the skull surface, thus diminishing the local deposition on the scalp and bone. In addition, the ability to drive the transducer elements individually with appropriately corrected phases and amplitudes allows for focal distortion effects to be partially compensated, but, especially at high acoustic frequencies, precise focusing is extremely important. A comprehensive review of analytical, numerical, and experimental skull-induced aberration correction techniques can be found in [28].

Moreover, the ability to predict and avoid unwanted secondary effects during the procedure, such as the development of standing pressure waves, especially in the case of low-frequency ultrasound or long sonications, or secondary hot-spots on the patient’s scalp and at bone-tissue and air-tissue interfaces, would be beneficial. Currently, this is feasible only through realistic acoustic and thermal numerical modeling of the entire procedure to derive the 3D pressure and temperature distributions.

The purpose of the current study is to introduce a novel integrated simulation framework that comprises both linear and nonlinear parallelized 3D full-wave acoustic solvers and powerful thermal solvers. These solvers have been tailored to simulate treatment setups involving heterogeneous anatomical models and are complemented by flexible geometric modeling, image segmentation, and post-processing tools that allow for rapid and precise acoustic and thermal modeling of FUS treatments. This framework was used to illustrate the feasibility of simulation-based optimization in a tcFUS therapy scenario. It was further used to investigate the impact of skull aberrations, compare different phase and amplitude correction approaches used to compensate for these aberrations, and achieve refocusing.

The simulated setup consists of an idealized model of the commercially available tcMRgFUS system ExAblate®; Neuro (InSightec, Haifa, Israel) and the detailed anatomical head model of a 34-year-old, healthy male; 22 targets, in various locations in the brain of the head model, were defined and examined. Phase and amplitude corrections for each of these targets were calculated according to four different approaches, and the impacts on the acoustic focus and the induced temperature increase were investigated. In addition, the effect of temperature-dependent perfusion on the temperature distribution was also explored.

Simulation framework

For the purposes of acoustic and thermal modeling of transcranial FUS propagation, an acoustic solver that allows for flexible 3D modeling of the simulation setup was developed and coupled to our existing thermal solver; the new solver has been integrated into our simulation platform SIM4LIFE (Zurich MedTech, Zürich, Switzerland) and optimized for anatomical model simulations, e.g., gridding and voxeling. In addition, an integrated medical image segmentation platform, iSEG (Zurich MedTech, Zürich, Switzerland), enabled the generation of patient-specific models and has been used to create high-quality reference anatomical surface models of healthy volunteers [29, 30]. The entire framework is supported by a versatile Python scripting interface, which allows the entire modeling, simulation, and post-processing procedure to be automated, and Visualization Toolkit (Kitware, New York, USA) based post-processing capabilities built into the platform to permit flexible visualization and evaluation of the results.

Acoustic solver

An acoustic solver based on the finite-difference time-domain (FDTD) method [3133] was developed to perform full-wave simulations.

Both linear and nonlinear variants, based on modified versions of the linear acoustic pressure wave equation (LAPWE) [34, 35] and the Westervelt-Lighthill equation (WLE) [36, 37], respectively, were developed. Both equations were extended with a density variation term to account for the change in density between neighboring voxels in the domain and to correctly capture the reflections at air-tissue and bone-tissue interfaces; FDTD stencils were tailored to nonuniform rectilinear grids to allow for flexible gridding of the computational domain. The LAPWE partial differential equation (PDE) is:

$$ \rho\nabla\frac{1}{\rho}\nabla p-\frac{1}{c^{2}}\frac{\partial^{2} p}{\partial t^{2}} - \frac{\stackrel{\sim}{a}}{c^{2}}\frac{\partial p}{\partial t}=0 $$

and the WLE PDE is:

$$ \rho\nabla\frac{1}{\rho}\nabla p-\frac{1}{c^{2}} \frac{\partial^{2} p}{\partial t^{2}}+\frac{\delta}{c^{4}} \frac{\partial^{3} p}{\partial t^{3}}+\frac{\beta}{2\rho c^{4}}\frac{\partial^{2} p^{2}}{\partial t^{2}}=0 $$

where ρ is the material density in kg m −3, p is the acoustic pressure in Pa, c is the speed of sound in m s −1, a is calculated as \(2a\sqrt {\frac {a^{2} c^{4}}{4\pi ^{2} f^{2}}+c^{2}}\) with a being the material attenuation coefficient in Np/m and f being the wave frequency in Hz. Parameter δ is the diffusivity of sound for a thermoviscous fluid and β is the nonlinearity coefficient. The nonlinearity coefficient β is related to the nonlinearity parameter of the fluid B/A, which is a pure number, by β=1+B/2A, while the diffusivity δ can be expressed as a function of the absorption coefficient of the fluid α, with \(\delta =\frac {2\alpha c^{3}}{4\pi ^{2} f^{2}}\) [36, 37].

Numerical truncation of the computational domain was achieved by fitting the solvers with heterogeneous perfectly matched layer (PML) boundaries [38] that allow the simulated domain to be restricted to the volume of interest without introducing reflections. The formulation was based on the stretched-coordinate approach first proposed in [38] and subsequently derived for scalar wave equations in [39].

The two solver variants were parallelized for both multi-core systems using OpenMP and GPU devices using NVIDIA’s CUDA, resulting in a speedup factor of up to 45 times in the case of the latter, in order to allow the simulation of large computational domains in viable time frames.

In order to ensure the sound operation of the acoustic solvers developed within the presented framework, analytical, numerical, and experimental validation in water-tank setups was performed. For the sake of brevity, however, only the analytical and numerical validation of the solvers will be presented. Experimental validation of water-tank setup measurements employing focus-distorting obstacles, as well as validation against experiments involving ex vivo human calvaria, will be presented in an upcoming study.

Thermal solver

A thermal solver tailored to biomedical applications has been previously developed, validated, and integrated into our simulation platform SIM4LIFE. The solver is based on a finite-difference implementation, with conformal corrections, of Pennes’ bioheat equation (BHE) [40]:

$$ \rho C\frac{\partial T}{\partial t} = \nabla \cdot \left(k\nabla T\right) + \rho Q + \rho S - \rho_{b} c_{b} \rho \omega \left(T - T_{b}\right) $$

where ρ is the material density, C is the specific heat capacity, T is the tissue temperature, k is the thermal conductivity, Q is the metabolic heat generation rate, ω is the perfusion rate, and ρ b , c b , and T b are the density, specific heat capacity, and temperature, respectively, of the blood. The term S denotes the time-averaged rate of heat generation by relaxation absorption in a tissue of a continuous sound field. The formula for this term is [41]:

$$ S = a\frac{p^{2}}{\rho c} $$

By introducing this term into Eq. 3, it is possible to couple the two solvers and calculate the temperature induced in the tissue due to exposure to the acoustic fields. The thermal solver has the ability to account for thermoregulation and vascular shutdown and has been augmented with a wide range of perfusion models, including the discrete vasculature (DIVA) [42] and Weinbaum-Jiji (WJ) [43] models, support for MRI perfusion maps, etc., as well as Arrhenius tissue damage and thermal dose models [44], thus permitting realistic modeling and assessment of thermal effects in the body. In addition, the thermal solver allows boundary conditions to be applied to selected interfaces between different tissues or regions. Three types of boundary conditions can be employed:

  • Dirichlet: A fixed temperature enforced at the interface (T=T boundary).

  • Neumann: A fixed thermal energy flux enforced at the interface \(\left (k\frac {dT}{dn} = F_{\text {boundary}}\right)\).

  • Mixed/Convenctive: Energy flux that depends on the local surface temperature and equilibrates it to the specified environment temperature T outside based on a heat transfer coefficient h in Wm −2 K −1, while in addition, a fixed heat flux can also be

    F boundary added \(\left (k\frac {dT}{dn} +h\left (T-T_{\text {outside}}\right) = F_{\text {boundary}}\right)\).

Further details on the implementation and validation of the thermal solver can be found in [45, 46].

Acoustic solver validation

Numerical validation

Numerical validation of the lossless and lossy LAPWE solvers was performed in simplified setups against the freely available software FOCUS [4749]. FOCUS employs an entirely different model of acoustic wave propagation, namely the fast near-field method (FNM) [4952], which is based on numerical approximations of analytical solutions and which has been extensively validated by its authors against the well-established software Field II [48].

Simplified transducers were initially simulated with the FNM module available in FOCUS. Subsequently, simulations of the same transducers were performed with the presented framework using the LAPWE model FDTD implementation. Three distinct transducers sonicating at 500kHz were simulated. These were a circular transducer with a radius of 10 mm, a rectangular transducer with dimensions of 20×20 mm, and a planar ring transducer with inner and outer radii of 7.5 and 10.0 mm, respectively. The transducers were embedded in an infinite homogeneous medium with acoustic properties akin to those of water, i.e., speed of sound c of 1500 m s −3 and density ρ of 1000 kg m −3. In order to separately validate the lossless and the lossy LAPWE solvers, two series of validations were performed with an attenuation coefficient α of 0.0 and 5.756 Np m −1 [53], respectively.

The truncated computational domains had dimensions of 40×40×90 mm and were discretized with a 0.25-mm grid step, which amounts to λ/12, where λ is the acoustic wavelength for the given frequency and medium. In the case of the LAPWE solver, these domains were truncated with 16 layers of PML in order to inhibit the manifestation of spurious reflections at the domain boundaries.

While the FNM model directly calculates the steady-state pressure distribution, the acoustic solvers presented in this work are explicit, time-domain solvers using the FDTD method. Thus, to assess the necessary number of simulated periods required to achieve steady-state, multiple simulations over 50–90 periods were performed. It was ascertained that 60 periods were sufficient to achieve steady-state, as longer durations resulted in less than 0.1 % difference in terms of absolute pressure.

The absolute pressure was plotted along the axis of propagation for both FOCUS and LAPWE in all six cases (see Fig. 1). Good agreement, with errors on the order of 1.9–3.5 % (see Table 1), was achieved between FOCUS and the LAPWE solver for all six comparison cases with only minor differences in the far-field regions, which are attributed to the cumulative phase dispersion errors inherent to the FDTD method [54]. In order to further quantify the agreement between the two pressure plots, the normalized standard deviation was calculated as the ratio of the 2 norm of the pressure difference between the two lines and the 2 norm of the FOCUS results:

$$ \begin{aligned} \text{Normalized standard deviation \%}\\ = 100\frac{{\ell}^{2}\left(\left|p^{\text{LAPWE}}\right|-\left|p^{\text{FOCUS}}\right|\right)}{{\ell}^{2}\left(\left|p^{\text{FOCUS}}\right|\right)} \end{aligned} $$
Fig. 1
figure 1

Validation of LAPWE against FOCUS. Absolute pressure line plots along the axis of propagation as calculated by both FOCUS and the LAPWE solver. The pressure comparisons within a lossless and a lossy medium can be seen for the cases of a circular (a), a rectangular (b), and a planar ring (c) transducers. Excellent agreement can be seen for all cases

Table 1 The normalized standard deviation values between the FOCUS and the LAPWE solver presented in this work. Good agreement between the two solvers can be seen for all transducer types and media examined

with 2 norm being the square root of the sum of squares of all absolute pressure values for the respective line. These results can be seen in Table 1.

Analytical validation

In order to validate the newly introduced density variation terms (see “Acoustic solver” section), which are not available in other numerical implementations, analytical validation based on the calculation of the reflection R and transmission T coefficients [55] at the interface between media with varying characteristic acoustic impedances Z was also performed. A 1-MHz acoustic plane wave was propagated through media with varying speeds of sound c and densities ρ at different incidence angles θ i . The reflection R and transmission T coefficients were calculated both analytically, through the characteristic acoustic impedances of the involved media and, numerically, through the pressure amplitudes of the incident, reflected, and transmitted waves at the interfaces of the different media based on the formulas below:

$$\begin{array}{@{}rcl@{}} R = \frac{p_{r}}{p_{i}} = \frac{\left(Z_{2}\cos\theta_{i} - Z_{1}\cos\theta_{t}\right)}{\left(Z_{2}\cos\theta_{i} + Z_{1}\cos\theta_{t}\right)} &&\\ T = \frac{p_{t}}{p_{i}} = \frac{\left(2Z_{2}\cos\theta_{i}\right)}{\left(Z_{2}\cos\theta_{i} + Z_{1}\cos\theta_{t}\right)} \end{array} $$

where p i , p r , and p t are the pressure amplitudes of the incident, reflected, and transmitted waves, respectively. The transmission angle θ t of a plane wave propagating from a medium A to a medium B with incidence angle θ i was calculated based on Snell’s law:

$$ \frac{\sin\theta_{i}}{c_{A}} = \frac{\sin\theta_{t}}{c_{B}} $$

where c A and c B are the speeds of sound in media A and B, respectively.

Nine such simulations were performed where the incidence of a plane wave from a medium with a speed of sound c of 1500 m s −1 and density ρ of 1000 kg m −3 on three different media at incidence angles of 0°, 15°, and 30° was modeled. Both the analytically (reference) and the numerically determined coefficients can be seen in Table 2. They are in excellent agreement within the numerical accuracy, thus verifying the correct implementation of the density variation and acoustic propagation term in the LAPWE solver.

Table 2 Analytical and numerical calculations of the reflection R and transmission T coefficients for the nine different cases. In the analytical case, the coefficients are calculated based on the acoustic impedance of the two media. The numerical calculations are based on the ratios between the pressure amplitudes of the reflected and transmitted waves to the amplitude of the incident wave

Transcranial FUS simulations

The simulation framework described above was used to simulate a clinically relevant tcFUS scenario, in which a detailed anatomical head model segmented from MRI data was placed within the transducer cavity of the ExAblate®; Neuro system. Four approaches, ranging from (semi-)analytical to simulation-based, were employed to calculate phase—and optionally amplitude—corrections of the skull-induced aberrations. The quality of the corrections was assessed both acoustically and thermally, and the different approaches were compared.

Simulation setup

Head model

The head model used in these simulations is an improved version of “Duke” of the “Virtual Population” collection of surface-based reference anatomical models segmented from MRI data of healthy volunteers [29, 30]. This improved version of the model, shown in Figs. 2 and 3, was segmented with iSEG (Zurich MedTech, Switzerland) from a re-sampled MRI dataset with a resolution of 0.5×0.5×0.5 mm, yielding over 50 individual tissues and anatomical structures in the head region alone. Additional details about the Virtual Population project can be found in [30].

Fig. 2
figure 2

“Duke” anatomical head model and ExAblate®; 4000 transducer array model. The “Duke” anatomical head model and the ExAblate®; 4000 transducer array were used in this study. The positioning of the head model within the transducer as well as the location of the geometric focus in the right thalamic ventral intermediate (VIM) nucleus can be seen

Fig. 3
figure 3

Sonication targets in the head model. The targets defined in the head model: a “cortex targets,” which all lie above the transducer’s focal plane, and b “structure targets,” where the red dashed line shows the transducer’s focal plane in relation to the targets, which heavily influences the achievable focusing. The initials L, R, A, M, and P stand for left, right, anterior, medial, and posterior, respectively, and the numbers in parentheses are the target distances in millimeter from the geometric focus of the transducer

Applicator model—ExAblate®; Neuro

The applicator model used in these simulations was based on the ExAblate®; Neuro (InSightec, Haifa, Israel), a tcMRgFUS system widely used in the FUS community, which is under evaluation for clinical safety and efficacy in functional neurosurgery, tumor ablation, and targeted drug delivery [1, 35, 8].

The ExAblate®; Neuro consists of a 30 cm diameter hemispherical phased-array transducer with 1024 elements operating at either 230 or 650 kHz. This device is coupled with a 1024-channel amplifier, which allows phase and amplitude control of each individual transducer element in the phased array.

A model applicator (operating at a frequency of 230 kHz), shown in Fig. 2, was generated to mimic the actual applicator, with the same number of transducer elements in similar groupings. All elements were modeled identically, with a surface area of 1 cm 2.

Patient positioning and target definition

To replicate clinically relevant setups, the head model was placed so that the right thalamic ventral intermediate (VIM) nucleus was located at the geometric center of the transducer, i.e., where the acoustic waves of all the elements, driven in phase with no source of phase aberrations—e.g., the presence of the skull—would converge according to the arrangement of the elements in the array, as shown in Fig. 2.

To evaluate the different refocusing approaches that are described in the next section and the ability of each to focus the acoustic waves at a given location further away from the transducer’s geometric focus than is currently considered feasible, two types of targets, “structure targets” and “cortex targets”, were defined in the anatomical model. Structure targets, regions of the human brain that have been associated with different neuropathic conditions, are very desirable targets for tcFUS neurosurgery; 14 such targets were defined and are shown in Fig. 3 b. In addition, eight cortex targets in different cortices of the brain were defined, as shown in Fig. 3 a, which, because of close proximity to the skull, are thought by the FUS community to be untreatable.

US simulations

Target focusing and aberration correction

To correct both amplitudes and phases for every element of the array, the LAPWE-based linear acoustic solver (see Eq. 1) was applied to each of the defined targets according to the following procedure:

  • A point source, driven with an arbitrary amplitude, was placed at the intended target location.

  • An inverse-propagation acoustic simulation of the entire head and applicator in which the waves from the point source were allowed to propagate for 1.0 ms, i.e., 230 periods at 230 kHz resulting in a propagation distance of ca. 1.5 m, was performed [56].

  • During the simulation, the transducer elements were used as receivers to record the complex pressure values, i.e., amplitude and phase, at the surface center of every element.

  • Four distinct focusing strategies were applied and evaluated as follows:

    1. 1.

      Distance-based phase corrections (DPC): Analytical phase corrections for every element were calculated based on the distance between each element’s surface center and the desired target, assuming that the transducer is in a homogeneous water medium, thus without wave distortion taken into account. The calculated distance-based phase corrections are applied to each element, the amplitudes of which are fixed to the appropriate pressure level for a given acoustic input power (see “Pressure levels” below).

    2. 2.

      Ray-tracing-based phase corrections (RTPC): As an extension of the DPC approach, a ray-tracing algorithm that takes the skull properties into account and allows calculation of improved effective distance-based phase corrections was devised. The algorithm calculates the skull entry/exit, i.e., intersection, points of the rays between each element’s surface center and the desired target and uses that information to calculate the thickness of the skull through which the waves from a given element propagate on the way to the target. Like in the DPC approach, the improved phase corrections and fixed pressure amplitudes are applied to each element.

    3. 3.

      Simulation-based phase corrections (SPC): The phases of the pressure phasors recorded during the inverse-propagation simulation are conjugated, and the amplitudes are fixed to the appropriate constant pressure level for a given acoustic input power.

    4. 4.

      Simulation-based phase and amplitude corrections (SPAC): The phases of the pressure phasors recorded during the inverse-propagation simulation are conjugated, and the recorded amplitudes, normalized to a given acoustic input power, are used (see “Amplitude normalization” below).

  • Subsequently, for each of the four focusing strategies, the calculated phase—and in the case of the SPAC approach, amplitude—corrections were applied to their respective elements and a forward-propagation acoustic simulation was performed to investigate the acoustic and thermal impact of those corrections.

It should be mentioned that the SPC/SPAC approaches are based on the “Virtual Source” time-reversal approach more information on which can be found in [28, 5759].

Gridding and voxeling

The simulation gridding was set up with at least 10 cells per minimum wavelength, which resulted in approximately 80×106 voxels for each simulation.

Pressure levels

The definition of the actual pressure amplitude level(s) for the array elements for a given acoustic input power is pivotal to obtaining realistic pressure values at the focus location and to calculating the induced temperature increase. The following formula was used to calculate the average element surface pressure for a given acoustic input power:

$$ p_{\text{element}} = \sqrt{\frac{P_{\text{total}} \cdot Z_{\text{water}}}{N_{\text{elements}} \cdot A_{\text{element}}}} $$

where p element is the pressure at the surface of an array element in Pa, P total is the total acoustic input power in W, Z water is the characteristic impedance of water in Rayl (or kg m −2 s −1 in SI units), N elements is the number of active elements in the array, and A element is the area of the element surface in m 2.

Assuming that water has a sound propagation speed of 1482.3 m s −1 [25] and a density of 1000 kg m −3, the resulting characteristic impedance Z water is 1.4823×106 kg m −1 s −1. Therefore, when all 1024 elements are active, sonicate at the same pressure, and have a surface area of A element= 1.0 cm 2, the pressure at the surface of every element for 1000 W acoustic input power is p element,1000 W≈ 120.31 kPa.

Amplitude normalization

In the SPC/SPAC approaches discussed prior, the complex pressure waves emanating from a point source with an arbitrary amplitude are captured, conjugated, and re-emitted in a subsequent simulation to achieve refocusing. While uniform pressure amplitudes are used across all elements in the SPC approach, in the case of the SPAC approach, the pressure amplitudes of these waves must be normalized. The factor f used for normalization of the captured complex pressure values is defined as:

$$ f= \sqrt{\frac{N_{\text{elements}} \cdot p^{2}_{\text{element}}}{ \sum_{i~=~1}^{N_{\text{elements}}} p^{2}_{\text{captured}_{i}}}} $$

where N elements is the number of active elements in the array, p element is the pressure amplitude calculated in the previous section, and \(p_{\text {captured}_{i}}\phantom {\dot {i}\!}\) is the pressure recorded at the ith element of the array during the inverse-propagation simulation.

Acoustic tissue properties

Although most artificial materials have been acoustically characterized for use in transducer manufacture or nondestructive testing (NDT), an extensive literature search revealed no comprehensive studies on the acoustic properties of human tissue, apart from the properties of the bone which were first investigated by Fry [24]. While promising projects dealing with the measurement of such properties are underway [60], the only existing literature is generally decades old and comprises empirically measured properties with large discrepancies between different datasets. Such sources can be found in [25, 55, 6165].

Consequently, most numerical studies on ultrasonic wave propagation involve simple modeling of a skull surrounded by water, with the argument that the ultrasonic tissue properties show little heterogeneity between different soft tissues.

In this study, however, accurate modeling of the entire head anatomy including all segmented tissues was necessary to allow for acoustic and, especially, thermal modeling of the procedure. To that end, the properties reported in [61] were used, while the grouping of soft tissue types based on physical composition was required to account for the lack of an extensive property database. The acoustic properties used in these simulations are summarized in Table 3. The material density ρ values were based on the IT’IS Foundation Tissue Properties Database [66]. Given that the anatomical model used was based on MRI data, it was not possible to acquire voxel-specific bone properties through Hounsfield units; thus, constant acoustic properties were assigned per tissue.

Table 3 The acoustic tissue properties used in this study where c is the speed of sound and a is the material attenuation coefficient [25, 61]

Thermal simulations

Following the acoustic simulations, the deposited acoustic energy was calculated for every voxel of the computational domain, using Eq. 4, and used as input in the thermal solver (see Eq. 3).

To realistically model the entire treatment setup, convective thermal boundary conditions were applied at the interfaces between tissues and the water-bolus surrounding the head as well as at air-tissue interfaces, both for the internal air in head cavities and the air surrounding the head. The water temperature was fixed to 16 °C [3, 4], and a heat-transfer coefficient h of 70 W m −2 K −1 was applied [45]. Similarly, both internal and external air were fixed to 25 °C with a heat-transfer coefficient h of 6 W m −2 K −1.

To properly account for the cooling effect of the water-bolus on the scalp, thermal simulations were performed for 30 min in the absence of sonication to allow the different tissues to reach thermal equilibrium, generally achieved after ca. 10 min. This was then followed by 20 s of sonication to calculate the temperature increase induced by the deposited acoustic energy during treatment. During these simulations, the effects of perfusion were taken into account, and the thermal properties for all tissues were based on [66].

These thermal simulations were repeated for every target in the head and aberration correction approach in order to calculate the induced temperature increase.

Vascular shutdown and temperature-dependent tissue perfusion

The above thermal simulations were performed with the assumption that thermal tissue properties are not temperature dependent. However, it is known that, during thermal ablation, the high temperatures induce vascular shutdown, thus eliminating perfusion in those locations and causing the temperature to increase more rapidly [67, 68]. To assess the importance of this effect, additional thermal simulations were performed for a few select targets where perfusion was assumed to start decreasing linearly when the tissue temperature exceeded 50 °C and cease entirely above 51 °C. Vascular dilation, which would cause an exponential increase in blood perfusion as a function of temperature, is not taken into account, as this effect manifests itself only after several minutes of exposure to increased local temperature [6971].


The four aberration correction approaches described under the “Thermal simulations” section were applied to all targets shown in Fig. 3. The resulting pressure distributions were used to calculate the temperature increase as described under the “Thermal simulations” section.

Due to the target location diversity, the targets were loosely categorized as cortex targets, structure targets at or above the transducer’s focal plane, and structure targets below the focal plane.

To analyze the acoustic and thermal performance of the different approaches for each target in a consistent manner, an automatized local maxima and connected-component analysis was performed on all calculated 3D pressure and temperature distributions. Firstly, the distributions were filtered, and all local pressure/temperature maxima were identified. Subsequently, the local maximum that was nearest to the intended target was detected and assumed to denote the focal region/primary lesion. It should be noted that henceforth the term “lesion” will be used to refer to the region of the thermal hot-spot, even though in some cases, the intensity of the hot-spot might not be sufficient to create an actual lesion.

Once that maximum and the peak absolute pressure or temperature increase were identified, the simulated fields were thresholded at 50 % of the peak pressure or temperature level, and the different connected components were analyzed. This yielded the full-width half-maximum (FWHM) size of the focal region or thermal lesion along the X, Y, and Z axes, the distance between that region’s center and the intended target, as well as the volume of the region, calculated as the sum of the voxel volumes belonging to the particular component. The results of the analyses are summarized in Table 4 for the acoustic pressure distributions and in Table 5 for the temperature increase distributions.

Table 4 The results, mean (standard deviation), of the connected-component analysis on the absolute acoustic pressure distributions for the four focusing approaches for the different target categories
Table 5 The results, mean (standard deviation), of the connected-component analysis on the temperature increase distributions for the four focusing approaches for the different target categories

A plot of the peak absolute pressure in these detected foci for all targets and approaches is shown in Fig. 4. The absolute pressure distribution obtained with the four approaches, in the case of the “Thalamic VIM (right)” target—which coincides with the geometric focus of the transducer—can be seen in Fig. 5. In addition, the absolute pressure distributions resulting from the use of the SPC approach in the case of four selected targets can be seen in Fig. 6.

Fig. 4
figure 4

Peak target pressure for different targets and focusing approaches. The peak absolute pressure (in MPa) achieved near each of the defined targets for all focusing approaches. It can be clearly seen that the SPC and SPAC approaches yield far stronger foci than the DPC and RTPC approaches

Fig. 5
figure 5

Pressure distribution comparison for different focusing approaches. Absolute pressure distribution in megapascal resulting from the use of the four approaches with the “Thalamic VIM (right)” target. Each distribution is plotted on the sagittal plane through the target, accompanied by the respective color map and scaled to the respective maximum absolute pressure. The distribution resulting from the DPC approach a shows a heavily distorted focal region and relatively low pressure amplitude at the target. Delineation of the focal region and amplitude improve slightly with the RTPC approach (b), while a significant improvement is seen in the case of the SPC (c) and SPAC (d) approaches

Fig. 6
figure 6

Pressure distribution comparison for four types of targets. Absolute pressure distribution in megapascal resulting from the use of the SPC approach with four selected targets. Each distribution is plotted on a plane through the target, accompanied by the respective color map and scaled to the respective maximum absolute pressure. Distributions ac are plotted on the sagittal plane, while distribution d is plotted on the coronal plane. The two structure targets above the transducer’s focal plane, “Thalamic VIM (left)” (a) and “Corpus Callosum (medial)” (b), show sharply delineated focal regions with very high pressure amplitudes. The resulting distribution for “Amygdala (left)” (c), a structure target below the focal plane, shows successful refocusing but significantly lower pressure amplitude. In the case of the cortex target “Auditory Cortex (left)” (d), a focus is visible at the target location but significant energy deposition in the patient’s scalp and skull is observed

A plot of the peak FUS-induced temperature increase achieved in each of the defined targets for all focusing approaches after 20 s of sonication is shown in Fig. 7.

Fig. 7
figure 7

Peak target temperature increase for different targets and focusing approaches. The peak FUS-induced temperature increase, not absolute temperature, achieved with each of the defined targets for all focusing approaches after 20 s of sonication. The red dashed line shows the 50 °C threshold, assumed to be the ablation threshold for a base tissue temperature of 37 °C


Acoustic distributions

As can be seen in Table 4, the acoustic pressure achieved in the focal regions nearest to the targets is lowest for the DPC approach—which neglects the impact of the skull—slightly higher for the RTPC approach (+10 % in average) and much higher for the SPC and SPAC approaches (+107 and +148 %, respectively). In terms of target location, structure targets above the transducer’s focal plane exhibited higher pressures than those below, while the achievable pressure in the cortex targets was the lowest.

The pressure level dependence on the focusing strategy can be mostly understood when the volume V of the different focal regions is considered. The phase corrections calculated through the DPC and RTPC approaches resulted in focal regions much larger, by a factor of 5.2 on average, than those obtained with simulation-based methods (see Fig. 5).

While the SPC and SPAC approaches yielded a nearly constant region size and shape for the different targets, focusing sharpness drastically decreased in the case of the DPC and RTPC approaches for structure targets below the transducer’s focal plane and even further for cortex targets, where dramatic size variations were observed. In general terms and regardless of the correction approach, focusing quality decreases as the distance between the intended target and the geometric focus of the array increases (see Fig. 6). As described by F shape, i.e., the ratio of the maximum to minimum dimensions of the focal region, in Table 4, the focal regions were more spherical for the simulation-based approaches than for the DPC and RTPC approaches, which yielded more elongated regions with substantial shape variations.

In terms of focal shift, precisions on the order of the discretization resolution were typically achieved with simulation-based approaches, as can be seen from D in Table 4, while the DPC and RTPC approaches showed an average shift of 3.7 mm, which was even more pronounced for cortex targets.

Thermal distributions

When considering the temperature increase results after 20 s of sonication (see Fig. 7 and Table 5), the DPC and RTPC approaches achieved ablative temperatures only in the case of structure targets above the transducer’s focal plane. As observed with the acoustic results, temperature rise was highest for structure targets above the focal plane.

The SPAC was clearly superior to the other approaches and achieved ablative temperatures in all investigated targets. Utilization of the SPC approach attained ablation in all structure targets above the transducer’s focal plane and the majority of targets below. In the case of cortical targets, even though focusing was achievable with this approach, the limited number of elements that could contribute to the focusing (due to the absence of a line-of-sight between many of the elements and the intended target) resulted in high energy deposition on the skull and scalp, while the temperature rise at the targets showed inverse proportionality to the distance between them and the transducer’s geometric focus. Effective thermal treatment of those targets would require for the anatomical model to be repositioned so that the desired targets would lie above the focal plane.

The thermal lesion and acoustic focus volumes exhibited similar behaviors, with the simulation-based approaches being clearly superior to the DPC and RTPC, especially for the cortex and structure targets below the transducer line where these approaches were unable to produce sharply demarcated lesions.

Skull heating

With the SPC approach, high acoustic energy deposition and subsequent thermal hotspots were observed near the skull surface for the majority of targets, which would result in significant heating of the scalp and skull. This phenomenon was partly alleviated for targets near the geometric focus with the use of RTPC, where improvements in thermal and focal gain were seen. This would suggest that, semi-analytical approaches, like the RTPC approach would be more appropriate for centralized targets, with longer sonication durations, or, for nonablative therapies.

In the case of the simulation-based approaches, these adverse effects were mostly observed for cortex targets and structure targets below the transducer level, where only a small number of elements could contribute to the focus, thus resulting in significant energy deposition on the scalp and skull bone (see Fig. 6). This trend was visible for both approaches, which leads us to conclude that these targets would benefit from further optimization of the steering parameters, e.g., deactivation of the nearby elements.

However, the pressure wave equation employed in this work does not capture shear-wave-related effects expected to occur near the skull, which is required to reliably predict the related secondary hot-spots. Therefore, further investigations are needed to extend the treatment envelope in clinical practice.

Impact of vascular shutdown

As discussed under the “Vascular shutdown and temperature-dependent tissue perfusion” section, additional thermal simulations were performed where the impact of vascular shutdown was considered. The phase corrections of the SPC approach were used in vascular shutdown simulations for the “Thalamic VIM (left)” and “Thalamic VIM (right)” targets. The impact of vascular shutdown was visible but minimal. Thermal simulations where the perfusion was assumed to decrease linearly after the tissue temperature exceeded 50 °C and cease entirely after 51 °C showed a steeper temperature increase, but after 20 s of sonication, only 1 °C of the total additional temperature increase was observed in tissues where vascular shutdown occurred. Hence, this effect was considered to have negligible impact on the formation of thermal lesions.

Impact of skull heterogeneity

The acoustic simulations performed in this study approximated the individual tissues as homogeneous, i.e., a single set of acoustic properties was assigned per tissue (see “Acoustic tissue properties” section). However, the particularly heterogeneous nature of the skull could have an impact on the degree of phase aberrations induced and on the quality of the subsequent compensation. To investigate this impact, additional acoustic simulations were performed utilizing the DPC and SPC approaches described in the “Thermal simulations” section in a setup comprising the applicator model discussed under the “Applicator model—ExAblate®; Neuro” section and two models of a human skull based on a CT dataset. For these simulations, a sonication target coinciding with the geometrical focus of the applicator was defined.

The CT dataset of an ex vivo scan with a resolution of 0.48×0.48×2.5 mm was segmented based on Hounsfield units (HU) thresholding where HU values greater than 700 were considered to signify bone structures [72]. The HU values were then used to calculate the porosity and subsequently the density and speed-of-sound for each voxel in the skull, using the relationships from [73].

One of the generated models was fully homogeneous whereas the other one was partially heterogeneous and consisted 20 HU bins, within which voxels with similar HU values were grouped into a single region. This binning led to a maximum deviation of 2 % from the original HU values and was necessary for the current solver implementation which does not support individual tissue properties per voxel.

The results obtained are within range of the values for structure targets (ST) above the focal plane in Table 4. Comparing the pressure distribution in the homogeneous and the inhomogeneous skull model shows that the SPC approach produces identical shifts (ca. 0.5 mm) and focal region volumes (ca. 50.0 mm 3), while considering heterogeneity reduces the peak pressure by 25 %. A similar reduction of 35 % in peak absolute pressure amplitude was also observed with the DPC approach. In addition, the DPC approach in the case of the heterogeneous skull exhibited both more pronounced focal region shifts (factor of 3) and defocusing (33 % volume decrease) when compared to the homogeneous case.

These results from a single case indicate that the prediction of the achievable SPC-based focusing (shift and volume) is unaffected by skull inhomogeneity, while prediction of aberration effects could be overestimated when employing the DPC approach and using a homogeneous model. The assumption of image data-based skull property distributions results in lowered predicted peak pressure, possibly due to the increased reflections and scattering at the skull.


A well-known limitation of tcFUS therapy is the skull-induced aberrations, which can induce focal shift and distortion as well as significant energy deposition on the patient’s skull and scalp, resulting in a significant decrease in the treatment’s focal and thermal gain. A numerical feasibility study was performed here to investigate the efficacy of four compensation techniques, ranging from (semi-)analytical to simulation-based, that aim to provide phase—and optionally amplitude—corrections to achieve refocusing, counter the aforementioned effects, and increase the treatment envelope of tcFUS therapy. To that end, an extensive, newly developed and validated simulation framework that allows for 3D full-wave, linear and nonlinear acoustic, and thermal simulations in large and complex clinically relevant setups, was utilized to perform simulations of a detailed anatomical head model sonicated with a model of the ExAblate®; Neuro applicator. The acoustic and thermal results of these correction approaches were ascertained for 22 distinct targets in various locations of the brain. Good overall agreement in focal pressure and temperature can be seen between this and similar, recently published, numerical studies [74].

Evaluation of the acoustic pressure, location, and size of the focal regions as well as the FUS-induced temperature increase and lesion volume/size suggests that simulation-based approaches provide far superior corrections than the analytical and (semi-)analytical ones. While the latter could be employed in thermal interventions for targets in the vicinity of the transducer’s geometric focus, their efficiency decreased dramatically when targeting more remote brain regions. Simulation-based approaches, on the other hand, appear to yield sharply demarcated focal regions and lesions at targets multiple centimeter from the geometric focus, as well as the ability to better utilize and control large transducer arrays. Employment of modern but affordable computer hardware, combined with state-of-the-art high-performance computing techniques such as those described for the numerical framework presented here, enables realistic acoustic and thermal simulations in complicated setups to be performed within minutes. Due to their increased refocusing efficiency and their ability to predict the acoustic and thermal effects of FUS therapies, if extensively validated, simulation-based correction approaches may eventually replace their analytical counterparts.

In future work, additional experimental validation of the presented aberration correction approaches should be carried out against measurements of ex vivo human calvaria. Furthermore, the impact of bone density heterogeneity on the quality of aberration correction should be further investigated through employing MR and CT image data, the impact of shear-wave-induced effects must be considered to reliably predict and exclude secondary hot-spots, and the analysis needs to be extended to higher frequencies, e.g., the 650-kHz system. It is necessary to develop a comprehensive uncertainty assessment and validation procedure to enable these techniques to be used in treatment planning.


  1. McDannold N, Clement G, Black P, Jolesz F, Hynynen K. Transcranial MRI-guided focused ultrasound surgery of brain tumors: initial findings in three patients. Neurosurgery. 2010; 66(2):323.

    Article  PubMed Central  PubMed  Google Scholar 

  2. Coluccia D, Fandino J, Schwyzer L, O’Gorman R, Remonda L, Anon J, et al.First noninvasive thermal ablation of a brain tumor with MR-guided focused ultrasound. J Ther Ultrasound. 2014; 2(1):17.

    Article  PubMed Central  PubMed  Google Scholar 

  3. Martin E, Jeanmonod D, Morel A, Zadicario E, Werner B. High-intensity focused ultrasound for noninvasive functional neurosurgery. Ann Neurol. 2009; 66(6):858–61.

    Article  PubMed  Google Scholar 

  4. Jeanmonod D, Werner B, Morel A, Michels L, Zadicario E, Schiff G, et al.Transcranial magnetic resonance imaging–guided focused ultrasound: noninvasive central lateral thalamotomy for chronic neuropathic pain. Neurosurg Focus. 2012; 32(1):1.

    Article  Google Scholar 

  5. Elias WJ, Huss D, Voss T, Loomba J, Khaled M, Zadicario E, et al.A pilot study of focused ultrasound thalamotomy for essential tremor. N Engl J Med. 2013; 369(7):640–8.

    Article  CAS  PubMed  Google Scholar 

  6. Moser D, Zadicario E, Schiff G, Jeanmonod D. Measurement of targeting accuracy in focused ultrasound functional neurosurgery: technical note. Neurosurg Focus. 2012; 32(1):2.

    Article  Google Scholar 

  7. Lipsman N, Schwartz ML, Huang Y, Lee L, Sankar T, Chapman M, et al.MR-guided focused ultrasound thalamotomy for essential tremor: a proof-of-concept study. Lancet Neurol. 2013; 12(5):462–8.

    Article  PubMed  Google Scholar 

  8. Monteith S, Sheehan J, Medel R, Wintermark M, Eames M, Snell J, et al.Potential intracranial applications of magnetic resonance-guided focused ultrasound surgery: a review. J Neurosurg. 2013; 118(2):215–21.

    Article  PubMed  Google Scholar 

  9. Alexandrov AV, Molina CA, Grotta JC, Garami Z, Ford SR, Alvarez-Sabin J, et al. Ultrasound-enhanced systemic thrombolysis for acute ischemic stroke. N Engl J Med. 2004; 351(21):2170–8.

    Article  CAS  PubMed  Google Scholar 

  10. Daffertshofer M, Gass A, Ringleb P, Sitzer M, Sliwka U, Els T, et al.Transcranial low-frequency ultrasound-mediated thrombolysis in brain ischemia increased risk of hemorrhage with combined ultrasound and tissue plasminogen activator: results of a phase II clinical trial. Stroke. 2005; 36(7):1441–6.

    Article  PubMed  Google Scholar 

  11. Selvaraj P, Okita K, Matsumoto Y, Voie A, Hoelscher T, Weiss H, et al.Effects of physical properties of the skull on high intensity focused ultrasound for transcranial sonothrombolysis. J Acoustical Soc Am. 2011; 130(4):2538.

    Article  Google Scholar 

  12. Weiss HL, Ahadi G, Hoelscher T, Szeri AJ. Cavitation damage during sonothrombolysis using high intensity focused ultrasound.J Acoustical Soc Am. 2011; 129(4):2576.

    Article  Google Scholar 

  13. Monteith SJ, Harnof S, Medel R, Popp B, Wintermark M, Lopes MBS, et al. Minimally invasive treatment of intracerebral hemorrhage with magnetic resonance-guided focused ultrasound: laboratory investigation. J Neurosurg. 2013; 118(5):1035–45.

    Article  PubMed  Google Scholar 

  14. Jordão JF, Ayala-Grosso CA, Markham K, Huang Y, Chopra R, McLaurin J, et al.Antibodies targeted to the brain with image-guided focused ultrasound reduces amyloid- β plaque load in the TgCRND8 mouse model of alzheimer’s disease. PLoS One. 2010; 5(5):10549.

    Article  Google Scholar 

  15. Hynynen K, McDannold N, Vykhodtseva N, Raymond S, Weissleder R, Jolesz FA, et al.Focal disruption of the blood-brain barrier due to 260-kHz ultrasound bursts: a method for molecular imaging and targeted drug delivery. J Neurosurg. 2006; 105(3):445–54.

    Article  CAS  PubMed  Google Scholar 

  16. Reinhard M, Hetzel A, Krüger S, Kretzer S, Talazko J, Ziyeh S, et al.Blood-brain barrier disruption by low-frequency ultrasound. Stroke. 2006; 37(6):1546–8.

    Article  PubMed  Google Scholar 

  17. Wagner T, Valero-Cabre A, Pascual-Leone A. Noninvasive human brain stimulation. Annu Rev Biomed Eng. 2007; 9:527–65.

    Article  CAS  PubMed  Google Scholar 

  18. Yoo SS, Bystritsky A, Lee JH, Zhang Y, Fischer K, Min BK, et al.Focused ultrasound modulates region-specific brain activity. Neuroimage. 2011; 56(3):1267–75.

    Article  PubMed Central  PubMed  Google Scholar 

  19. Bystritsky A, Korb AS, Douglas PK, Cohen MS, Melega WP, Mulgaonkar AP, et al.A review of low-intensity focused ultrasound pulsation. Brain Stimul. 2011; 4(3):125–36.

    Article  PubMed  Google Scholar 

  20. Tufail Y, Yoshihiro A, Pati S, Li MM, Tyler WJ. Ultrasonic neuromodulation by brain stimulation with transcranial ultrasound. Nat Protoc. 2011; 6(9):1453–70.

    Article  CAS  PubMed  Google Scholar 

  21. Legon W, Sato TF, Opitz A, Mueller J, Barbour A, Williams A, et al.Transcranial focused ultrasound modulates the activity of primary somatosensory cortex in humans. Nat Neurosci. 2014; 17(2):322–9.

    Article  CAS  PubMed  Google Scholar 

  22. Hynynen K. Mri-guided focused ultrasound treatments. Ultrasonics. 2010; 50(2):221–9.

    Article  PubMed  Google Scholar 

  23. Hayat A. Tumors of the Central Nervous System, V.3, Pt.1. New York: Springer; 2011.

    Book  Google Scholar 

  24. Fry F, Barger J. Acoustical properties of the human skull. J Acoust Soc Am. 1978; 63:1576.

    Article  CAS  PubMed  Google Scholar 

  25. Szabó TL. Diagnostic ultrasound imaging: inside out. Waltham: Academic Press; 2004.

    Google Scholar 

  26. 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. 2011; 56(1):219.

    Article  PubMed Central  PubMed  Google Scholar 

  27. Clement G, Hynynen K. A non-invasive method for focusing ultrasound through the human skull. Phys Med Biol. 2002; 47(8):1219.

    Article  CAS  PubMed  Google Scholar 

  28. Kyriakou A, Neufeld E, Werner B, Paulides MM, Szekely G, Kuster N. A review of numerical and experimental compensation techniques for skull-induced phase aberrations in transcranial focused ultrasound. Int J Hyperthermia. 2014; 30(1):36–46.

    Article  PubMed  Google Scholar 

  29. Christ A, Kainz W, Hahn EG, Honegger K, Zefferer M, Neufeld E, et al. The virtual family—development of surface-based anatomical models of two adults and two children for dosimetric simulations. Phys Med Biol. 2010; 55(2):23.

    Article  Google Scholar 

  30. Gosselin MC, Neufeld E, Moser H, Huber E, Farcito S, Gerber L, et al. Development of a new generation of high-resolution anatomical models for medical device evaluation: the Virtual Population 3.0. Phys Med Biol. 2014; 59(18):5287.

    Article  PubMed  Google Scholar 

  31. Hallaj IM, Cleveland RO. Fdtd simulation of finite-amplitude pressure and temperature fields for biomedical ultrasound. J Acoust Soc Am. 1999; 105(5):7–12.

    Article  Google Scholar 

  32. Ginter S, Liebler M, Steiger E, Dreyer T, Riedlinger RE. Full-wave modeling of therapeutic ultrasound: nonlinear ultrasound propagation in ideal fluids. J Acoust Soc Am. 2002; 111:2049.

    Article  PubMed  Google Scholar 

  33. Pinton GF, Trahey GE. A comparison of time-domain solutions for the full-wave equation and the parabolic wave equation for a diagnostic ultrasound transducer. IEEE Trans Ultrason, Ferroelectr, Freq Control. 2008; 55(3):730–3.

    Article  Google Scholar 

  34. Mast TD, Hinkelman LM, Metlay LA, Orr MJ, Waag RC. Simulation of ultrasonic pulse propagation, distortion, and attenuation in the human chest wall. J Acoust Soc Am. 1999; 106:3665.

    Article  CAS  PubMed  Google Scholar 

  35. Liebler M, Ginter S, Dreyer T, Riedlinger RE. Full wave modeling of therapeutic ultrasound: efficient time-domain implementation of the frequency power-law attenuation. J Acoust Soc Am. 2004; 116:2742.

    Article  PubMed  Google Scholar 

  36. Westervelt PJ. Parametric acoustic array. J Acoust Soc Am. 1963; 35:535.

    Article  Google Scholar 

  37. 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–67.

    Article  PubMed  Google Scholar 

  38. Chew WC, Weedon WH. A 3D perfectly matched medium from modified Maxwell’s equations with stretched coordinates. Microw Opt Technol Lett. 1994; 7(13):599–604.

    Article  Google Scholar 

  39. Zhou D, Huang W, Xu C, Fang D, Chen B. The perfectly matched layer boundary condition for scalar finite-difference time-domain method. IEEE Photon Technol Lett. 2001; 13(5):454–6.

    Article  Google Scholar 

  40. Pennes HH. Analysis of tissue and arterial blood temperatures in the resting human forearm. J Appl Physiol. 1948; 1(2):93–122.

    CAS  PubMed  Google Scholar 

  41. Nyborg WL. Heat generation by ultrasound in a relaxing medium. J Acoust Soc Am. 1981; 70:310.

    Article  Google Scholar 

  42. Kotte A, Van Leeuwen G, Lagendijk J. Modelling the thermal impact of a discrete vessel tree. Phys Med Biol. 1999; 44(1):57.

    Article  CAS  PubMed  Google Scholar 

  43. Song WJ, Weinbaum S, Jiji LM. A theoretical model for peripheral tissue heat transfer using the bioheat equation of Weinbaum and Jiji. J Biomech Eng. 1987; 109(1):72–8.

    Article  CAS  PubMed  Google Scholar 

  44. Sapareto SA, Dewey WC. Thermal dose determination in cancer therapy. Int J Radiat Oncol. 1984; 10(6):787–800.

    Article  CAS  Google Scholar 

  45. Neufeld E. High resolution hyperthermia treatment planning. Konstanz: Hartung-Gorre Verlag; 2008.

    Google Scholar 

  46. Neufeld E, Kühn S, Szekely G, Kuster N. Measurement, simulation and uncertainty assessment of implant heating during MRI. Phys Med Biol. 2009; 54(13):4151.

    Article  CAS  PubMed  Google Scholar 

  47. Hlawitschka M, McGough RJ, Ferrara KW, Kruse DE. Fast ultrasound beam prediction for linear and regular two-dimensional arrays. IEEE Trans Ultrason, Ferroelectr, Freq Control. 2011; 58(9):2001–12.

    Article  Google Scholar 

  48. Zhu Y, Szabo TL, McGough RJ. A comparison of ultrasound image simulations with FOCUS and FieldII. In: Ultrasonics Symposium (IUS), 2012 IEEE International. Dresden: IEEE: 2012. p. 1694–1697.

    Google Scholar 

  49. McGough RJ. Rapid calculations of time-harmonic nearfield pressures produced by rectangular pistons. J Acoust Soc Am. 2004; 115:1934.

    Article  PubMed Central  PubMed  Google Scholar 

  50. Kelly JF, McGough RJ. A time-space decomposition method for calculating the nearfield pressure generated by a pulsed circular piston. IEEE Trans Ultrason, Ferroelectr, Freq Control. 2006; 53(6):1150–9.

    Article  Google Scholar 

  51. Chen D, Kelly JF, McGough RJ. A fast near-field method for calculations of time-harmonic and transient pressures produced by triangular pistons. J Acoust Soc Am. 2006; 120:2450.

    Article  PubMed Central  PubMed  Google Scholar 

  52. Chen D, McGough RJ. A 2D fast near-field method for calculating near-field pressures generated by apodized rectangular pistons. J Acoust Soc Am. 2008; 124(3):1526–37.

    Article  PubMed Central  PubMed  Google Scholar 

  53. Zeng X, McGough RJ. Optimal simulations of ultrasonic fields produced by large thermal therapy arrays using the angular spectrum approach. J Acoust Soc Am. 2009; 125(5):2967–77.

    Article  PubMed Central  PubMed  Google Scholar 

  54. Taflove A, Hagness S. Computational electrodynamics: the finite-difference time-domain method, 3rd edn. London: Artech House; 2005.

    Google Scholar 

  55. Azhari H. Basics of biomedical ultrasound for engineers. Hoboken: Wiley; 2010.

    Book  Google Scholar 

  56. Song J, Pulkkinen A, Huang Y, Hynynen K. Investigation of standing-wave formation in a human skull for a clinical prototype of a large-aperture, transcranial mr-guided focused ultrasound (MRgFUS) phased array: an experimental and simulation study. IEEE Trans Biomed Eng. 2012; 59(2):435–44.

    Article  PubMed Central  PubMed  Google Scholar 

  57. Fink M, Prada C, Wu F, Cassereau D. Self focusing with time reversal mirror in inhomogeneous media. In: Proc. IEEE Ultrason. Symp. Montreal, Que: 1989. p. 681–6.

  58. Fink M. Time reversal of ultrasonic fields. I. Basic principles. IEEE Trans Ultrason, Ferroelectr, Freq Control. 1992; 39(5):555–66.

    Article  CAS  Google Scholar 

  59. Fink M, Montaldo G, Tanter M. Time-reversal acoustics in biomedical engineering. Annu Rev Biomed Eng. 2003; 5(1):465–97.

    Article  CAS  PubMed  Google Scholar 

  60. El-Brawany M, Nassiri D, Terhaar G, Shaw A, Rivens I, Lozhken K. Measurement of thermal and ultrasonic properties of some biological tissues. J Med Eng Technol. 2009; 33(3):249–56.

    Article  CAS  PubMed  Google Scholar 

  61. NCRP. Biological effects of ultrasound: mechanisms and clinical implications. In: Report No. 074 - Biological Effects of Ultrasound: Mechanisms and Clinical Implications. Bethesda, Maryland, USA: National Council on Radiation Protection & Measurements: 1983.

  62. Bamber J, Hill C. Ultrasonic attenuation and propagation speed in mammalian tissues as a function of temperature. Ultrasound Med Biol. 1979; 5(2):149–57.

    Article  CAS  PubMed  Google Scholar 

  63. Duck FA. Physical properties of tissue: a comprehensive reference book. Waltham: Academic Press; 1990.

    Google Scholar 

  64. Goss S, Johnston R, Dunn F. Comprehensive compilation of empirical ultrasonic properties of mammalian tissues. J Acoust Soc Am. 1978; 64:423.

    Article  CAS  PubMed  Google Scholar 

  65. Goss S, Frizzell L, Dunn F. Ultrasonic absorption and attenuation in mammalian tissues. Ultrasound Med Biol. 1979; 5(2):181–6.

    Article  CAS  PubMed  Google Scholar 

  66. Hasgall PA, Di Gennaro F, Baumgartner C, Neufeld E, Gosselin MC, Payne D, Klingenböck A, Kuster N. IT’IS Database for thermal and electromagnetic parameters of biological tissues. Version 2.6, January 13th. 2015.

  67. He X, Bischof JC. Quantification of temperature and injury response in thermal therapy and cryosurgery. Crit Rev Biomed Eng. 2003; 31(5–6):355–422.

    Article  PubMed  Google Scholar 

  68. Stauffer P, Goldberg S. Introduction: thermal ablation therapy. Int J Hyperthermia. 2004; 20(7):671–7.

    Article  CAS  PubMed  Google Scholar 

  69. Emery AF, Sekins KM. The use of heat transfer principles in designing optimal diathermy and cancer treatment modalities. Int J Heat Mass Tran. 1982; 25(6):823–34.

    Article  Google Scholar 

  70. Laakso I, Hirata A. Dominant factors affecting temperature rise in simulations of human thermoregulation during RF exposure. Phys Med Biol. 2011; 56(23):7449.

    Article  PubMed  Google Scholar 

  71. Murbach M, Neufeld E, Capstick M, Kainz W, Brunner DO, Samaras T, et al.Thermal tissue damage model analyzed for different whole-body SAR and scan durations for standard MR body coils. Magn Reson Med. 2014; 71(1):421–31.

    Article  PubMed  Google Scholar 

  72. Suetens P. Fundamentals of medical imaging: Cambridge university press; 2009. ISBN: 9780521519151.

  73. Aubry JF, Tanter M, Pernot M, Thomas JL, Fink M. Experimental demonstration of noninvasive transskull adaptive focusing based on prior computed tomography scans. J Acoust Soc Am. 2003; 113:84.

    Article  CAS  PubMed  Google Scholar 

  74. Pulkkinen A, Werner B, Martin E, Hynynen K. Numerical simulations of clinical focused ultrasound functional neurosurgery. Phys Med Biol. 2014; 59(7):1679.

    Article  PubMed Central  PubMed  Google Scholar 

Download references


This work was financed by the Swiss National Center of Competence in Research (NCCR) CO-ME. The authors would also like to thank Dr. Andreas Christ, for his valuable feedback and critical assessment of this work.

Author information

Authors and Affiliations


Corresponding author

Correspondence to Adamos Kyriakou.

Additional information

Competing interests

The authors declare that they have no competing interests.

Authors’ contributions

AK was the PhD student responsible for developing the acoustic solvers, performing the computational studies, data analysis, and writing the majority of the manuscript. EN contributed several subparts of the acoustic solvers’ code, conceived the amplitude correction scheme used in the SPAC approach, developed the thermal solver used for thermal simulations, as well as the software used to segment the MR data and create the anatomical model, and contributed to the analysis and manuscript writing. BW provided clinical background, in-depth technical knowledge on FUS applicator systems, currently employed aberration-correction approaches, information on patient positioning, target selection, and helped in the interpretation of the resulting pressure and temperature distributions. GS contributed to the development of the medical image-based modeling framework and the generation of the anatomical models. NK led the design of all solver validation approaches and provided insight in the interpretation of the collected data. All authors meticulously and critically reviewed the manuscript, amended different sections pertaining to their respective fields of expertise, and approved of its final version.

Rights and permissions

This is an Open Access article distributed under the terms of the Creative Commons Attribution License (, which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly credited. The Creative Commons Public Domain Dedication waiver ( applies to the data made available in this article, unless otherwise stated.

Reprints and Permissions

About this article

Verify currency and authenticity via CrossMark

Cite this article

Kyriakou, A., Neufeld, E., Werner, B. et al. Full-wave acoustic and thermal modeling of transcranial ultrasound propagation and investigation of skull-induced aberration correction techniques: a feasibility study. J Ther Ultrasound 3, 11 (2015).

Download citation

  • Received:

  • Accepted:

  • Published:

  • DOI:


  • Acoustic modeling
  • Thermal modeling
  • Transcranial focused ultrasound
  • Focusing
  • Treatment envelope
  • Aberration correction
  • Treatment planning framework