Rapid fullwave phase aberration correction method for transcranial highintensity focused ultrasound therapies
 Scott Almquist^{1}Email authorView ORCID ID profile,
 Dennis L. Parker^{2, 3} and
 Douglas A. Christensen^{4, 5}
DOI: 10.1186/s4034901600747
© The Author(s). 2016
Received: 27 November 2015
Accepted: 13 October 2016
Published: 8 December 2016
Abstract
Background
Noninvasive highintensity focused ultrasound (HIFU) can be used to treat a variety of disorders, including those in the brain. However, the differences in acoustic properties between the skull and the surrounding soft tissue cause aberrations in the path of the ultrasonic beam, hindering or preventing treatment.
Methods
We present a method for correcting these aberrations that is fast, fullwave, and allows for corrections at multiple treatment locations. The method is simulationbased: an acoustic model is built based on highresolution CT scans, and simulations are performed using the hybrid angular spectrum (HAS) method to determine the phases needed for correction.
Results
Computation of corrections for clinically applicable resolutions can be achieved in approximately 15 min. Experimental results with a plastic model designed to mimic the aberrations caused by the skull show that the method can recover 95 % of the peak pressure obtained using hydrophonebased timereversal methods. Testing using an ex vivo human skull flap resulted in recovering up to 70 % of the peak pressure at the focus and 61 % when steering (representing, respectively, a 1.52 and 1.19fold increase in the peak pressure over the uncorrected case). Additionally, combining the phase correction method with rapid HAS simulations allows evaluation of such treatment metrics as the effect of misregistration on resulting pressure levels.
Conclusions
The method presented here is able to rapidly compute phases required to improve ultrasound focusing through the skull at multiple treatment locations. Combining phase correction with rapid simulation techniques allows for evaluation of various treatment metrics such as the effect of steering on pressure levels. Since the method computes 3D pressure patterns, it may also be suitable for predicting offfocus hot spots during treatments—a primary concern for transcranial HIFU. Additionally, the plasticskull method presented here may be a useful tool in evaluating the effectiveness of phase correction methods.
Keywords
MRgHIFU Transcranial HIFU Phase aberrations Ultrasound simulation methodsBackground
Highintensity focused ultrasound (HIFU) is becoming prominent as a treatment method for a large number of diseases. The lack of ionizing radiation and noninvasive nature of HIFU makes it well suited as a primary or supplementary treatment for uterine fibroids [1], prostate cancer [2, 3], bone metastasis [4], and other afflictions [5]. In the brain, transcranial HIFU is being investigated for neurostimulation [6], to treat movement disorders [7] and gliomas [8], and may be useful in other treatments [9].
Transcranial HIFU, however, faces challenges not present for treatments in other areas in the body. The large attenuation (10 dB on average at 0.5 MHz and 20 dB at 1.5 MHz [10]) caused by the skull requires the use of largeaperture transducers to evenly spread the power over the skull to prevent skin burns [11]. Additionally, the large difference in acoustic properties of the skull and surrounding soft tissues—notably the speed of sound, which is approximately 2740 m/s in the skull [10] and near 1500 m/s in most soft tissues [12]—introduces phase aberrations that cause the focus to be distorted, displaced, and of lower intensity. Nevertheless, aberrations can be corrected by phasing the elements of a phasedarray transducer such that constructive interference is achieved at the intended treatment location.
There exist several methods for determining the phases required to correct for aberrations [13]. The methods can be broadly classified into timereversal, energybased, and simulationbased methods. One version of the timereversal method [14] relies on an implantable hydrophone to measure the phases from each element in order to determine the inverse phase for correction. Although this offers excellent correction, the invasive nature makes it inappropriate for HIFU treatments. Energybased methods use imaging to determine the effect of the skull on the ultrasonic beam. For example, cavitation can create shock waves that can be imaged noninvasively (“ultrasonic stars” [15]), but the energy required to create cavitation and the bubbles themselves may pose a risk to the patient. Methods based on MR acoustic radiation force imaging [16] avoid cavitation but may require more time to converge to the correct phases. Simulationbased methods predict the phases through numerical simulation. Often, these simulations are based on acoustic properties derived from a CT scan [17], although using MR images of the skull obtained using UTE sequences may be possible [18]. Finitedifference timedomain (FDTD) simulations took about 2 h to compute phase corrections [19]. Faster simulations can be achieved by simplifying the acoustic models used [20], but the resulting phases will have reduced accuracy. It has been shown that fullwave simulation methods have an improved accuracy over ray tracing and can be computed within minutes [21]; however, the technique used there was limited to a single point of correction and did not include experimental results. Another work [22] has also shown fullwave methods to produce superior results compared to ray tracing at the expense of computational time (hours vs. seconds). A review of phase correction methods can be found in [13].
In this paper, we introduce a simulationbased fullwave phase aberration correction method. The method leverages the speed of the hybrid angular spectrum (HAS) method of acoustic simulation [23] in order to rapidly determine the phases required. HAS is implemented on a graphics processing unit (GPU) in order to quickly determine the acoustic pressure field in a 3D volume for each element, allowing for phase correction at multiple treatment locations with negligible increase in computational time.
We demonstrate the method’s effectiveness experimentally using two models. First, we correct the phase aberrations for a 3Dprinted plastic aberrator model that is designed to simulate the phase aberrations caused by a human skull, but for which the acoustic properties can be accurately determined. Second, phase corrections are performed using a section of an ex vivo human skull. We further expand on the capabilities of this method by combining it with rapid HAS simulations to demonstrate and evaluate the effects of transducer/model misregistration and the effects of steering on pressure levels with reasonable calculation times.
Methods
Simulation technique
The phase aberration correction method in this paper simulates acoustic pressure patterns using the HAS method [23]. HAS alternates between the space and spatialfrequency domains as it propagates a steady state wave through transverse sections of a model to simulate 3D ultrasonic beam patterns; it simulates refraction, reflection, and absorption. HAS demonstrates a computational speed increase over other simulation methods. For example, a 301 × 301 × 300voxel model with 0.15mm isotropic resolution with three tissue types computed on a 2.67GHz i7Core Windows desktop with 12 GB of RAM using MATLAB 7.10 took 46 s in HAS versus 467 min using a FDTD technique using the same input source [23], though this did not exploit parallelism for either HAS or the FDTD technique. The HAS method models radiating boundary conditions, which was also modeled in the FDTD simulation.
where α and β are the directional cosines (α = λf _{ x } and β = λf _{ y }), b _{ n } ′ is an average propagation constant across the plane, and Δz is the propagation distance from one plane to the next. The inverse Fourier transform of A _{ n } will be the pressure pattern at the face of the next plane, p _{ n }. The process can be repeated for each slice along the model.
HAS is able to account for inhomogeneities at a voxelbyvoxel level and, if needed, is able to account for reflections by calculating the impedance mismatch between voxels during the spacedomain step. Although HAS is capable of modeling scattering [24], it was not employed for these simulations. HAS does not currently model nonlinearities or shear waves.
For fast computation that would be required in a clinical setting, each element’s pressure pattern is computed in parallel on a NVIDIA Tesla GPU (NVIDIA, Santa Clara, CA) using Jacket (ArrayFire, Atlanta, GA) and MATLAB 2012 (MathWorks, Natick, MA). For small acoustic models (229 × 159 × 182 voxels with six tissue types and 256 transducer elements), pressure patterns and phase corrections can be computed in approximately 45 s. Models such as the more clinically relevant skull flap model (421 × 648 × 170 voxels with 3000 distinct types of bone and 256 elements) take approximately 15 min for pressure patterns and corrections.
Aberrator model
Acoustic properties of phase aberrator model
Speed of sound  Attenuation 

2492 m/s  4.72 dB/(cm MHz) 
Ex vivo skull model
To assess what level of corrections could be expected in an actual skull, an ex vivo section of human skull approximately 15 cm × 10 cm in size was evaluated. The skull flap was discarded, deidentified, and obtained from a deceased patient, and therefore IRB exempt. It was cleaned and frozen several months before being imaged in a clinical CT scanner at a resolution of 0.46 × 0.46 × 0.3 mm. An acoustic model of the skull was built using the CT images and a published method of mapping CT Hounsfield units (HU) to acoustic properties (speed of sound, density, and attenuation). The method, described in Pichardo et al. [27], maps a linear relationship between density and HU. The speed of sound and attenuation are determined from HU using a series of curves that were created by optimizing simulations to match experimental data. Before phase correction experiments, the skull flap was rehydrated overnight inside a degasser to remove air.
Hydrophone scans
Results
Aberrator model
Figure 5, middle and bottom rows, shows the corrected experimental and simulated pressure patterns though the aberrator model using phases found from both the hydrophone measurements and the simulationbased method.
Phase correction improvements in the aberrator model
Ratio of maximum pressure corrected to uncorrected  Distance from focus (mm)  

Uncorrected  1  1.90 
Corrected, simulationbased  1.41  0.56 
Corrected hydrophonebased  1.50  0.25 
Skull flap
After applying simulationbased phase corrections, significant improvements in the beam focus were observed. The computation time to compute the phases required for correction using simulations was approximately 15 min.
Phase correction improvements in ex vivo skull
Ratio of maximum pressure corrected to uncorrected  Distance from focus (mm)  

Uncorrected  1  1.77 
Corrected, simulationbased  1.51  0.71 
Corrected, hydrophonebased  2.17  0.25 
Further simulationbased analysis
The HAS method combined with phase correction allows for rapid evaluation of a number of different parameters. Of particular clinical interest are the patientspecific responses to treatments—including the ability to electronically steer the beam while maintaining adequate heating capability and the sensitivity of the system to errors in misregistration.
Although the results here cannot be generalized to other skulls, they provide a proofofconcept of the ability of HAS in combination with phase correction procedures to rapidly simulate possible treatment scenarios. Each graph in Figs. 8 and 9 represents 529 individual simulations and can be created in less than 4 h using a Windows computer with an i5core processor and 16 GB of RAM using MATLAB 2015a. Note that these simulations were not run in parallel—doing so would provide additional speed. Each individual simulation took approximately 20 s for the skull flap and 25 s for the aberrator model.
Discussion
The phase correction method introduced here was able to rapidly calculate the phases required to correct for aberrations. In the case where the acoustic properties were fully known (for the 3Dprinted phase aberrator model), the simulationbased correction was able to recover 95 % of the peak pressure achieved with hydrophone time reversal for targeting the geometric focus (as seen in Table 2). The relatively small difference between the simulationbased and hydrophonebased corrections for the aberrator model (95 % of the peak pressure was recovered) is most likely due primarily to registration errors. The simulations require close matching between the simulated and experimental placement of the aberrator with respect to the transducer elements. Misregistration will lead to poorer corrections. Figure 9 demonstrates this effect with translational errors in a single plane. Relatively small misregistrations can cause a significant decrease in the pressures seen. For example, simulations on the misregistration for the aberrator model show that a 1mm translation of the model results in a 1.1–7.8 % loss in the peak pressure. Misregistration causes less of an effect for the skull model, which has smoother variations over the area through which the beam travels. Figure 10 further demonstrates this fact. When there is longer correlation distance for the phase lengths, one would expect misregistration in that direction to be less disruptive, a factor for phase correction. The correlation pattern for the skull flap shows longer correlation distances along the yaxis, which agrees with the data in Fig. 9 showing more tolerance to misregistration along that axis. Furthermore, the skull flap displays broader areas of correlation than the aberrator model, again agreeing with the results in Fig. 9 that show that the skull flap is more tolerant to misregistration overall. The 50 % contour in the results for the skull flap encompasses an area of 205 mm^{2} while the area for the aberrator model covers 24 mm^{2}.
The hydrophone timereversal method is performed experimentally with the aberrator in place and is not subject to the registration requirements. However, there may be slight errors in the location of the transducer element positions or the power output of each element. We estimate a misregistration error of less than 0.25 mm in each direction for these experiments. Data from the aberrator model, where the acoustic properties are homogeneous and known, agree with this: 95 % of the pressure is recovered, which also agrees with the simulated misregistration data presented in Fig. 9.
The 3Dprinted aberrator model introduced here is a useful tool for assessing the ability of phase correction methods in the absence of acoustic parameter uncertainties. In the skull flap, where acoustic parameter uncertainties are present, simulationbased corrections were only able to recover 70 % of the pressure achieved using hydrophone timereversal when targeting to the geometric focus (and 61 % when steering) compared to 95 % in the aberrator model. The aberrator model could be used to evaluate other possible methods of phase correction (ray tracing or FDTD timereversal, for example) and to compare using different transducers or setups without needing to account for possible errors in acoustic modeling.
The difference between the percentage of pressure recovered using the simulationbased and hydrophone timereversal methods of phase correction in the aberrator model (95 %) and the human skull flap (70 %) demonstrates the necessity for accurate acoustic parameter estimation of the skull. In transcranial HIFU, a large barrier to treatments is the necessary power required to achieve sufficient heating at the focal location without causing damage to unintended locations. Although the acoustic model of the skull flap used in this paper showed a 1.52fold increase in peak power over uncorrected results, the hydrophone timereversal phase correction demonstrated a 2.17fold improvement. This suggests that there may be benefit in further research into acoustic parameter estimation of the skull. The aberrator model parameters were determined via a throughtransmission test on a rectangular block of the same material used to 3D print the model, while the bone was modeled using the method presented in [27], which assigned values based on CT images. However, there is some evidence that the method presented in [27] may not be accurate in all cases as simulations using the method deviated in temperature by 24 % compared to the observed data, and the peak focal point distance between the simulations and experimental measurements was off by an average of 1.6 mm [28].
Other noninvasive methods have shown recovery of 14–58 % of peak intensity of hydrophonebased corrections using human skull flaps [29] compared to this method’s recovery of 49 % of the peak intensity (70 % of the peak pressure) at the geometric focus. However, caution must be used in directly comparing results for different transducer setups and skulls. Larger aperture size will usually cause the ultrasound to travel through a larger, more inhomogeneous region, yielding a greater need for phase correction. Additionally, a change in frequency of the transducer will change the precision required for phase correction. Unquantified variations in acoustic properties between human skulls also make direct comparisons difficult. Ideally, phase correction methods should be evaluated on the same system using an aberrator model similar to the one used here to evaluate effectiveness absent these disparities.
The method introduced here is more computationally efficient when compared to other methods presented in the literature. Simulations on a GPU for the skull flap modeled with clinically relevant resolution were performed in approximately 15 min. This computational time is for the full 3D volume for each individual element, allowing for corrections at multiple treatment locations. This computational time represents an eightfold improvement over a similar FDTDbased timereversal simulation that was performed for only a single treatment location [19].
The method described has the benefit of simulating a full 3D volume. While we have used the phase information to demonstrate the possibility for multiple treatment locations, it would also be possible to use the amplitude information for predicting offfocus hot spots. For transcranial treatments, this could allow risk assessment of skin burns or heating in important brain structures. Combining the phase correction method with rapid simulations could allow for many clinically valuable insights, including determining which patients may not be suitable for transcranial HIFU treatments, treatment envelope evaluation, establishing limits for the maximum allowable power, or treatment planning.
Conclusions
We have presented a simulationbased method for phase correction that allows for improved focusing in transcranial HIFU applications. The method is rapid; phases can be computed in approximately 15 min for the resolution required for transcranial applications, and it creates corrections for a 3D volume allowing for treatment planning. The method was tested on both a 3Dprinted acoustic aberrator and an ex vivo human skull flap where it resulted in increased pressure, decreased beam width, and reduction in side lobes. Corrections were achieved for both targeting to the geometric focus and steering to a target 5 mm away from the focus. Additionally, the method of corrections was combined with fast ultrasound simulations to evaluate the effects of steering and transducer/model misregistration on pressure levels. Since simulations are performed for a 3D volume, future work will include predicting maximum safe power usage and offfocus hot spots produced during treatments.
Abbreviations
 FDTD:

Finitedifference timedomain
 HAS:

Hybrid angular spectrum
 HIFU:

Highintensity focused ultrasound
 HU:

Hounsfield units
Declarations
Acknowledgements
The authors wish to thank Henrik Odéen and Allison Payne for their insight and comments relating to this work as well as Robb Merrill for his help in designing and manufacturing the holder used for registration.
Funding
Funding for this study was provided by the Focused Ultrasound Foundation and NIH R01 grants EB013433 and CA172787.
Availability of data and materials
Please contact author for data requests.
Authors’ contributions
SA performed the experiments and drafted the manuscript. SA and DAC developed and implemented the phase correction method. SA, DAC, and DLP conceived of and designed the study. All authors read and approved the manuscript.
Competing interests
The authors declare that they have no competing interests.
Consent for publication
Not applicable
Ethics approval and consent to participate
Tissues used in this study were deidentified and obtained from a deceased patient, and therefore IRB exempt.
Open AccessThis article is distributed under the terms of the Creative Commons Attribution 4.0 International License (http://creativecommons.org/licenses/by/4.0/), which permits unrestricted use, distribution, and reproduction in any medium, provided you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons license, and indicate if changes were made. The Creative Commons Public Domain Dedication waiver (http://creativecommons.org/publicdomain/zero/1.0/) applies to the data made available in this article, unless otherwise stated.
Authors’ Affiliations
References
 Okada A, Morita Y, Fukunishi H, Takeichi K, Murakami T. Noninvasive magnetic resonanceguided focused ultrasound treatment of uterine fibroids in a large Japanese population: impact of the learning curve on patient outcome. Ultrasound Obstet Gynecol. 2009;34:579–83.View ArticlePubMedGoogle Scholar
 Murat FJ, Poissonnier L, Rabilloud M, Belot A, Bouvier R, Rouviere O, Chapelon JY, Gelet A. Midterm results demonstrate salvage highintensity focused ultrasound (HIFU) as an effective and acceptably morbid salvage treatment option for locally radiorecurrent prostate cancer. Eur Urol. 2009;55:640–9.View ArticlePubMedGoogle Scholar
 Poissonnier L, Chapelon JY, Rouvière O, Curiel L, Bouvier R, Martin X, Dubernard JM, Gelet A. Control of prostate cancer by transrectal HIFU in 227 patients. Eur Urol. 2007;51:381–7.View ArticlePubMedGoogle Scholar
 Li C, Zhang W, Fan W, Huang J, Zhang F, Wu P. Noninvasive treatment of malignant bone tumors using highintensity focused ultrasound. Cancer. 2010;116:3934–42.View ArticlePubMedGoogle Scholar
 Jolesz FA, McDannold N. Current status and future potential of MRIguided focused ultrasound surgery. J Magn Reson Imaging. 2008;27:391–9.View ArticlePubMedGoogle Scholar
 King RL, Brown JR, Newsome WT, Pauly KB. Effective parameters for ultrasoundinduced in vivo neurostimulation. Ultrasound Med Biol. 2013;39:312–31.View ArticlePubMedGoogle Scholar
 Elias WJ, Huss D, Voss T, Loomba J, Khaled M, Zadicario E, Frysinger RC, Sperling SA, Wylie S, Monteith SJ, Druzgal J, Shah BB, Harrison M, Wintermark M. A pilot study of focused ultrasound thalamotomy for essential tremor. N Engl J Med. 2013;369:640–8.View ArticlePubMedGoogle Scholar
 McDannold N, Clement G, Black P, Jolesz F, Hynynen K. Transcranial MRIguided focused ultrasound surgery of brain tumors: initial findings in three patients. Neurosurgery. 2010;66:323–32.View ArticlePubMedPubMed CentralGoogle Scholar
 Monteith S, Sheehan J, Medel R, Wintermark M, Eames M, Snell J, Kassell NF, Elias WJ. Potential intracranial applications of magnetic resonance–guided focused ultrasound surgery: a review. J Neurosurg. 2013;118:215–21.View ArticlePubMedGoogle Scholar
 Fry FJ, Barger JE. Acoustical properties of the human skull. J Acoust Soc Am. 1978;63:1576–90.View ArticlePubMedGoogle Scholar
 Clement GT, White J, Hynynen K. Investigation of a largearea phased array for focused ultrasound surgery through the skull. Phys Med Biol. 2000;45:1071.View ArticlePubMedGoogle Scholar
 Fairhead AC. ICRU report 61: providing reference data for tissue properties. J Acoust Soc Am. 1999;105:1324.View ArticleGoogle Scholar
 Kyriakou A, Neufeld E, Werner B, Paulides MM, Szekely G, Kuster N. A review of numerical and experimental compensation techniques for skullinduced phase aberrations in transcranial focused ultrasound. Int J Hyperthermia. 2014;30:36–46.View ArticlePubMedGoogle Scholar
 Fink M, Cassereau D, Derode A, Prada C, Roux P, Tanter M, Thomas JL, Wu F. Timereversed acoustics. Rep Prog Phys. 2000;63:1933.View ArticleGoogle Scholar
 Pernot M, Montaldo G, Tanter M, Fink M. “Ultrasonic stars” for timereversal focusing using induced cavitation bubbles. Appl Phys Lett. 2006;88:34102134102–3.View ArticleGoogle Scholar
 Vyas U, Kaye E, Pauly KB. Transcranial phase aberration correction using beam simulations and MRARFI. Med Phys. 2014;41:32901.View ArticleGoogle Scholar
 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–93.View ArticlePubMedGoogle Scholar
 Miller GW, Eames M, Snell J, Aubry JF. Ultrashort echotime MRI versus CT for skull aberration correction in MRguided transcranial focused ultrasound: In vitro comparison on human calvaria. Med Phys. 2015;42:2223–33.View ArticlePubMedGoogle Scholar
 Marquet F, Pernot M, Aubry JF, Montaldo G, Marsac L, Tanter M, Fink M. Noninvasive transcranial ultrasound therapy based on a 3D CT scan: protocol validation and in vitro results. Phys Med Biol. 2009;54:2597.View ArticlePubMedGoogle Scholar
 Clement GT, Hynynen K. Correlation of ultrasound phase with physical skull properties. Ultrasound Med Biol. 2002;28:617–24.View ArticlePubMedGoogle Scholar
 Kyriakou A, Neufeld E, Werner B, Székely G, Kuster N. Fullwave acoustic and thermal modeling of transcranial ultrasound propagation and investigation of skullinduced aberration correction techniques: a feasibility study. J Ther Ultrasound. 2015;3:11.View ArticlePubMedPubMed CentralGoogle Scholar
 Jones RM, Hynynen K. Comparison of analytical and numerical approaches for CTbased aberration correction in transcranial passive acoustic imaging. Phys Med Biol. 2016;61:23.View ArticlePubMedGoogle Scholar
 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:1093–100.View ArticlePubMedGoogle Scholar
 Christensen D, Almquist S. Incorporating tissue absorption and scattering in rapid ultrasound beam modeling. Proc SPIE. 2013;8584:85840X–1–7.View ArticleGoogle Scholar
 Gélat P, ter Haar G, Saffari N. A comparison of methods for focusing the field of a HIFU array transducer through human ribs. Phys Med Biol. 2014;59:3139–71.View ArticlePubMedGoogle Scholar
 Kak AC, Dines KA. Signal processing of broadband pulsed ultrasound: measurement of attenuation of soft biological tissues. IEEE Trans Biomed Eng. 1978;BME25:321–44.View ArticleGoogle Scholar
 Pichardo S, Sin VW, Hynynen K. Multifrequency characterization of the speed of sound and attenuation coefficient for longitudinal transmission of freshly excised human skulls. Phys Med Biol. 2011;56:219.View ArticlePubMedGoogle Scholar
 Pulkkinen A, Werner B, Martin E, Hynynen K. Numerical simulations of clinical focused ultrasound functional neurosurgery. Phys Med Biol. 2014;59:1679.View ArticlePubMedPubMed CentralGoogle Scholar
 Clement GT, Hynynen K. A noninvasive method for focusing ultrasound through the human skull. Phys Med Biol. 2002;47:1219–36.View ArticlePubMedGoogle Scholar