Compensating for Magnetic Field Inhomogeneity in
Multigradient-Echo-Based MR Thermometry
Frank F.J. Simonis,
1* Esben T. Petersen,
1Lambertus W. Bartels,
2Jan J. W. Lagendijk,
1and Cornelis A.T. van den Berg
1Purpose: MR thermometry (MRT) is a noninvasive method for measuring temperature that can potentially be used for radio frequency (RF) safety monitoring. This application requires measuring absolute temperature. In this study, a multigradient-echo (mGE) MRT sequence was used for that purpose. A drawback of this sequence, however, is that its accuracy is affected by background gradients. In this article, we present a method to minimize this effect and to improve absolute tem-perature measurements using MRI.
Theory: By determining background gradients using a B0map
or by combining data acquired with two opposing readout directions, the error can be removed in a homogenous phan-tom, thus improving temperature maps.
Methods: All scans were performed on a 3T system using eth-ylene glycol-filled phantoms. Background gradients were var-ied, and one phantom was uniformly heated to validate both compensation approaches. Independent temperature record-ings were made with optical probes.
Results: Errors correlated closely to the background gradients in all experiments. Temperature distributions showed a much smaller standard deviation when the corrections were applied (0.21C vs.
0.45C) and correlated well with thermo-optical probes.
Conclusion: The corrections offer the possibility to measure RF heating in phantoms more precisely. This allows mGE MRT to become a valuable tool in RF safety assessment. Magn Reson Med 73:1184–1189, 2015.VC 2014 The Authors
Mag-netic Resonance in Medicine published by Wiley Periodi-cals, Inc. on behalf of International Society for Magnetic Resonance in Medicine. This is an open access article under the terms of the Creative Commons Attribution-NonCommercial-NoDerivs License, which permits use and distribution in any medium, provided the original work is properly cited, the use is non-commercial and no modifi-cations or adaptations are made.
Key words: MR thermometry; multigradient-echo; background gradient; absolute temperature; RF safety monitoring
INTRODUCTION
MR thermometry (MRT) is a noninvasive method for mon-itoring thermotherapy interventions, for example, high-intensity focused ultrasound, laser-induced thermother-apy, and radio frequency (RF) ablation (1). MRT can moni-tor temperature distributions inside the body, enabling MR guidance of thermal therapy. A commonly used mea-sure to quantify the therapeutic effect is thermal dose, which can be expressed in the amount of cumulative equivalent minutes at 43C (CEM
43) (2,3). Therefore,
abso-lute temperature measurements are very valuable in pre-dicting response and monitoring surrounding tissue.
MRT can also be useful for RF safety assessments in MR (4–7). Simulation studies have shown that RF deposition quantified by specific absorption rate (SAR) can spatially be very inhomogeneous, leading to hotspots that can exceed temperature limits set in regulations (8,9). To ensure safety, RF power is restricted to comply with SAR constraints as defined by International Electrotechnical Commission guidelines (10). However, tissue damage is not directly related to global or local SAR but to the duration of absolute temperature levels, similar to the thermal dose concept in thermotherapy. Thus, knowledge of absolute tis-sue temperature should lie at the core of every RF safety assessment. This has been pursued in several studies by means of the thermal modeling of the human body with assumptions on patient-specific thermal properties such as tissue perfusion and discrete vasculature (11–16). A major improvement would be to directly measure absolute tem-perature inside the patient exposed to an intense RF field.
However, there are some prerequisites for such meas-urements: The precision and accuracy should be high (60.2C), as well as the spatial resolution, in order to
detect even small local temperature increases that are constrained by the boundaries set in the guidelines (10). A method that combines all desired characteristics has not been presented to date.
There are several temperature-dependent properties in MRI, but proton resonance frequency shift (PRFS) methods are most commonly used in MRT (17). These methods are based on the principle that the nuclei of hydrogen atoms in hydrogen-bonds have a temperature-dependent electron screening. Therefore, their resonance frequency is also slightly temperature-dependent (a ¼ 0.01 ppm=C); and the
differ-ence in phase of an MR signal acquired with a gradient echo before and after heating is thus a direct measure for the change
1
Department of Radiotherapy, Imaging Division, University Medical Center Utrecht, Utrecht, The Netherlands.
2
Image Sciences Institute, Imaging Division, University Medical Center Utrecht, Utrecht, The Netherlands.
*Correspondence to: Frank F.J. Simonis, M.Sc., Department of Radiother-apy, University Medical Center Utrecht, Heidelberglaan 100, 3584 CX, Utrecht, The Netherlands. E-mail: F.F.J.Simonis@umcutrecht.nl
The copyright line for this article was changed on 1 November 2016 after original online publication.
Received 11 October 2013; revised 18 February 2014; accepted 20 February 2014
DOI 10.1002/mrm.25207
Published online 24 March 2014 in Wiley Online Library (wileyonlinelibrary.com). VC 2014 The Authors Magnetic Resonance in Medicine published by Wiley Periodicals, Inc. on behalf of International Society for Magnetic Resonance in Medicine. This is an open access article under the terms of the Creative Commons Attribution-NonCommercial-oDerivs License, which permits use and distribution in any medium, provided the original work is properly cited, the use is non-commercial and no modifications or adaptations are made.
Magnetic Resonance in Medicine 73:1184–1189 (2015)
in temperature. A large disadvantage of PRFS thermometry is that it does not provide absolute temperature—only a change in temperature relative to a baseline scan.
An alternative MRT technique is based on a multigra-dient-echo sequence (mGE) and can measure absolute temperatures (18). It relies on a signal containing a temperature-dependent frequency fT (e.g., water) and a
temperature-independent frequency of a reference com-pound fref (e.g., fat). When both compounds are present
and mixed in the same voxel, a beating in the signal enve-lope of the echo train occurs at the difference frequency (Df) of fT and fref. The reference compound makes the
method insensitive to magnet drift and movement, thus enabling calculation of an absolute temperature map.
A drawback of mGE MRT, however, is that it is ham-pered by background gradients caused by susceptibility variations and magnet imperfections (19). These can be reduced by precise shimming, but residuals will always be present—especially at higher magnetic field strengths. Background gradients lead to spatially varying, systematic errors in the temperature estimation, which make the method unsuitable for measuring small temperature changes caused by RF heating. Therefore, in this study a solution for the problems caused by background gradients is proposed and tested in phantoms at a field strength of 3T, both at a constant temperature and during a heating experiment.
THEORY
Because of the chemical shift between the reference and temperature-dependent signal, the complete signal from one voxel will originate from two slightly different spa-tial locations. After each echo time (TE), this will have led to an additional phase difference and therefore to a different observed frequency. This means that back-ground gradients cause a difference between the true fre-quency Df and observed frefre-quency Df0(19):
Df0¼ G
G þ G0Df [1]
in which G0 is the background gradient and G is the
readout gradient. The error in frequency can then be described as the difference between the true and observed frequency: eDf ¼ Df0 Df ¼ Df G0 G þ G0 [2] with eDfbeing the error in frequency. The error will be a global
shift in frequency for a linear G0, but will vary spatially for
higher orders. eDf is dependent on the direction of the
fre-quency encoding gradient G, so it changes when the gradient direction is reversed over time during readout. Therefore, odd echoes and even echoes will result in a different Df0. As a
result, the echoes can no longer be combined into one-time series for fitting (20). This problem can be solved by using flyback gradients for rewinding and only acquiring in one readout direction. However, this lowers the spectral band-width because DTE increases. Furthermore, although eDf is
equal for all echoes, the absolute temperature is still incorrect. A possible compensation for the error lies in the use of two scans. A second scan with a readout gradient oppo-site to the first results in a complementary series, which
has an opposite gradient for each echo. By combining all echoes with a positive gradient direction in one series and all the echoes with a negative gradient direction in another, two composite series can be obtained.
The effects of the background gradient on the fre-quency estimation are the opposite in both series. The series in which the readout direction is in the same direction as the background gradient will lead to fre-quency Df0, as seen in Eq. [1]. However, in the other
series the equation becomes: Df00¼ G
G G0Df [3]
This leads to two equations with two unknowns, which makes it possible to express the true frequency Df as a combination of the two shifted frequencies:
Df ¼2Df
0
Df00
Df0þ Df00 [4]
When G0<<G, which will be the case in most situations,
the true frequency can also by approximated by taking the average of Df0and Df00:
Df0þ Df00 2 ¼ Df 1 þG0 G þ Df 1 G0 G 2 1 þG0 G 1 G0 G ¼ Df 1 G0 G 2 Df [5] METHODS
All scans were performed on a 3T system (Achieva, Phi-lips Healthcare, Best, The Netherlands) using ethylene glycol-filled phantoms. Ethylene glycol is a well suited fluid for mGE MRT because the frequency difference between its two resonances is directly related to tempera-ture (21–23). A plastic sphere (10-cm diameter) filled with ethylene glycol placed inside a 16-elements knee coil (Philips SENSE T/R knee coil, Philips Healthcare, Best, The Netherlands) was used. Shimming gradients were intentionally all set to zero and then varied between scans with 60.1 mT=m in both readout and phase-encoding direction. This was done to mimic the effect of linear background gradients and to demonstrate the related artifacts. Flyback was used to fix the readout direction with respect to the background gradients. In order to achieve a higher temporal resolution and high bandwidth, two dynamics were acquired in which the DTE was kept constant but the first TE was shifted DTE=2. This resulted in a virtual DTE that was one-half the original DTE. The most relevant parameters were: two-dimensional (2D) multigradient-echo with flyback; flip angle ¼ 30; 32 echoes; TE ¼ 1.4 þ N*2.4 ms and 2.6
þ N*2.4 ms; repetition time (TR) ¼ 80 ms; voxel size 2 2 5 mm3; matrix size 64 64; total scan time 10.7 s; and readout gradient 27.4 mT=m, bandwidth 833 Hz.
Subsequently, a heating experiment was performed using a different setup. Flexible surface receive coils (Flex L, Philips Healthcare) were wrapped around a glass cylin-der filled with water (see Fig. 1). The water was homoge-neously heated by pumping water from a heater through a glass spiral inside the cylinder and was kept stable at three different temperatures (24C, 29C, and 32C). Inside
the cylinder, a test tube (height 7 cm; diameter 2 cm) filled with ethylene glycol was placed, containing two optical
probes at different heights for temperature verification (Luxtron m3300 Biomedical Lab Kit Fluoroptic Thermom-eter [Santa Clara, CA]; 1 Hz sampling).
To test the compensation method, a similar mGE sequence to the above was used—except the flyback was switched off—enabling a lower DTE and TR. To achieve sufficient bandwidth and avoid aliasing, two dynamic scans were done with opposite readout gradients. By combining the even echoes of the first dynamic scan with the odd echoes of the second dynamic scan and vice versa, a composite series was obtained for each readout direction. The most relevant parameters were: 2D multigradient-echo; flip angle¼ 30; 32 echoes; TE ¼
1.5þN*1.25 ms; TR¼ 45 ms; voxel size 2 2 5 mm3; matrix size ¼ 64 64; total scan time ¼ 6 s; readout gra-dient 27.4 mT=m; and bandwidth 800 Hz.
All data processing was done with MATLAB (The MathWorks Inc., Natick, MA). The first echo of each acquisition was removed because it disrupted the fit. Subsequently, the complex signal of the echoes of each voxel was fit in a least-squares manner to a two Gaussian peak model (24):
S ¼ A1 expð2pif1t iw1Þ expðd1t2Þ þ A2
expð2pif2t iw2Þ expðd2t2Þ [6]
with S being the total complex signal and t denoting time. For both peaks, the amplitude A, frequency f,
Gaussian decay factor d, and phase offset f were included. The peaks of a zero-filled Fourier spectrum were used as initial values for the frequencies. The amplitude ratio of the peaks was fixed to 1:2, which is known from the molecular structure of ethylene glycol. In order to compensate for background gradients, the fre-quencies resulting from the scans without flyback were averaged according to Eq. [5].
Several similar relations for linking the frequency dif-ference between both peaks of ethylene glycol to temper-ature can be found in litertemper-ature (21–23). However, the resulting temperatures of these equations differ up to 1.3 degrees, which is much more than the intended accuracy of this measurement. The discrepancy can be related to field strength, other temperature references, or conduc-tivity (25). Therefore, for this research the median Df of a 5 5 voxel region of interest (ROI) around one probe was fit against the temperature in MATLAB with R2 ¼
0.9994, resulting in this equation: T½C ¼ 200:2 105:9 Df
gB0
[7] with B0the amplitude of the static field and g the
gyro-magnetic ratio of hydrogen nuclei in MHz=T. This result is within the range of the other relationships described in literature.To quantify background gradients caused by the susceptibility distribution of the setup, a B0 map was constructed out of the scans without
fly-back. This was achieved by taking the difference of two phase images with equal readout direction (the first and third) and dividing by 2DTE. The gradient in the read-out direction of this B0 map was determined
numeri-cally in MATLAB. Df then followed from Eq. [1], and a map of the expected temperature could be calculated using Eq. [2] (19).
RESULTS
Temperature maps of the mGE MRT with flyback, that is, only one readout direction, show the direct influence of changing background gradients (see Fig. 2). When the shimming gradients vary in phase-encoding direction, no significant effect on the temperature maps can be seen. However, when gradients in the readout direction are increased, a global offset in temperature is observed. Because the gradient strength is known, the real fre-quency can be retrieved using Eq. [1]. This calculated FIG. 1. Schematic of the setup of the heating experiment.
FIG. 2. Temperature maps with different shims: all shim gradients set to zero (a), increase in the phase direction of 0.1 mT=m (b), increase in the readout gradient of 0.1 mT=m (c), adjusted temperature map using Eq. [1] (d). The gradient in the phase direction does not affect the temperature, whereas an increase in the readout direction gives a global offset. When this is known, however, it can be compensated.
frequency results in an adjusted temperature map equal to the map without added background gradients.
The background gradients are nonlinear in the heating setup because they originate from susceptibility varia-tions (Fig. 3a and 3b). This nonlinearity leads to a vari-able temperature error across the phantom (Fig. 3c). The uncorrected temperature map shows a clear gradient in vertical direction, although the temperature was constant over the phantom according to the optical probes. When the temperature error was subtracted, a more homogene-ous temperature distribution was observed inside the phantom.
The correction based on two readout directions is shown in the second row of Figure 3. Both readout direc-tions individually show a large variation in temperature across the phantom and have an average standard devia-tion (SD) of 0.45C inside a large ROI (9 21 voxels).
However, when the temperatures of both directions are averaged, as suggested in Eq. [5], the temperature inside the ROI is almost constant (SD ¼ 0.21C).
During the experiment, the optical probes were used as a gold standard. The temperature over time, given by the probes and MRT, is shown in Figure 4. The tempera-ture given by the MRT is defined as the mean and SD of the ROI mentioned above, which covers the location of both probes.
DISCUSSION
Absolute temperature estimation by mGE MRT is highly affected by background gradients in the magnetic field. Uncorrected data showed a clear temperature variation in vertical direction, matching with background gradient amplitudes. In this work, we presented a correction method leading to a much more homogeneous tempera-ture distribution in our phantom. This matched with our heating setup, which was specifically designed to pro-vide homogenous heating. Furthermore, both optical probes showed the same temperature at different heights in the phantom, which did not match with the uncor-rected temperature distribution.
Our correction method can compensate for errors caused by background gradients in a robust manner. In principle, as others have also demonstrated (19), a cor-rection can be obtained by acquiring a B0map and
deter-mining the gradient in the readout direction by a finite difference scheme. However, this has some serious draw-backs. First of all, a finite difference derivation of a measured B0 map is challenging—especially for in vivo
data—because of increased noise and artifacts. More important, it assumes that the background gradients dur-ing the B0 map and MRT acquisition are identical. An
advantage of our method is that it can cope with varying B0fields over time, which might occur due to heating or
FIG. 3. A B0map (a), its resulting background gradient map (b), a temperature map without any correction (c), and a corrected
temper-ature map (d). By removing the tempertemper-ature error estimated from the gradient of the B0map, the homogeneity of the temperature map
is improved. The second row shows a temperature map with readout in anterior direction = up (e), the same map with readout in poste-rior direction = down (f), and the average (g). The average temperature map shows the expected constant temperature inside the phan-tom. The errors in temperature at the top and bottom of the phantom are caused by signal drops due to large local susceptibility differences.
system drift, as long as the corresponding time constants are longer than the acquisition time of two series.
In this work, we acquired two series in order to have sufficient bandwidth to sample the beating frequency at 3T. In principle, one series in which the odd and even echo would be separately fitted should suffice. The use of a lower field strength (e.g., 1.5T) would relax the require-ment with respect to the echo-spacing, making such a single-series approach potentially feasible. At higher field strengths, the errors become more serious, as can be seen from Eq. [1], because larger frequency differences and more pronounced background gradients occur. This makes the compensation more valuable. Therefore, the require-ments for the sequence also become more strict, with the need for a higher readout gradient and even shorter TEs.
After correction, the deviations between the MRT and probe temperature readings were small. However, this was the result of a frequency difference-temperature cali-bration that was performed with the same optical probes. After correction, we found that the precision of the measurements was sufficient for measuring RF heating. The SD over the large ROI inside the phantom after cor-rection (0.21C) was twice smaller than before correction
(0.45C). Furthermore, the high spatial resolution of the
scans makes it possible to detect even the smallest hot-spots (2 2 mm).
This method can also be useful in vivo; even after shimming, the residual background gradients will be pres-ent. Their amplitude depends on the target that is scanned and the direction of the readout gradient. For
example, background gradients due to susceptibility changes have been shown to be in the order of 0.1 mT=m for the human breast—and can reach up to 0.5 mT=m for the human head at 3T (26,27). Because water is the main component in tissue and fat is the most obvious choice for a reference, the resulting frequency difference becomes twice as large as in ethylene glycol. This will not only require a higher bandwidth but will also increase the error in temperature estimation. This underlines the importance of the presented correction scheme. Using fat as a reference will also result in additional difficulties: It consists of multiple peaks; the spectral position can be orientation-dependent in muscle tissue (28); and its mag-netic susceptibility is temperature-dependent (29). Fur-thermore, fat and water should be simultaneously present and well-mixed in a voxel in measureable quantities. This is because separated quantities of water and fat will lead to susceptibility errors, which are not compensated in the presented method. The limiting factors for in vivo use were not addressed in this study. However, we believe that by using the presented correction scheme, a large sys-tematic error related to background gradients can be elim-inated, significantly advancing this technique for the use of RF safety monitoring.
CONCLUSION
In RF safety assessment, absolute temperatures are an important feature. Therefore, mGE MRT can be a valua-ble tool, but these scans are known to be affected by FIG. 4. The lines in the top graph represent the temperature measured by the probes over time, the markings show the mean and standard deviation (SD) of the MR thermometry region of interest (ROI) in both directions (blue and red) and combined (green).The bot-tom graph shows the difference in temperature between the probes and the thermometry measurement. The SD of the green ROI is smallest due to the increased homogeneity.
background gradients. However, by combining two acquisitions with opposite readout direction, the temper-ature maps greatly improve. This method works during heating experiments and is stable over time. It improves absolute MRT of ethylene glycol and offers the possibil-ity to measure RF heating in phantoms more precisely. ACKNOWLEDGMENT
This research is financially supported by ZonMw (The Hague, The Netherlands) through the Electromagnetic Fields & Health program: Measuring EMF induced tissue heating and physiological changes in-vivo (85300005). REFERENCES
1. De Senneville BD, Mougenot C, Quesson B, Dragonu I, Grenier N, Moonen CTW. MR thermometry for monitoring tumor ablation. Eur Radiol 2007;17:2401–2410.
2. Sapareto SA, Dewey WC. Thermal dose determination in cancer ther-apy. Int J Radiat Oncol Biol Phys 1984;10:787–800.
3. Yarmolenko PS, Moon EJ, Landon C, et al. Thresholds for thermal damage to normal tissues: an update. Int J Hyperthermia 2011;27: 320–343.
4. Shapiro E, Borthakur A. Fast MRI of RF heating via phase difference mapping. Magn Reson Med 2002;498:492–498.
5. Cline H, Mallozzi R, Li Z, McKinnon G, Barber W. Radiofrequency power deposition utilizing thermal imaging. Magn Reson Med 2004; 51:1129–1137.
6. Oh S, Webb AG, Neuberger T, Park B, Collins CM. Experimental and numerical assessment of MRI-induced temperature change and SAR distributions in phantoms and in vivo. Magn Reson Med 2010;63: 218–223.
7. Seifert F, W€ubbeler G, Junge S, Ittermann B, Rinneberg H. Patient safety concept for multichannel transmit coils. J Magn Reson Imaging 2007;26:1315–1321.
8. Van Lier ALHMW, Kotte ANTJ, Raaymakers BW, Lagendijk JJW, van den Berg CAT. Radiofrequency heating induced by 7T head MRI: ther-mal assessment using discrete vasculature or Pennes’ bioheat equation. J Magn Reson Imaging 2012;35:795–803.
9. Murbach M, Neufeld E, Capstick M, et al. Thermal tissue damage model analyzed for different whole-body SAR and scan durations for standard MR body coils. Magn Reson Med 2013;000:1–11.
10. International Electrotechnical Commission. International standard, medical equipment–IEC 60601-2-33: particular requirements for the safety of Magnetic resonance equipment, 3rd ed. Geneva, Switzer-land: IEC; 2010.
11. Laakso I, Hirata A. Dominant factors affecting temperature rise in simulations of human thermoregulation during RF exposure. Phys Med Biol 2011;56:7449–471.
12. Foster KR, Adair ER. Modeling thermal responses in human subjects following extended exposure to radiofrequency energy. Biomed Eng Online 2004;3:4.
13. Nelson DA, Curran AR, Nyberg HA, Marttila EA, Mason PA, Ziriax JM. High-resolution simulations of the thermophysiological effects of human exposure to 100 MHz RF energy. Phys Med Biol 2013;58: 1947–1968.
14. Lagendijk J. Hyperthermia treatment planning. Phys Med Biol 2000; 45:61–76.
15. Van Leeuwen GM, Hand JW, Lagendijk JJ, Azzopardi D V, Edwards a D. Numerical modeling of temperature distributions within the neo-natal head. Pediatr Res 2000;48:351–356.
16. Flyckt VMM, Raaymakers BW, Lagendijk JJ. Modelling the impact of blood flow on the temperature distribution in the human eye and the orbit: fixed heat transfer coefficients versus the Pennes bioheat model versus discrete blood vessels. Phys Med Biol 2006;51:5007–5021. 17. Rieke V, Butts Pauly K. MR thermometry. J Magn Reson Imaging
2008;27:376–390.
18. Sprinkhuizen SM, Bakker CJG, Bartels LW. Absolute MR thermometry using time-domain analysis of multi-gradient-echo mag-nitude images. Magn Reson Med 2010;64:239–48.
19. Sprinkhuizen S. The influence of background gradients in multi gradient-echo MR thermometry. In Proceedings of the 17th Annual Meeting of ISMRM, Honolulu, Hawaii, 2009. p. 4406.
20. Du W, Du YP, Fan X, Zamora MA, Karczmar GS. Reduction of spec-tral ghost artifacts in high-resolution echo-planar spectroscopic imag-ing of water and fat resonances. Magn Reson Med 2003;49:1113–20. 21. Van Geet AL. Calibration of the methanol and glycol nuclear
mag-netic resonance thermometers with a static thermistor probe. Anal Chem 1968;40:2227–2229.
22. Ammann C, Meier P, Merbach A. A simple multinuclear NMR ther-mometer. J Magn Reson 1982;46:319–321.
23. Raiford DS, Fisk CL, Becker ED. Calibration of methanol and ethylene glycol nuclear magnetic resonance thermometers. Anal Chem 1979; 51:2050–2051.
24. Vanhamme L, Sundin T. MR spectroscopy quantitation: a review of timedomain methods. NMR Biomed 2001;14:233–246.
25. Peters RD, Henkelman RM. Proton-resonance frequency shift MR thermometry is affected by changes in the electrical conductivity of tissue. Magn Reson Med 2000;43:62–71.
26. Sprinkhuizen SM. Thesis: MR Thermometry for guidance of thermal therapy. Utrecht, The Netherlands: University of Utrecht; 2010. 27. Baudrexel S, Volz S, Preibisch C, et al. Rapid single-scan
T2*-map-ping using exponential excitation pulses and image-based correction for linear background gradients. Magn Reson Med 2009;62:263–268. 28. Boesch C, Slotboom J, Hoppeler H, Kreis R. In vivo determination of
intra-myocellular lipids in human muscle by means of localized 1H-MR-spectroscopy. Magn Reson Med 1997;37:484–493.
29. Sprinkhuizen SM, Bakker CJG, Ippel JH, Boelens R, Viergever M a, Bartels LW. Temperature dependence of the magnetic volume suscep-tibility of human breast fat tissue: an NMR study. MAGMA 2012;25: 33–39.
