- Open Access
MR thermometry for focused ultrasound monitoring utilizing model predictive filtering and ultrasound beam modeling
Journal of Therapeutic Ultrasoundvolume 4, Article number: 23 (2016)
A major challenge in using magnetic resonance temperature imaging (MRTI) to monitor focused ultrasound (FUS) applications is achieving high spatio-temporal resolution over a large field of view (FOV). This is important to accurately monitor all ultrasound (US) power depositions. Magnetic resonance (MR) subsampling in conjunction with thermal model-based reconstruction of the MRTI utilizing Pennes bioheat transfer equation (PBTE) is one promising approach. The thermal properties used in the thermal model are often estimated from a pre-treatment, low-power sonication.
In this proof-of-concept study we investigate the use of US simulations computed using the hybrid angular spectrum (HAS) method to estimate the US power deposition density Q, thereby avoiding the pre-treatment sonication and any potential tissue damage. MRTI reconstructions are performed using a thermal model-based reconstruction method called model predictive filtering (MPF). Experiments are performed in a homogeneous gelatin phantom and in a gelatin phantom with embedded plastic skull. MPF reconstructions are compared to separate sonications imaged with fully sampled data over a smaller FOV. Temperature root-mean-square errors (RMSE) and focal spot positions and shapes are evaluated.
HAS simulations accurately predict the location of the focal spot (to within 1 mm) in both phantoms. Accurate temperature maps (RMSE below 1 °C), where the location of the focal spot agrees well with fully sampled “truth” (to within 1 mm), are also achieved in both phantoms.
HAS simulations can be used to accurately predict the focal spot location in homogeneous media and when focusing through an aberrating plastic skull. The HAS simulated power deposition (Q) patterns can be used in the MPF thermal model-based reconstruction to obtain accurate temperature maps with high spatio-temporal resolution over large FOVs.
Focused ultrasound (FUS) treatments are currently being performed under both ultrasound (US) and magnetic resonance imaging (MRI) guidance. One of the main advantages of MRI guidance is the ability to accurately measure and monitor the temperature evolution in real-time using MR temperature imaging (MRTI) [1, 2]. MRTI enables evaluation of the treatment target site in terms of measuring temperature rise and thermal dose, as well as safety monitoring of the surrounding healthy tissues in the US near and far field. The current gold standard in MRTI is the proton-resonance frequency shift (PRFS) method, which relies on the temperature-dependent shielding of the hydrogen nucleus from the surrounding electron cloud. The shielding effect increases with temperature as hydrogen bonds between water molecules stretch, bend, and break, resulting in a lower resonance frequency with increased temperature [3–6]. The main advantages of the PRFS method are its relatively high sensitivity, that it can be monitored with readily available MRI pulse sequences, and that it is to a large degree tissue type independent for most aqueous-based tissues . One drawback of the PRFS method is that it does not work in lipid-based tissues since these lack the necessary hydrogen bondings.
One of the main challenges with MRTI is to achieve accurate temperature imaging with sufficient temporal and spatial resolution  over a large enough field of view (FOV) to monitor all surrounding and intervening tissues for unwanted heating. This is especially challenging in FUS applications since the US energy can be focused deep into the body. Attempts to increase the acquisition speed and/or volume coverage for MRTI have included using fast pulse sequences such as echo planar imaging (EPI), where multiple k-space lines are acquired after each radiofrequency (RF) excitation. Both single- and multi-slice 2D as well as fully 3D versions, with both single-shot and segmented phase encodings, have been investigated [9–15]. Other fast pulse sequences that have been investigated include balanced steady-state free precession (bSSFP)-type sequences which have a relatively short repetition time (TR) [16, 17], as well as utilizing signal-to-noise ratio (SNR)-efficient non-Cartesian approaches [18–21]. Speedups have also been achieved by combining subsampled k-space data acquisition with subsampled data reconstruction methods such as parallel imaging [22–26], compressed sensing [27–29], and temporally constrained reconstruction approaches [30, 31]. Filtering approaches, such as Kalman filters, and fitting and modeling methods, such as thermal modeling utilizing the Pennes bioheat transfer equation (PBTE), have been used to reconstruct temperature images based upon limited k-space data samples or as a mean of noise reduction in fully sampled data [19, 32–35].
Model predictive filtering (MPF) is a previously published thermal model-based method for reconstructing subsampled MRTI data for FUS applications. MPF uses thermal modeling through the PBTE to forward predict temperature change. The temperature change is converted to a change in the MR phase image and then used to fill in missing data in k-space from the subsampled acquisition. MPF can achieve high spatio-temporal resolution over large FOVs, but the thermal modeling requires estimates of both tissue thermal and acoustic properties. These properties, such as thermal conductivity, k (in W/m/°C), and deposited US power density, Q (in W/m3), can be estimated from a pre-treatment low-temperature rise sonication [36, 37]. While these test sonications currently serve purposes to localize the focal spot and ensure good US coupling, they also result in additional, potentially harmful tissue heating [38–40] and may also prolong the total treatment time.
In this study, we use US modeling with the previously published hybrid angular spectrum (HAS) method [41–46] to obtain the US power deposition (Q) required for MPF [32, 47, 48]. This avoids the potentially harmful pre-treatment heating and can result in shorter total treatment time since the treatment can be planned and Q can be calculated beforehand. In this proof-of-concept study, experiments are performed in a homogeneous gelatin phantom and in a gelatin phantom with an embedded plastic skull. In the skull phantom, HAS-based US simulations are also used for phase aberration correction using a previously described method .
MR imaging and image reconstruction
All MR imaging was performed on a 3T MRI scanner (TIM Trio, Siemens Medical Solutions, Erlangen, Germany) with imaging parameters summarized in Table 1. A 3D gradient recalled echo (GRE) pulse sequence with segmented EPI readout was used for all MRTI. k-space was subsampled with a Cartesian variable-density equally spaced sampling scheme , where the sampling in the kz-slice encode direction was varied to achieve denser sampling over the k-space center while maintaining a constant echo-spacing down the echo train. The mask used to segment out the skull for the phase aberration correction (see below) was based on a higher resolution standard 3D GRE scan (i.e., no EPI read out). All MR data were zero-filled interpolated to 0.5-mm isotropic voxel spacing to minimize partial volume effects .
Two in-house-built RF receive-only single-channel loop coils (schematically shown in green in Fig. 1) with different diameters (12 cm for the homogeneous phantom and 18 cm for the skull phantom, respectively) were used for signal detection.
The MRTI data was reconstructed with both the MPF method and using a previously published compressed sensing-like temporally constrained reconstruction (TCR) algorithm for comparison [15, 30, 52]. Full descriptions of the MPF and TCR algorithms are given elsewhere [30, 32], and only a brief description as background is given here. The MPF algorithm utilizes the PBTE  given by the following:
where ρ = tissue density (kg/m3); C t , C b = specific heat of tissue, blood (J/kg/°C); T, T blood = tissue, arterial blood temperature (°C); k = thermal conductivity (W/m/°C); W = Pennes perfusion parameter (kg/m3/s); Q = US power deposition density (W/m3). This equation is used to forward predict the temperature change from one dynamic MR time frame, (n), to the next, (n + 1), in 3D. The 3D temperature change map is converted to a corresponding phase change map using the PRFS equation (see Eq. 2) and is then combined with the magnitude image from the previous dynamic time frame (n) to create a complex image. The complex image is Fourier transformed into k-space and used to fill in the parts of k-space that were not sampled during time frame (n + 1). This updated k-space is finally inverse Fourier transformed back to image space, where updated temperature maps can be calculated using the PRFS equation.
The TCR algorithm iteratively adjusts predicted images to minimize a cost function consisting of a fidelity term and a constraint term. The fidelity term ensures that the solution is faithful to the acquired subsampled k-space data, and the constraint term applies an appropriate temporal smoothness constraint, in this work utilizing the L1-norm.
All MR temperature maps were calculated with the PRFS method, where the phase difference between two images, Δφ, is scaled to yield the temperature difference, ΔT, according to:
where γ is the gyromagnetic ratio (Hz/T), α is a thermal coefficient (here assumed to be −0.01 ppm/°C), B 0 is the main magnetic field strength (T), and TE is the echo time (s). Since all FUS sonications were relatively short (<30 s, see below), no field drift correction was applied. All image reconstruction and post-processing was performed using MATLAB (R2015a, The MathWorks Inc., Natick, MA).
The subsampled MRTI data for both phantom experiments were reconstructed in three ways: (1) with MPF using Q from HAS US simulations (see below) and previously published phantom properties (ρ, C t , k) ; (2) with MPF using Q and k determined from fitting a pre-treatment, low-power sonication (i.e., low-temperature rise) using a previously published analytical solution of the heating [36, 37, 55] together with the previously published phantom properties (ρ, C t ) ; and (3) using the TCR algorithm. The three methods are referred to as MPF HAS, MPF LP analytical fit, and TCR, respectively, in the remainder of this paper. The three reconstruction methods for the subsampled data were compared to an estimate of the “true” temperature rise obtained using identical US sonications monitored with fully sampled k-space data over a smaller FOV (Table 1). The short-duration, low-temperature rise test sonications used in MPF LP analytical fit to analytically obtain Q and k were imaged with the same MR parameters and fully sampled data as the “truth” sonications (Table 1).
Both phantoms were fabricated from 250-bloom gelatin  [speed of sound = 1553 ± 10 m/s (mean ± standard deviation), US attenuation = 0.06 ± 0.01 Np/(cm*MHz), density ρ = 1057 ± 44 kg/m3, specific heat capacity C t = 3635 ± 88 J/(kg*°C), thermal conductivity k = 0.55 ± 0.01 W/(m*°C)] using previously published methods . The mean values listed were used for the US simulations (see below). The plastic skull used in this study was commercially purchased (model A20, 3B Scientific, Tucker, GA, USA) and composed of homogeneous PVC plastic with varying thickness. Tabular acoustic properties for the PVC material were used in the HAS simulations [speed of sound = 2376 m/s, US attenuation = 1.50 Np/(cm*MHz), density ρ = 1200 kg/m3 [56, 57]].
FUS sonications were performed using an MR-compatible phased-array ultrasound transducer (256 elements, 1-MHz frequency, 13-cm radius of curvature, focal spot size full width at half maximum 2 × 2 × 8 mm, Imasonic, Voray-sur-l’Ognon, France), accompanying hardware and software for mechanical positioning and electronic beam steering (Image Guided Therapy, Pessac, France), and in-house developed hardware and software (code written in MATLAB) for applying the phase aberration correction. The transducer was coupled to the phantoms with a bath of de-ionized and de-gassed water as indicated in Fig. 1.
In the homogeneous phantom, three continuous-wave sonications at 40 W for 28.7 s were performed for both the subsampled and the fully sampled truth cases, for a total of six sonications. For experimental Q and k identification using the LP analytical fit method, a low-temperature rise continuous-wave sonication at 14 W for 28.7 s was performed. Similar experiments were performed in the skull phantom, at 125 W for 23.9 s (three continuous-wave sonications each for fully sampled truth and subsampled k-space, for a total of six sonications) and at 100 W for 14.3 s (LP analytical fit).
US simulation and phase aberration correction
The US simulations were performed with the HAS method [41–43]. HAS is a full-wave simulation method that takes into account US refraction, absorption, and first-order reflection. Inhomogeneous tissues are modeled by segmenting the anatomy into a 3D Cartesian grid and assigning unique speed-of-sound, absorption coefficient, and density values for the tissue type within each voxel. The hybrid angular spectrum approach leap-frogs between the space and spatial-frequency domains, utilizing the space domain to account for individual voxel attenuation and speed-of-sound properties and the spatial-frequency domain to linearly propagate the US wave from plane-to-plane throughout the model. One of the main advantages of HAS is the short computation time, which can be up to two orders of magnitude shorter than with, e.g., finite-difference time-domain (FDTD) simulations .
Due to differences in the speed of sound in the skull and gelatin, the US waves from the different elements in the phased-array transducer will arrive at the US focus out of phase with each other. To correct these phase aberrations, HAS calculates the individual pressure pattern of each of the 256 elements in the phased-array transducer at the intended focal spot location, with an initial assumption of all elements having zero phase and the same amplitude . Maximum constructive interference at the focal spot is then achieved by compensating for the phase offsets of each element. The computation for each of the 256 transducer elements was performed in parallel on a GPU (Nvidia Tesla, Nvidia, Santa Clara, CA) for maximum computation speed, and the full phase correction simulation took approximately 5 min. Once the correct phases for all 256 elements have been found, a single HAS simulation is done to find the Q pattern, and this takes approximately 10 s (this can be compared to the time it takes to perform the analytical fit for MPF LP analytical fit, which is on the order of minutes ). Since HAS is a full-wave simulation method, US amplitude and phases are calculated in a full 3D volume, so treatment planning and phase aberration correction for multiple and arbitrary positions in the full volume can be achieved with minimal added computation time.
For the homogeneous phantom in this study, a two-compartment segmentation consisting of water and phantom material (gelatin) was performed based on 3D segmented EPI magnitude images. Based on this segmentation and the distance from the transducer to the bottom of the phantom, measured in the same images, HAS was applied to simulate the 3D Q-pattern inside the phantom. For the skull phantom study, the segmentation was based on a high-resolution 3D GRE scan (i.e., no EPI readout) and a three-compartment segmentation consisting of water, phantom material (gelatin), and skull was performed.
Evaluation of accuracy and precision
The root-mean-square error (RMSE) was used to compare the accuracy of the three different reconstruction methods. The mean of the three fully sampled datasets was used as the reference truth for comparison to all other reconstruction methods (the three truths were acquired interleaved with the three subsampled datasets). RMSEs were calculated for the hottest voxel (HV) in each sonication as well as for a global error around the focal spot, defined as all voxels having a temperature rise greater than 15 % of the hottest truth voxel, and for a more local error around the focal spot, defined to be all voxels having a temperature rise greater than 85 % of the hottest truth voxel. Percentage temperature rise relative to truth was chosen since the temperature rise in the homogeneous phantom was almost twice as high as in the skull phantom, and it was deemed that this would be a more comparable measure between the two phantoms than the RMSE for all voxels above an absolute temperature limit. Since the temporal resolution of the subsampled scans was twice as high as truth (2.4 vs. 4.8 s, see Table 1), and the temperature measured with MRTI should be assigned to the time when the center of k-space is sampled (which for segmented EPI pulse sequences occurs half way through the acquisition), linear interpolation between two consecutive subsampled scans was used to estimate the temperature at the time corresponding to when the fully sampled acquisition sampled the center of k-space. The linearly interpolated values were used when calculating the RMSE. The mean and standard deviation (SD) of the RMSE for the three repeated sonications were calculated for both phantoms.
For all sonications, the focal spot center and the full width at half maximum (FWHM) of the temperature profiles were computed and the MPF and TCR values were compared to truth. All FWHM values were found by linear interpolation within the zero-filled interpolated voxels.
Q-patterns as obtained from the HAS simulations and the analytical fits to the LP pre-treatment sonications were also calculated and compared for the two phantoms. Maximum values and FWHM were compared between the two methods (i.e., HAS and LP analytical fit), and positions were compared to truth heatings.
The mean and SD of the temperature rise of the hottest voxel versus time for the two phantoms can be seen in Fig. 2a, b, respectively. To be able to calculate the RMSE, the subsampled temperatures were linearly interpolated to agree in time with the fully sampled truth as described above, and Fig. 2a, b plots these, lower temporal resolution, interpolated values. Figure 2c, d shows errors for the two phantoms, as compared to the fully sampled truth. The TCR reconstruction can be seen to perform slightly better than the two MPF reconstructions, and this is also seen in Fig. 3, which shows the mean and SD of the RMSE for this hottest voxel for the three reconstruction methods. Figure 3 further shows the RMSE for all voxels with temperature rises greater than 15 and 85 % of the hottest truth voxels in the two phantoms. In general, TCR can be seen to be more accurate than the two MPF methods. Comparing the two MPF methods, MPF HAS is generally more accurate in the homogeneous phantom whereas MPF LP analytical fit is more accurate in the skull phantom. All RMSE for all reconstruction methods were below 1.1 °C. Figure 4 shows the spatial error over the focal spot in one of the three repeated sonications for the three reconstruction methods, at the time of maximum temperature rise.
Figure 5 shows the FWHM of the temperature profiles in the three orthogonal encoding directions for truth and the three reconstruction methods. For the homogeneous phantom, the largest in-plane (i.e., in the FE-PE-plane) difference between truth and any of the reconstruction methods is below 0.3 mm (largest for MPF LP analytical fit in FE) and the largest through-plane (i.e., in the SE-direction) difference is 1.0 mm (MPF HAS). It can be noted that all these are within the acquired voxel size (i.e., before zero filling, Table 1). For the skull phantom, the widths in all three orthogonal directions also agree well with truth, and the largest differences are 0.2 mm in-plane (MPF HAS in FE) and 0.3 mm through-plane (MPF HAS). It can further be seen that the focal spot is wider in the skull phantom than in the homogeneous phantom, which is to be expected due to phase aberration when focusing through the skull.
The distance between the hottest truth voxel and the focal spot center for the MPF and TCR temperature reconstructions is shown in Fig. 6a. Figure 6b shows the distance between the hottest truth temperature voxel and the position of the hottest Q-voxel for HAS prediction and LP analytical fit. For most cases, the focal spot positions coincide, and in the cases when the distance is non-zero, it is still within the acquired voxel dimensions.
Figure 7 shows three orthogonal views of the spatial distribution of Q for the two methods in determining Q (i.e., through HAS simulations and through LP analytical fit) for the two phantoms, and corresponding line plots through the maximum Q values. The maximum Q values and the FWHM are listed in Table 2. As can be expected, the Q in the skull phantom is considerably lower (approximately an order of magnitude) than in the homogeneous phantom since the US is transmitted through the skull, which has higher US absorption than the gelatin. The Q prediction from the LP analytical fit is also higher in both phantoms, 25 % higher in the homogeneous phantom and 38 % higher in the skull phantom. In the homogeneous phantom, the widths of the two Q-patterns agree to within the measured accuracy, but HAS simulated Q is approximately 1 mm longer in the through-plane direction. In the skull phantom, the differences are greater with the HAS-simulated Q being approximately 1 mm narrower in-plane and 1.2 mm shorter through-plane.
For all MPF reconstructions using HAS-simulated Q, a tabular value of k = 5.49 × 10−4 W/(m*°C) was used for thermal conductivity , whereas the analytical fit of the low temperature rise heatings for the homogeneous and skull phantom resulted in thermal conductivity values of k = 6.50 × 10−4 W/(m*°C) and k = 7.06 × 10−4 W/(m*°C), respectively.
The effect of the phase aberration correction can be seen in Fig. 8, where a comparison of temperature maps from identical sonications monitored with fully sampled acquisitions are shown. The FWHM of the temperature distribution in the non-corrected data was 4.8, 4.2, and 14.0 mm (in FE, PE, and SE directions), which can be compared to the narrower FWHM for the aberration corrected truth seen in Fig. 5 (4.2, 3.9, and 13.8 mm). The maximum temperature rise also increased by 8 %, from 11.0 to 11.9 °C, when the aberration correction was applied.
In this proof-of-concept study, we have investigated the use of US simulations to estimate the power deposition density Q for use in the thermal modeling of the MPF method. The temperature maps reconstructed with the HAS-simulated Q were, in general, as accurate as those reconstructed with MPF using the empirically determined Q from the LP analytical fit, and all temperature RMSE were below 1 °C. Both the position and the shape of the focal spots from the HAS-based temperature maps agreed well with fully sampled truth, and were all within the size of one acquisition voxel (before zero-filling). These results demonstrate the feasibility of using US simulations in conjunction with model-based reconstruction to achieve high spatio-temporal resolution, large-FOV MRTI.
Figures 2, 3, and 4 show that MPF HAS performs as well as or better than MPF LP analytical fit for the case of the homogeneous phantom. This mainly seems to be due to a smaller overestimation of the temperatures (Fig. 2c), as Figs. 5 and 6a show very similar widths and positions of the focal spots compared to truth. There is also a relatively small difference in the widths of the predicted Q-patterns (Table 2), and their locations coincide with within the zero-filled resolution of 0.5 mm. There is, however, a relatively large difference of 25 % in the maximum Q value (Table 2), but the higher Q estimated for the LP analytical fit solution is counteracted by an 18 % higher estimated k than the tabular value used in MPF HAS (6.50 compared to 5.49 W/(m*°C)). Comparing the offsets in temperature and in position of Q (Fig. 6a compared to Fig. 6b), it can further be observed that the temperatures more accurately predict the focal spot location since they are based on a combination of the model (which depends on Q) and the acquired MR data.
For the more challenging case of the skull phantom, MPF HAS slightly underestimates the temperatures during the heating whereas MPF LP analytical fit overestimates the temperatures (Fig. 2d), and in general, MPF LP analytical fit is more accurate (Fig. 3). This seems to be mainly due to a narrower focal spot in all three directions for MPF HAS compared to MPF LP analytical fit and truth (Fig. 5) and not due to a mismatch in the actual position of the focal spot, where all three reconstruction methods are within 0.5 mm (Fig. 6a). The underestimation in the width of the MPF HAS focal spot can be attributed to the narrower Q from the HAS simulations: approximately 1 mm in all three directions as compared to the LP analytical fit solution (Table 2). Just as for the homogeneous phantom, MPF LP analytical fit estimates higher Q than HAS (38 %) and the LP analytical fit solution correspondingly also results in a compensating higher k, in this case, approximately 29 % higher.
The current study intended to compare MPF HAS to MPF LP analytical fit and since the analytical fit method assumes a Gaussian-shaped focal spot, the comparison and error calculations were focused on the focal spot (i.e., not including off-focal regions). Although no truth for the US near- and far-field (e.g. along the “brain” surface in the skull phantom) was available in this study, and hence, RMSEs could not be calculated, it can be noted that outside of the focal spot region where Q was set to zero (and hence, the thermal modeling contribution will be zero), the MPF method will reconstruct temperature maps using a “sliding window” or “view sharing” of the subsampled k-space. This ensures that no off-focal heating will go undetected, and as long as any off-focal heating is relatively slow (i.e., low dT/dt), it can be accurately monitored with the sliding window reconstruction. For the sliding window reconstruction to work properly, the perfusion and thermal conductivity (i.e., W and k) are also set to zero for the off-focal regions, so that no heat dissipation is modeled. It can further be noted that HAS predicts power deposition for off-focal areas which could be used for MPF reconstructions of, e.g., the US near-field. The accuracy and precision of this will be investigated in future studies.
The accuracy of the position of the Q-patterns determined with HAS makes it a promising approach for pre-treatment localization of the focal spot. In many current applications, such as transcranial treatments, initial low-power sonications monitored with MRTI are performed to evaluate the exact position of the focal spot [59, 60]. With accurate US simulations, this could possibly be performed without the potentially tissue-damaging pre-treatment sonications. The achieved accuracy of approximately 0.5–1.0 mm in Fig. 6 is on the order of the current MRTI voxel size and can be deemed adequate. Alternative ways to localize the focal spot are through MR acoustic radiation force imaging (AFRI) or possibly using so called MR tracker coils. In MR-ARFI, low-energy US sonications are used to encode tissue displacement into the MR phase image, which can then be used to localize the focal spot [61–64]. Just as for MRTI, the sonication needed for ARFI can potentially be tissue damaging [65, 66]. For situations where the ultrasound transducer can be moved relative to the patient, MR tracker coils on the transducer can be used to triangulate the intended focal spot position without the need for an actual sonication [67, 68]. However, unlike the alternative methods (i.e., US simulation, MRTI, and MR-ARFI), tracker coils will not work when focusing through an aberrating medium such as the skull bone. US simulation-based focal spot localization hence has the potential to both avoid the pre-treatment sonication and work in, e.g., transcranial applications.
The fact that the magnitude of the Q-patterns can differ between 25 and 38 % for the HAS and the LP analytical fit method, as seen in Fig. 7, and that accurate temperature maps can still be reconstructed highlights the robustness of the MPF method. As was previously shown, the accuracy depends on both the amount of k-space subsampling (R) and the rate of temperature increase (dT/dt) . A relatively high subsampling factor of R = 7 was used in this study (the previous study showed that subsampling factors up to R = 12 produced accurate temperature measurements ex vivo), and a lower R could be used to achieve more accurate temperature measurements.
As expected, the plastic skull attenuated the US considerably, and both the HAS-simulated and analytically derived Q were at least an order of magnitude lower in the skull phantom than in the homogeneous phantom. A more than three times increase in US power (125 versus 40 W used in the skull and homogeneous phantoms, respectively) resulted in only about half the temperature rise (max ΔT of 11.9 versus 22.0 °C in the skull and homogeneous US phantoms, respectively). This difference of approximately 6× (3× the power resulting in 0.5× the temperature rise) agrees well with previous hydrophone studies showing an approximate 85 % drop in US intensity when focusing through the skull, compared to focusing in water .
The fact that HAS simulations did not include US scattering may be part of the reason for the narrower simulated Q compared to the Q from the LP analytical fit solution. Other factors affecting the larger discrepancies seen in the skull phantom can include uncertainties in the tabular acoustic values used for the skull, where, for example, an underestimation of the speed of sound would result in an underestimated wavelength and hence beam width. The accuracy of the skull segmentation, which for in vivo applications routinely is performed based on high-resolution computer tomography (CT) images with sub-millimeter resolution , can also be assumed to play a role.
This study was performed in relatively simple phantoms, and further in vivo validation studies are needed. Possible challenges related to application in vivo include inhomogeneous tissues and more complex anatomies. Both HAS and MPF can readily handle inhomogeneous tissues, but the effect on MRTI accuracy will have to be investigated. For in vivo applications, blood perfusion will also have to be taken into account, and even though perfusion is relatively low in the brain and (resting) muscle tissue, it can act as a significant heat sink in tumor tissue, near blood vessels, and in organs such as the kidney. For transcranial applications, the challenge of more complex composition of the skull bone will also have to be accounted for (the plastic skull used in this study had varying thickness, but real bone will also have varying density between cortical and trabecular bone). Current clinical transcranial treatments make use of high-resolution CT scans  for phase aberration correction, but recent studies have also shown that MR imaging from both ultrashort echo time (UTE) and T1-weighted pulse sequences may achieve similar accuracy [70, 71].
Despite the challenges that arise when going to in vivo applications, recent studies using HAS to simulate the US field from a clinical tcMRgFUS system inside the human skull, in addition to using these simulations for thermal modeling utilizing PBTE, show promising results [72, 73]. This, together with results presented here and previously  showing that the MPF method is relatively robust to errors in Q and k makes it seem possible that accurate in vivo temperature maps can be achieved.
This proof-of-concept study has shown that it is possible to utilize US simulations from HAS to accurately estimate the US power density deposition Q, and further use the estimated Q to reconstruct high spatio-temporal resolution temperature maps covering large FOVs. The reconstructed temperature maps were accurate for both a homogeneous phantom and when focusing the US through a plastic skull. Future studies will aim at improving the HAS simulations through inclusion of US scattering and in vivo validation of the described methods.
Balanced steady state free precession
Echo planar imaging
Echo train length
Field of view
Full width at half maximum
Gradient recalled echo
Hybrid angular spectrum
Model predictive filtering
Magnetic resonance imaging
Magnetic resonance temperature imaging
Pennes bioheat transfer equation
Proton-resonance frequency shift
- tacq :
Ultrashort echo time
Rieke V, Butts Pauly K. MR thermometry. J Magn Reson Imaging. 2008;27:376–90.
Quesson B, de Zwart JA, Moonen CT. Magnetic resonance temperature imaging for guidance of thermotherapy. J Magn Reson Imaging. 2000;12:525–33.
Hindman JC. Proton resonance shift of water in the gas and liquid states. J Chem Phys. 1966;44:4582–92.
Kuroda K, Miki Y, Nakagawa N, Tsutsumi S, Ishihara Y, Suzuki Y, Sato K. Non-invasive temperature measurement by means of NMR parameters—use of proton chemical shift with spectral estimation technique. Med Biol Eng Comput. 1991;29:902.
Ishihara Y, Calderon A, Watanabe H, Okamoto K, Suzuki Y, Kuroda K. A precise and fast temperature mapping using water proton chemical shift. Magn Reson Med. 1995;34:814–23.
De Poorter J, De Wagter C, De Deene Y, Thomsen C, Ståhlberg F, Achten E. Noninvasive MRI thermometry with the proton resonance frequency (PRF) method: in vivo results in human muscle. Magn Reson Med. 1995;33:74–81.
McDannold N. Quantitative MRI-based temperature mapping based on the proton resonant frequency shift: review of validation studies. Int J Hyperth. 2005;21:533–46.
Todd N, Vyas U, de Bever J, Payne A, Parker DL. The effects of spatial sampling choices on MR temperature measurements. Magn Reson Med. 2011;65:515–21.
Weidensteiner C, Quesson B, Caire-Gana B, Kerioui N, Rullier A, Trillaud H, Moonen CTW. Real-time MR temperature mapping of rabbit liver in vivo during thermal ablation. Magn Reson Med. 2003;50:322–30.
Stafford RJ, Price RE, Diederich CJ, Kangasniemi M, Olsson LE, Hazle JD. Interleaved echo-planar imaging for fast multiplanar magnetic resonance temperature imaging of ultrasound thermal ablation therapy. J Magn Reson Imaging. 2004;20:706–14.
Köhler MO, Mougenot C, Quesson B, Enholm J, Le Bail B, Laurent C, Moonen CTW, Ehnholm GJ. Volumetric HIFU ablation under 3D guidance of rapid MRI thermometry. Med Phys. 2009;36:3521–35.
Kickhefel A, Roland J, Weiss C, Schick F. Accuracy of real-time MR temperature mapping in the brain: a comparison of fast sequences. Phys Med. 2010;26:192–201.
Roujol S, Ries M, Quesson B, Moonen C, Denis de Senneville B. Real-time MR-thermometry and dosimetry for interventional guidance on abdominal organs. Magn Reson Med. 2010;63:1080–7.
Quesson B, Laurent C, Maclair G, de Senneville BD, Mougenot C, Ries M, Carteret T, Rullier A, Moonen CTW. Real-time volumetric MRI thermometry of focused ultrasound ablation in vivo: a feasibility study in pig liver and kidney. NMR Biomed. 2011;24:145–53.
Todd N, Vyas U, de Bever J, Payne A, Parker DL. Reconstruction of fully three-dimensional high spatial and temporal resolution MR temperature maps for retrospective applications. Magn Reson Med. 2012;67:724–30.
Scheffler K. Fast frequency mapping with balanced SSFP: theory and application to proton-resonance frequency shift thermometry. Magn Reson Med. 2004;51:1205–11.
Paliwal V, El-Sharkawy A-M, Du X, Yang X, Atalar E. SSFP-based MR thermometry. Magn Reson Med. 2004;52:704–8.
Stafford RJ, Hazle JD, Glover GH. Monitoring of high-intensity focused ultrasound-induced temperature changes in vitro using an interleaved spiral acquisition. Magn Reson Med. 2000;43:909–12.
Fielden S, Zhao L, Miller W, Feng X, Wintermark M, Pauly KB, Meyer C. Accelerating 3D spiral MR thermometry with Kalman filter. In: Int. Soc. Magn. Reson. Med. Milan: 2014, 22:2346.
Marx M, Plata J, Butts Pauly K. Toward volumetric MR thermometry with the MASTER sequence. IEEE Trans Med Imaging. 2014;62:1–9.
Vallo S, Eichler K, Kelly K, Schulz B, Bartsch G, Haferkamp A, Vogl TJ, Zangos S. MR-guided laser-induced thermotherapy in ex vivo porcine kidney: comparison of four different imaging sequences. Lasers Surg Med. 2014;46:558–62.
Weidensteiner C, Kerioui N, Quesson B, Denis de Senneville B, Trillaud H, Moonen CTW. Stability of real-time MR temperature mapping in healthy and diseased human liver. J Magn Reson Imaging. 2004;19:438–46.
Bankson JA, Stafford RJ, Hazle JD. Partially parallel imaging with phase-sensitive data: increased temporal resolution for magnetic resonance temperature imaging. Magn Reson Med. 2005;53:658–65.
Guo J-Y, Kholmovski EG, Zhang L, Jeong E-K, Parker DL. k-space inherited parallel acquisition (KIPA): application on dynamic magnetic resonance imaging thermometry. Magn Reson Imaging. 2006;24:903–15.
Streicher MN, Schäfer A, Müller D, Kögler C, Reimer E, Dhital B, Trampel R, Rivera D, Pampel A, Ivanov D, Turner R. Frequency-selective asymmetric spin-echo EPI with parallel imaging for fast internally referenced MR thermometry. Int Soc Magn Reson Med. 2011;19:529.
Mei CS, Panych LP, Yuan J, McDannold NJ, Treat LH, Jing Y, Madore B. Combining two-dimensional spatially selective RF excitation, parallel imaging, and UNFOLD for accelerated MR thermometry imaging. Magn Reson Med. 2011;66:112–22.
Leonard P, Chopra R, Nachman A. Compressed sensing for accelerated MR thermometry in MRI-controlled transurethral ultrasound therapy. In: Proc. 20th Sci. Meet. ISMRM, Melb. Melbourne: 2012:2918.
Marx M, Butts Pauly K. Use of Compressed sensing for acceleration of volumetric MR thermometry. In: 9th Int. Interv. MRI Symp. Boston: 2012, c:110.
Cao Z, Oh S, Otazo R, Sica CT, Griswold MA, Collins CM. Complex difference constrained compressed sensing reconstruction for accelerated PRF thermometry with application to MRI-induced RF heating. Magn Reson Med. 2015;73:1420–31.
Todd N, Adluru G, Payne A, DiBella EVR, Parker D. Temporally constrained reconstruction applied to MRI temperature data. Magn Reson Med. 2009;62:406–19.
Todd N, Prakash J, Odéen H, de Bever J, Payne A, Yalavarthy P, Parker DL. Toward real-time availability of 3D temperature maps created with temporally constrained reconstruction. Magn Reson Med. 2014;71:1394–404.
Todd N, Payne A, Parker DL. Model predictive filtering for improved temporal resolution in MRI temperature imaging. Magn Reson Med. 2010;63:1269–79.
Roujol S, Denis de Senneville B, Hey S, Moonen C, Ries M. Robust adaptive extended Kalman filtering for real time MR-thermometry guided HIFU interventions. IEEE Trans Med Imaging. 2012;31:533–42.
Fuentes D, Yung J, Hazle JD, Weinberg JS, Stafford RJ. Kalman filtered MR temperature imaging for laser induced thermal therapies. IEEE Trans Med Imaging. 2012;31:984–94.
Gaur P, Grissom WA. Accelerated MRI thermometry by direct estimation of temperature from undersampled k-space data. Magn Reson Med. 2015;73:1914–25.
Dillon CR, Vyas U, Payne A, Christensen DA, Roemer RB. An analytical solution for improved HIFU SAR estimation. Phys Med Biol. 2012;57:4527–44.
Dillon CR, Payne A, Christensen DA, Roemer RB. The accuracy and precision of two non-invasive, magnetic resonance-guided focused ultrasound-based thermal diffusivity estimation methods. Int J Hyperth. 2014;30:362–71.
McDannold NJ, King RL, Jolesz F a, Hynynen KH. Usefulness of MR imaging-derived thermometry and dosimetry in determining the threshold for tissue damage induced by thermal surgery in rabbits. Radiology. 2000;216:517–23.
McDannold N, Vykhodtseva N, Jolesz FA, Hynynen K. MRI investigation of the threshold for thermally induced blood-brain barrier disruption and brain tissue damage in the rabbit brain. Magn Reson Med. 2004;51:913–23.
Paeng D-G, Xu Z, Snell J, Quigg AH, Eames M, Jin C, Everstine AC, Sheehan JP, Lopes BS, Kassell N. Thermal dose effects by MR- guided focused ultrasound on the pig brain tissue—preliminary result. In: Int Symp Ther Ultrasound. 2016;61.
Vyas U, Christensen DA. Ultrasound beam propagation using the hybrid angular spectrum method. In: Proc. IEEE Eng. Med. Biol. Soc. Vancouver: 2008:2526–2529.
Vyas U, Christensen DA. Extension of the angular spectrum method to calculate pressure from a spherically curved acoustic source. J Acoust Soc Am. 2011;130:2687.
Vyas U, Christensen DA. Ultrasound beam simulations in inhomogeneous tissue geometries using the hybrid angular spectrum method. IEEE Trans Ultrason Ferroelectr Freq Control. 2012;59:1093–100.
Ratcliffe JA. Some aspects of diffraction theory and their application to the ionosphere. Rep Prog Phys. 1956;19:188–267.
Christopher PT, Parker KJ. New approaches to the linear propagation of acoustic fields. J Acoust Soc Am. 1991;90:507–21.
Goodman JW. Introduction to Fourier optics. 3rd ed. Greenwood Village: Roberts and Co.; 2004.
Odéen H, Todd N, Dillon C, Payne A, Parker DL. Model predictive filtering MR thermometry: Effects of model inaccuracies, k-space reduction factor, and temperature increase rate. Magn Reson Med. 2015, ePub ahead.
Odéen H, Roberts J, de Bever J, Parker DL. Real-time online reconstruction of 3D MR thermometry data for MRgFUS applications. In: Int. Soc. Magn. Reson. Med. Toronto: 2015, 1279:3716.
Almquist S, de Bever J, Merrill R, Parker D, Christensen D. A full-wave phase aberration correction method for transcranial high-intensity focused ultrasound brain therapies. In: IEEE Eng Med Biol Soc. 2014, 2014:310–13.
Odéen H, Todd N, Diakite M, Minalga E, Payne A, Parker DL. Sampling strategies for subsampled segmented EPI PRF thermometry in MR guided high intensity focused ultrasound. Med Phys. 2014;41:092301.
Du YP, Parker DL, Davis WL, Cao G. Reduction of partial-volume artifacts with zero-filled interpolation in three-dimensional MR angiography. J Magn Reson Imaging. 1994;4:733–41.
Adluru G, Awate SP, Tasdizen T, Whitaker RT, Dibella EVR. Temporally constrained reconstruction of dynamic cardiac perfusion MRI. Magn Reson Med. 2007;57:1027–36.
Pennes HH. Analysis of tissue and arterial blood temperatures in the resting human forearm. J Appl Physiol. 1948;1:93–122.
Farrer AI, Odéen H, de Bever J, Coats B, Parker DL, Payne A, Christensen DA. Characterization and evaluation of tissue-mimicking gelatin phantoms for use with MRgFUS. J Ther Ultrasound. 2015;3:9.
Dillon CR, Todd N, Payne A, Parker DL, Christensen DA, Roemer RB. Effects of MRTI sampling characteristics on estimation of HIFU SAR and tissue thermal diffusivity. Phys Med Biol. 2013;58:7291–307.
Mark JE. Physical properties of polymers handbook. New York: Springer Science & Business Media; 2007.
The Engineering Toolbox. [http://www.engineeringtoolbox.com/]. Accessed 05 Oct 2015.
Shi YC, Parker DL, Dillon CR. Sensitivity of tissue properties derived from MRgFUS temperature data to input errors and data inclusion criteria: ex vivo study in porcine muscle. Phys Med Biol. 2016;61:N373–85.
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.
Ghanouni P, Pauly KB, Elias WJ, Henderson J, Sheehan J, Monteith S, Wintermark M. Transcranial MRI-guided focused ultrasound: a review of the technologic and neurologic applications. Am J Roentgenol. 2015;205:150–9.
McDannold N, Maier SE. Magnetic resonance acoustic radiation force imaging. Med Phys. 2008;35:3748.
Bitton RR, Kaye E, Dirbas FM, Daniel BL, Pauly KB. Toward MR-guided high intensity focused ultrasound for presurgical localization: focused ultrasound lesions in cadaveric breast tissue. J Magn Reson Imaging. 2012;35:1089–97.
Kaye EA, Pauly KB. Adapting MRI acoustic radiation force imaging for in vivo human brain focused ultrasound applications. Magn Reson Med. 2013;69:724–33.
de Bever JT, Odéen H, Todd N, Farrer AI, Parker DL. Evaluation of a three-dimensional MR acoustic radiation force imaging pulse sequence using a novel unbalanced bipolar motion encoding gradient. Magn Reson Med. 2016;76:803–13.
Auboiroux V, Viallon M, Roland J, Hyacinthe JN, Petrusca L, Morel DR, Goget T, Terraz S, Gross P, Becker CD, Salomir R. ARFI-prepared MRgHIFU in liver: simultaneous mapping of ARFI-displacement and temperature elevation, using a fast GRE-EPI sequence. Magn Reson Med. 2012;68:932–46.
de Bever JT, Odéen H, Todd N, Farrer AI, Parker DL. Evaluation of a three-dimensional MR acoustic radiation force imaging pulse sequence using a novel unbalanced bipolar motion encoding gradient. Magn Reson Med. 2015. ePub ahead of print.
Coutts G, Gilderdale DJ, Chui M, Kasuboski L, DeSouza NM. Integrated and interactive position tracking and imaging of interventional tools and internal devices using small fiducial receiver coils. Magn Reson Med. 1998;40:908–13.
Svedin BT, Beck MJ, Hadley JR, Merrill R, Bolster BD, Parker DL. Focal position determination in breast MRgHIFU using 3 tracking coils. In: Proc. Int. Soc. Magn. Reson. Med. Toronto: 2015:1641.
Odéen H, de Bever J, Almquist S, Farrer A, Todd N, Payne A, Snell JW, Christensen DA, Parker DL. Treatment envelope evaluation in transcranial magnetic resonance-guided focused ultrasound utilizing 3D MR thermometry. J Ther Ultrasound. 2014;2:19.
Miller GW, Eames M, Snell J, Aubry J-F. Ultrashort echo-time MRI versus CT for skull aberration correction in MR-guided transcranial focused ultrasound: in vitro comparison on human calvaria. Med Phys. 2015;42:2223–33.
Wintermark M, Tustison NJ, Elias WJ, Patrie JT, Xin W, Demartini N, Eames M, Sumer S, Lau B, Cupino A, Snell J, Hananel A, Kassell N, Aubry J-F. T1-weighted MRI as a substitute to CT for refocusing planning in MR-guided focused ultrasound. Phys Med Biol. 2014;59:3599–614.
Vyas U, Webb T, Bitton R, Pauly B, Ghanouni P. Acoustic and thermal simulations of tcMRgFUS in patient specific models: validation with experiments. J Ther Ultrasound. 2015;3:P35.
Almquist S, Parker D, Christensen D. Simulation of hemispherical transducers for transcranial HIFU treatments using the hybrid angular spectrum approach. J Ther Ultrasound. 2015;3:P9.
The authors appreciate helpful contributions from collaborators at the University of Utah.
This work was supported by The Focused Ultrasound Surgery Foundation, Siemens Healthcare, The Ben B. and Iris M. Margolis Foundation, and NIH grants R01s EB013433, and CA134599.
Availability of data and material
Please contact author for data requests.
All authors helped design the study. HO and SA performed the MRgFUS experiments under DLP’s supervision. HO conducted the data reconstruction and analysis, and drafted the manuscript. SA and DAC developed the ultrasound modeling and phase aberration correction, which was performed by SA. JdB designed, built, and implemented hardware and software to trigger the FUS from the MR scanner and to apply the phase aberration correction. All authors read and approved the manuscript.
The authors declare that they have no competing interests.
Consent for publication
Ethics approval and consent to participate