 Research
 Open Access
 Published:
MR thermometry for focused ultrasound monitoring utilizing model predictive filtering and ultrasound beam modeling
Journal of Therapeutic Ultrasoundvolume 4, Article number: 23 (2016)
Abstract
Background
A major challenge in using magnetic resonance temperature imaging (MRTI) to monitor focused ultrasound (FUS) applications is achieving high spatiotemporal 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 modelbased 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 pretreatment, lowpower sonication.
Methods
In this proofofconcept 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 pretreatment sonication and any potential tissue damage. MRTI reconstructions are performed using a thermal modelbased 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 rootmeansquare errors (RMSE) and focal spot positions and shapes are evaluated.
Results
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.
Conclusions
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 modelbased reconstruction to obtain accurate temperature maps with high spatiotemporal resolution over large FOVs.
Background
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 realtime 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 protonresonance frequency shift (PRFS) method, which relies on the temperaturedependent 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 aqueousbased tissues [7]. One drawback of the PRFS method is that it does not work in lipidbased 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 [8] 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 kspace lines are acquired after each radiofrequency (RF) excitation. Both single and multislice 2D as well as fully 3D versions, with both singleshot and segmented phase encodings, have been investigated [9–15]. Other fast pulse sequences that have been investigated include balanced steadystate free precession (bSSFP)type sequences which have a relatively short repetition time (TR) [16, 17], as well as utilizing signaltonoise ratio (SNR)efficient nonCartesian approaches [18–21]. Speedups have also been achieved by combining subsampled kspace 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 kspace 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 modelbased 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 kspace from the subsampled acquisition. MPF can achieve high spatiotemporal 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/m^{3}), can be estimated from a pretreatment lowtemperature 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 pretreatment heating and can result in shorter total treatment time since the treatment can be planned and Q can be calculated beforehand. In this proofofconcept study, experiments are performed in a homogeneous gelatin phantom and in a gelatin phantom with an embedded plastic skull. In the skull phantom, HASbased US simulations are also used for phase aberration correction using a previously described method [49].
Methods
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. kspace was subsampled with a Cartesian variabledensity equally spaced sampling scheme [50], where the sampling in the kzslice encode direction was varied to achieve denser sampling over the kspace center while maintaining a constant echospacing 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 zerofilled interpolated to 0.5mm isotropic voxel spacing to minimize partial volume effects [51].
Two inhousebuilt RF receiveonly singlechannel 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 sensinglike 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 [53] given by the following:
where ρ = tissue density (kg/m^{3}); 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/m^{3}/s); Q = US power deposition density (W/m^{3}). 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 kspace and used to fill in the parts of kspace that were not sampled during time frame (n + 1). This updated kspace 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 kspace data, and the constraint term applies an appropriate temporal smoothness constraint, in this work utilizing the L1norm.
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 postprocessing 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) [54]; (2) with MPF using Q and k determined from fitting a pretreatment, lowpower sonication (i.e., lowtemperature rise) using a previously published analytical solution of the heating [36, 37, 55] together with the previously published phantom properties (ρ, C _{ t } ) [54]; 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 kspace data over a smaller FOV (Table 1). The shortduration, lowtemperature 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).
Phantoms
Both phantoms were fabricated from 250bloom gelatin [54] [speed of sound = 1553 ± 10 m/s (mean ± standard deviation), US attenuation = 0.06 ± 0.01 Np/(cm*MHz), density ρ = 1057 ± 44 kg/m^{3}, specific heat capacity C _{ t } = 3635 ± 88 J/(kg*°C), thermal conductivity k = 0.55 ± 0.01 W/(m*°C)] using previously published methods [54]. 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/m^{3} [56, 57]].
FUS sonications
FUS sonications were performed using an MRcompatible phasedarray ultrasound transducer (256 elements, 1MHz frequency, 13cm radius of curvature, focal spot size full width at half maximum 2 × 2 × 8 mm, Imasonic, Voraysurl’Ognon, France), accompanying hardware and software for mechanical positioning and electronic beam steering (Image Guided Therapy, Pessac, France), and inhouse 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 deionized and degassed water as indicated in Fig. 1.
In the homogeneous phantom, three continuouswave 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 lowtemperature rise continuouswave 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 continuouswave sonications each for fully sampled truth and subsampled kspace, 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 fullwave simulation method that takes into account US refraction, absorption, and firstorder reflection. Inhomogeneous tissues are modeled by segmenting the anatomy into a 3D Cartesian grid and assigning unique speedofsound, absorption coefficient, and density values for the tissue type within each voxel. The hybrid angular spectrum approach leapfrogs between the space and spatialfrequency domains, utilizing the space domain to account for individual voxel attenuation and speedofsound properties and the spatialfrequency domain to linearly propagate the US wave from planetoplane 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., finitedifference timedomain (FDTD) simulations [49].
Due to differences in the speed of sound in the skull and gelatin, the US waves from the different elements in the phasedarray 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 phasedarray transducer at the intended focal spot location, with an initial assumption of all elements having zero phase and the same amplitude [49]. 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 [58]). Since HAS is a fullwave 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 twocompartment 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 Qpattern inside the phantom. For the skull phantom study, the segmentation was based on a highresolution 3D GRE scan (i.e., no EPI readout) and a threecompartment segmentation consisting of water, phantom material (gelatin), and skull was performed.
Evaluation of accuracy and precision
The rootmeansquare 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 kspace 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 kspace. 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 zerofilled interpolated voxels.
Qpatterns as obtained from the HAS simulations and the analytical fits to the LP pretreatment 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.
Results
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 inplane (i.e., in the FEPEplane) 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 throughplane (i.e., in the SEdirection) 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 inplane (MPF HAS in FE) and 0.3 mm throughplane (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 Qvoxel for HAS prediction and LP analytical fit. For most cases, the focal spot positions coincide, and in the cases when the distance is nonzero, 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 Qpatterns agree to within the measured accuracy, but HAS simulated Q is approximately 1 mm longer in the throughplane direction. In the skull phantom, the differences are greater with the HASsimulated Q being approximately 1 mm narrower inplane and 1.2 mm shorter throughplane.
For all MPF reconstructions using HASsimulated Q, a tabular value of k = 5.49 × 10^{−4} W/(m*°C) was used for thermal conductivity [54], 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 noncorrected 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.
Discussion
In this proofofconcept 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 HASsimulated 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 HASbased temperature maps agreed well with fully sampled truth, and were all within the size of one acquisition voxel (before zerofilling). These results demonstrate the feasibility of using US simulations in conjunction with modelbased reconstruction to achieve high spatiotemporal resolution, largeFOV 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 Qpatterns (Table 2), and their locations coincide with within the zerofilled 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 Gaussianshaped focal spot, the comparison and error calculations were focused on the focal spot (i.e., not including offfocal regions). Although no truth for the US near and farfield (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 kspace. This ensures that no offfocal heating will go undetected, and as long as any offfocal 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 offfocal regions, so that no heat dissipation is modeled. It can further be noted that HAS predicts power deposition for offfocal areas which could be used for MPF reconstructions of, e.g., the US nearfield. The accuracy and precision of this will be investigated in future studies.
The accuracy of the position of the Qpatterns determined with HAS makes it a promising approach for pretreatment localization of the focal spot. In many current applications, such as transcranial treatments, initial lowpower 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 tissuedamaging pretreatment 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 MRARFI, lowenergy 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 MRARFI), tracker coils will not work when focusing through an aberrating medium such as the skull bone. US simulationbased focal spot localization hence has the potential to both avoid the pretreatment sonication and work in, e.g., transcranial applications.
The fact that the magnitude of the Qpatterns 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 kspace subsampling (R) and the rate of temperature increase (dT/dt) [47]. 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 HASsimulated 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 [69].
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 highresolution computer tomography (CT) images with submillimeter resolution [60], 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 highresolution CT scans [60] for phase aberration correction, but recent studies have also shown that MR imaging from both ultrashort echo time (UTE) and T1weighted 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 [47] 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.
Conclusions
This proofofconcept 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 spatiotemporal 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.
Abbreviations
 3D:

Three dimensional
 bSSFP:

Balanced steady state free precession
 BW:

Bandwidth
 EPI:

Echo planar imaging
 ES:

Echo spacing
 ETL:

Echo train length
 FA:

Flip angle
 FDTD:

Finitedifference timedomain
 FE:

Frequency encoding
 FOV:

Field of view
 FUS:

Focused ultrasound
 FWHM:

Full width at half maximum
 GRE:

Gradient recalled echo
 HAS:

Hybrid angular spectrum
 HV:

Hottest voxel
 LP:

Low power
 MPF:

Model predictive filtering
 MR:

Magnetic resonance
 MRI:

Magnetic resonance imaging
 MRTI:

Magnetic resonance temperature imaging
 PBTE:

Pennes bioheat transfer equation
 PE:

Phase encoding
 PRFS:

Protonresonance frequency shift
 R:

Reduction factor
 RF:

Radiofrequency
 RMSE:

Rootmeansquare error
 SD:

Standard deviation
 SE:

Slice encoding
 SNR:

Signaltonoise ratio
 t_{acq} :

Acquisition time
 TE:

Echo time
 TR:

Repetition time
 US:

Ultrasound
 UTE:

Ultrashort echo time
References
 1.
Rieke V, Butts Pauly K. MR thermometry. J Magn Reson Imaging. 2008;27:376–90.
 2.
Quesson B, de Zwart JA, Moonen CT. Magnetic resonance temperature imaging for guidance of thermotherapy. J Magn Reson Imaging. 2000;12:525–33.
 3.
Hindman JC. Proton resonance shift of water in the gas and liquid states. J Chem Phys. 1966;44:4582–92.
 4.
Kuroda K, Miki Y, Nakagawa N, Tsutsumi S, Ishihara Y, Suzuki Y, Sato K. Noninvasive temperature measurement by means of NMR parameters—use of proton chemical shift with spectral estimation technique. Med Biol Eng Comput. 1991;29:902.
 5.
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.
 6.
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.
 7.
McDannold N. Quantitative MRIbased temperature mapping based on the proton resonant frequency shift: review of validation studies. Int J Hyperth. 2005;21:533–46.
 8.
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.
 9.
Weidensteiner C, Quesson B, CaireGana B, Kerioui N, Rullier A, Trillaud H, Moonen CTW. Realtime MR temperature mapping of rabbit liver in vivo during thermal ablation. Magn Reson Med. 2003;50:322–30.
 10.
Stafford RJ, Price RE, Diederich CJ, Kangasniemi M, Olsson LE, Hazle JD. Interleaved echoplanar imaging for fast multiplanar magnetic resonance temperature imaging of ultrasound thermal ablation therapy. J Magn Reson Imaging. 2004;20:706–14.
 11.
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.
 12.
Kickhefel A, Roland J, Weiss C, Schick F. Accuracy of realtime MR temperature mapping in the brain: a comparison of fast sequences. Phys Med. 2010;26:192–201.
 13.
Roujol S, Ries M, Quesson B, Moonen C, Denis de Senneville B. Realtime MRthermometry and dosimetry for interventional guidance on abdominal organs. Magn Reson Med. 2010;63:1080–7.
 14.
Quesson B, Laurent C, Maclair G, de Senneville BD, Mougenot C, Ries M, Carteret T, Rullier A, Moonen CTW. Realtime volumetric MRI thermometry of focused ultrasound ablation in vivo: a feasibility study in pig liver and kidney. NMR Biomed. 2011;24:145–53.
 15.
Todd N, Vyas U, de Bever J, Payne A, Parker DL. Reconstruction of fully threedimensional high spatial and temporal resolution MR temperature maps for retrospective applications. Magn Reson Med. 2012;67:724–30.
 16.
Scheffler K. Fast frequency mapping with balanced SSFP: theory and application to protonresonance frequency shift thermometry. Magn Reson Med. 2004;51:1205–11.
 17.
Paliwal V, ElSharkawy AM, Du X, Yang X, Atalar E. SSFPbased MR thermometry. Magn Reson Med. 2004;52:704–8.
 18.
Stafford RJ, Hazle JD, Glover GH. Monitoring of highintensity focused ultrasoundinduced temperature changes in vitro using an interleaved spiral acquisition. Magn Reson Med. 2000;43:909–12.
 19.
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.
 20.
Marx M, Plata J, Butts Pauly K. Toward volumetric MR thermometry with the MASTER sequence. IEEE Trans Med Imaging. 2014;62:1–9.
 21.
Vallo S, Eichler K, Kelly K, Schulz B, Bartsch G, Haferkamp A, Vogl TJ, Zangos S. MRguided laserinduced thermotherapy in ex vivo porcine kidney: comparison of four different imaging sequences. Lasers Surg Med. 2014;46:558–62.
 22.
Weidensteiner C, Kerioui N, Quesson B, Denis de Senneville B, Trillaud H, Moonen CTW. Stability of realtime MR temperature mapping in healthy and diseased human liver. J Magn Reson Imaging. 2004;19:438–46.
 23.
Bankson JA, Stafford RJ, Hazle JD. Partially parallel imaging with phasesensitive data: increased temporal resolution for magnetic resonance temperature imaging. Magn Reson Med. 2005;53:658–65.
 24.
Guo JY, Kholmovski EG, Zhang L, Jeong EK, Parker DL. kspace inherited parallel acquisition (KIPA): application on dynamic magnetic resonance imaging thermometry. Magn Reson Imaging. 2006;24:903–15.
 25.
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. Frequencyselective asymmetric spinecho EPI with parallel imaging for fast internally referenced MR thermometry. Int Soc Magn Reson Med. 2011;19:529.
 26.
Mei CS, Panych LP, Yuan J, McDannold NJ, Treat LH, Jing Y, Madore B. Combining twodimensional spatially selective RF excitation, parallel imaging, and UNFOLD for accelerated MR thermometry imaging. Magn Reson Med. 2011;66:112–22.
 27.
Leonard P, Chopra R, Nachman A. Compressed sensing for accelerated MR thermometry in MRIcontrolled transurethral ultrasound therapy. In: Proc. 20th Sci. Meet. ISMRM, Melb. Melbourne: 2012:2918.
 28.
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.
 29.
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 MRIinduced RF heating. Magn Reson Med. 2015;73:1420–31.
 30.
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.
 31.
Todd N, Prakash J, Odéen H, de Bever J, Payne A, Yalavarthy P, Parker DL. Toward realtime availability of 3D temperature maps created with temporally constrained reconstruction. Magn Reson Med. 2014;71:1394–404.
 32.
Todd N, Payne A, Parker DL. Model predictive filtering for improved temporal resolution in MRI temperature imaging. Magn Reson Med. 2010;63:1269–79.
 33.
Roujol S, Denis de Senneville B, Hey S, Moonen C, Ries M. Robust adaptive extended Kalman filtering for real time MRthermometry guided HIFU interventions. IEEE Trans Med Imaging. 2012;31:533–42.
 34.
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.
 35.
Gaur P, Grissom WA. Accelerated MRI thermometry by direct estimation of temperature from undersampled kspace data. Magn Reson Med. 2015;73:1914–25.
 36.
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.
 37.
Dillon CR, Payne A, Christensen DA, Roemer RB. The accuracy and precision of two noninvasive, magnetic resonanceguided focused ultrasoundbased thermal diffusivity estimation methods. Int J Hyperth. 2014;30:362–71.
 38.
McDannold NJ, King RL, Jolesz F a, Hynynen KH. Usefulness of MR imagingderived thermometry and dosimetry in determining the threshold for tissue damage induced by thermal surgery in rabbits. Radiology. 2000;216:517–23.
 39.
McDannold N, Vykhodtseva N, Jolesz FA, Hynynen K. MRI investigation of the threshold for thermally induced bloodbrain barrier disruption and brain tissue damage in the rabbit brain. Magn Reson Med. 2004;51:913–23.
 40.
Paeng DG, 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.
 41.
Vyas U, Christensen DA. Ultrasound beam propagation using the hybrid angular spectrum method. In: Proc. IEEE Eng. Med. Biol. Soc. Vancouver: 2008:2526–2529.
 42.
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.
 43.
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.
 44.
Ratcliffe JA. Some aspects of diffraction theory and their application to the ionosphere. Rep Prog Phys. 1956;19:188–267.
 45.
Christopher PT, Parker KJ. New approaches to the linear propagation of acoustic fields. J Acoust Soc Am. 1991;90:507–21.
 46.
Goodman JW. Introduction to Fourier optics. 3rd ed. Greenwood Village: Roberts and Co.; 2004.
 47.
Odéen H, Todd N, Dillon C, Payne A, Parker DL. Model predictive filtering MR thermometry: Effects of model inaccuracies, kspace reduction factor, and temperature increase rate. Magn Reson Med. 2015, ePub ahead.
 48.
Odéen H, Roberts J, de Bever J, Parker DL. Realtime online reconstruction of 3D MR thermometry data for MRgFUS applications. In: Int. Soc. Magn. Reson. Med. Toronto: 2015, 1279:3716.
 49.
Almquist S, de Bever J, Merrill R, Parker D, Christensen D. A fullwave phase aberration correction method for transcranial highintensity focused ultrasound brain therapies. In: IEEE Eng Med Biol Soc. 2014, 2014:310–13.
 50.
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.
 51.
Du YP, Parker DL, Davis WL, Cao G. Reduction of partialvolume artifacts with zerofilled interpolation in threedimensional MR angiography. J Magn Reson Imaging. 1994;4:733–41.
 52.
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.
 53.
Pennes HH. Analysis of tissue and arterial blood temperatures in the resting human forearm. J Appl Physiol. 1948;1:93–122.
 54.
Farrer AI, Odéen H, de Bever J, Coats B, Parker DL, Payne A, Christensen DA. Characterization and evaluation of tissuemimicking gelatin phantoms for use with MRgFUS. J Ther Ultrasound. 2015;3:9.
 55.
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.
 56.
Mark JE. Physical properties of polymers handbook. New York: Springer Science & Business Media; 2007.
 57.
The Engineering Toolbox. [http://www.engineeringtoolbox.com/]. Accessed 05 Oct 2015.
 58.
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.
 59.
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.
 60.
Ghanouni P, Pauly KB, Elias WJ, Henderson J, Sheehan J, Monteith S, Wintermark M. Transcranial MRIguided focused ultrasound: a review of the technologic and neurologic applications. Am J Roentgenol. 2015;205:150–9.
 61.
McDannold N, Maier SE. Magnetic resonance acoustic radiation force imaging. Med Phys. 2008;35:3748.
 62.
Bitton RR, Kaye E, Dirbas FM, Daniel BL, Pauly KB. Toward MRguided high intensity focused ultrasound for presurgical localization: focused ultrasound lesions in cadaveric breast tissue. J Magn Reson Imaging. 2012;35:1089–97.
 63.
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.
 64.
de Bever JT, Odéen H, Todd N, Farrer AI, Parker DL. Evaluation of a threedimensional MR acoustic radiation force imaging pulse sequence using a novel unbalanced bipolar motion encoding gradient. Magn Reson Med. 2016;76:803–13.
 65.
Auboiroux V, Viallon M, Roland J, Hyacinthe JN, Petrusca L, Morel DR, Goget T, Terraz S, Gross P, Becker CD, Salomir R. ARFIprepared MRgHIFU in liver: simultaneous mapping of ARFIdisplacement and temperature elevation, using a fast GREEPI sequence. Magn Reson Med. 2012;68:932–46.
 66.
de Bever JT, Odéen H, Todd N, Farrer AI, Parker DL. Evaluation of a threedimensional MR acoustic radiation force imaging pulse sequence using a novel unbalanced bipolar motion encoding gradient. Magn Reson Med. 2015. ePub ahead of print.
 67.
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.
 68.
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.
 69.
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 resonanceguided focused ultrasound utilizing 3D MR thermometry. J Ther Ultrasound. 2014;2:19.
 70.
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.
 71.
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 JF. T1weighted MRI as a substitute to CT for refocusing planning in MRguided focused ultrasound. Phys Med Biol. 2014;59:3599–614.
 72.
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.
 73.
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.
Acknowledgements
The authors appreciate helpful contributions from collaborators at the University of Utah.
Funding
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.
Authors’ contributions
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.
Competing interests
The authors declare that they have no competing interests.
Consent for publication
Not applicable
Ethics approval and consent to participate
Not applicable
Author information
Rights and permissions
Open Access This 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.
About this article
Received
Accepted
Published
DOI
Keywords
 MR thermometry
 Ultrasound modeling
 Model predictive filtering
 Hybrid angular spectrum
 Protonresonance frequency shift