Improved magnetic resonance temperature imaging through pulse sequence and reconstruction techniques

This dissertation focuses on improving magnetic resonance (MR) thermometry techniques. The goal is to develop and improve data acquisition and image reconstruction methods such that MR can be used to effectively monitor a variety of thermal therapy procedures, with high intensity focused ultrasound (HIFU) treatments being of particular interest. These treatments selectively heat a small volume of tissue rapidly, making the temperature monitoring aspect both vital and challenging. Magnetic resonance imaging can measure and map temperature changes in tissue based on the temperature dependence of the resonant frequency of water protons. The first set of results given quantifies the effects that different spatial sampling properties have on these MR temperature maps during HIFU treatments. The study indicates that accurate measurement of a single point HIFU induced temperature distribution requires a spatial resolution of 1x1x3 mm or less. Covering the entire volume of interest at this resolution would take an excessively long time to acquire using conventional MR thermometry means. The remainder of the work presented in this dissertation concerns two reconstruction techniques that allow MRI temperature maps to be acquired at an accelerated rate. The first method is known as temporally constrained reconstruction (TCR). This algorithm generates images from undersampled data by iteratively minimizing a cost function.

IMPROVED MAGNETIC RESONANCE TEMPERATURE IMAGING THROUGH PULSE SEQUENCE AND RECONSTRUCTION TECHNIQUES by Nicholas Todd A dissertation submitted to the faculty of The University of Utah in partial fulfillment of the requirements for the degree of Doctor of Philosophy Department of Physics and Astronomy The University of Utah August 2010 Copyright © Nicholas Todd 2010 All Rights Reserved The University of Utah Graduate School STATEMENT OF DISSERTATION APPROVAL The dissertation of Nicholas Todd has been approved by the following supervisory committee members: Orest Symko , Chair 03/15/2010 Date Approved Wayne Springer , Member 03/15/2010 Date Approved Werner Gellerman , Member 03/15/2010 Date Approved Eun-Kee Jeong , Member 03/09/2010 Date Approved Dennis L. Parker , Member 03/12/2010 Date Approved and by David Kieda , Chair of the Department of Physics and Astronomy and by Charles A. Wight, Dean of The Graduate School. ABSTRACT This dissertation focuses on improving magnetic resonance (MR) thermometry techniques. The goal is to develop and improve data acquisition and image reconstruction methods such that MR can be used to effectively monitor a variety of thermal therapy procedures, with high intensity focused ultrasound (HIFU) treatments being of particular interest. These treatments selectively heat a small volume of tissue rapidly, making the temperature monitoring aspect both vital and challenging. Magnetic resonance imaging can measure and map temperature changes in tissue based on the temperature dependence of the resonant frequency of water protons. The first set of results given quantifies the effects that different spatial sampling properties have on these MR temperature maps during HIFU treatments. The study indicates that accurate measurement of a single point HIFU induced temperature distribution requires a spatial resolution of 1x1x3 mm or less. Covering the entire volume of interest at this resolution would take an excessively long time to acquire using conventional MR thermometry means. The remainder of the work presented in this dissertation concerns two reconstruction techniques that allow MRI temperature maps to be acquired at an accelerated rate. The first method is known as temporally constrained reconstruction (TCR). This algorithm generates images from undersampled data by iteratively minimizing a cost function. The cost function has terms that penalize deviations from the acquired data and sharp changes in time. Results demonstrate that three- to fourfold data acquisition accelerations can be achieved while keeping a temperature root mean square error (RMSE) below 1.0 °C. The second method, Model Predictive Filtering (MPF), uses a more sophisticated model based on the Pennes bioheat equation. Predictions from this model are combined with the undersampled data to create temperature maps in near real time. Results for the MPF technique demonstrate that data acquisition acceleration rates of 5 to 12 can be achieved while keeping a temperature RMSE below 1.0 °C. Each technique can be used with a number of different MR thermometry techniques and in combination with any pulse sequence. The acceleration factors achieved by these algorithms represent a large step towards realizing accurate guidance and control of thermal therapies. iv CONTENTS ABSTRACT ………………….…………………………………………………………iii LIST OF FIGURES .………….………………………………………………………...viii LIST OF TABLES………………….…………………………………………………...xiv LIST OF ACRONYMS…………….…………………………………………………....xv ACKNOWLEDGEMENTS …………………………………………………………....xvii CHAPTERS 1. PRINCIPLES OF MAGNETIC RESONANCE IMAGING................................ 1 1.1. Introduction................................................................................................... 1 1.2. Nuclear Magnetic Resonance Background................................................... 3 1.3. Signal Generation - Quantum View............................................................. 5 1.3.1. Spin in a Magnetic Field.................................................................... 5 1.3.2. The Zeeman Effect and Net Magnetization....................................... 8 1.4. Signal Generation - Classical View............................................................. 10 1.5. Signal Detection............................................................................................ 16 1.6. Signal Localization....................................................................................... 22 1.6.1. Frequency Encoding.......................................................................... 22 1.6.2. Phase Encoding.................................................................................. 23 1.6.3. 2-D and 3-D Imaging......................................................................... 24 1.6.4. A k-space Interpretation..................................................................... 26 1.7. Image Reconstruction................................................................................... 28 2. MR GUIDED HIGH INTENSITY FOCUSED ULTRASOUND 30 2.1. Ultrasound Physics........................................................................................ 31 2.2. HIFU Challenges........................................................................................... 34 2.3. Need for Temperature Monitoring................................................................. 36 3. MRI THERMOMETRY....................................................................................... 38 3.1. Temperature Dependence of the T1 Relaxation Time of Water Protons....... 38 3.2. Temperature Dependence of the Diffusion Coefficient of Water................. 44 3.3. Temperature Dependence of the Proton Resonance Frequency.................... 46 3.4. Challenges with the PRF Method.................................................................. 52 3.5. Overcoming MR Temperature Imaging Limitations..................................... 55 4. SPATIAL SAMPLING CONSIDERATIONS FOR ACCURATE MR TEMPERATURE MEASUREMENTS................................................................ 63 4.1. Introduction.................................................................................................... 63 4.2. Methods.......................................................................................................... 65 4.2.1. Theory................................................................................................. 65 4.2.2. Simulations.......................................................................................... 70 4.2.3. Experiments........................................................................................ 73 4.2.4. Data Analysis...................................................................................... 74 4.3. Results........................................................................................................... 75 4.3.1. Simulations.......................................................................................... 75 4.3.2. Experiments........................................................................................ 79 4.4. Discussion...................................................................................................... 82 5. TEMPORALLY CONSTRAINED RECONSTRUCTION APPLIED TO MRI TEMPERATURE DATA..................................................................................... 88 5.1. Introduction.................................................................................................... 88 5.2. Methods............................................................................….......................... 91 5.2.1. TCR Theory........................................................................................ 91 5.2.2. TCR Applied to PRF Temperature Data............................................. 93 5.2.3. MRI Acquisition and Heating Set Up................................................. 95 5.2.4. Data Analysis...................................................................................... 97 5.3. Results............................................................................................................ 99 5.3.1. Optimization........................................................................................ 99 5.3.2. Maximum Reduction Factor….......................................................... 102 5.3.3. Comparison to Sliding Window and Low Resolution Reconstructions...................................................................................... 111 5.4. Discussion..................................................................................................... 113 6. MODEL PREDICTIVE FILTERING FOR IMPROVED TEMPORAL RESOLUTION IN MRI TEMPERATURE IMAGING.................................... 120 6.1. Introduction................................................................................................. 120 6.2. Methods....................................................................................................... 124 6.2.1. The Thermal Model......................................................................... 124 6.2.1.1. Parameter Determination........................................................ 125 6.2.1.2. Model Implementation........................................................... 126 6.2.2. The Model Predictive Filtering Algorithm...................................... 126 vi 6.2.3. MRI Acquisition and HIFU Heating Set Up................................... 128 6.2.4. Experiments..................................................................................... 130 6.2.4.1. MPF with Retrospectively Undersampled Data Sets.............. 130 6.2.4.2. MPF with Actually Undersampled Data Sets.......................... 131 6.2.4.3. MPF with Motion...................................................................... 132 6.2.4.4. Data Analysis.......................................................................... 133 6.3. Results......................................................................................................... 134 6.3.1. Model Performance.......................................................................... 134 6.3.2. Heating Repeatability with Fully Sampled Data Sets….................. 135 6.3.3. MPF with Retrospectively Undersampled Data Sets....................... 135 6.3.4. MPF with Actually Undersampled Data Sets.................................. 138 6.3.5. MPF with Motion............................................................................ 141 6.3.6. MPF Robustness.............................................................................. 143 6.4. Discussion................................................................................................... 145 7. CONCLUSIONS................................................................................................ 151 REFERENCES................................................................................................... 155 vii LIST OF FIGURES Figure 1.1 The Zeeman effect. A spin ½ particle placed in an external field will have two energy levels............................................................................................................. 8 1.2 Relaxation behavior for a spin system with T2 = 50 ms and T1 = 600 ms that has been tipped into the transverse plane with a 90° RF pulse. The transverse component decays to zero while the longitudinal component regrows to its original value............................................................................................................ 16 1.3 FID signals. A) FID signal from a single spin isochromat. B) FID signal from a system with two different spin isochromats.............................................................. 18 1.4 Schematic showing the behavior of the magnetization during the spin echo signal acquisition process.................................................................................................... 19 1.5 Timing diagram for the spin echo acquisition, showing the signal decaying at a rate determined by T2 * and being refocused to an amplitude determined by T2.......................................................................................................................... 20 1.6 Timing diagram for the gradient echo showing the signal decaying at a rate faster than T2 * and being refocused to an amplitude determined by T2 *............. 21 1.7 Pulse sequence timing diagram. A) Timing diagram for gradient echo sequence showing the RF pulses, all gradient waveforms, and ADC signal acquisition. B) Cartesian k-space path determined by x- and y gradients........... 27 1.8 MR k-space and image. A) Fully sampled k-space data set from a gradient echo sequence for a 128 x 128 image. B) Magnitude image created from k-space data............................................................................................................. 29 2.1 An ultrasound wave, VI, incident on an interface where the two different media has different impedances, Z, speeds of sound, c, and absorptions, μ. Some of the ultrasound wave will be reflected, and the transmitted portion will be refracted................................................................................................... 32 2.2 Schematic of the set up for a HIFU tumor ablation. Ultrasound energy is used to selectively heat the tumor to the point of necrosis. A number of potential challenges exist, including 1) rapid temperature rises at the focus, 2) near field heating, 3) beam distortion due to fat, 4) heating at tissue/bone interfaces, and 5) motion of the target organ............................................................................... 35 3.1 Local field effects T1. A) Randomly varying local magnetic field at the site of a water proton. B) Frequency spectrum density for fields varying at different rates....................................................................................................... 40 3.2 Plot of T1 values against temperature. The temperature dependence of T1 shows good linearity over a small range of temperatures around body temperature.......................................................................................................... 41 3.3 Apparent diffusion coefficient plotted against temperature. The diffusion coefficient also shows reasonable linearity in the range of temperatures around body temperature..................................................................................... 45 3.4 A schematic depicting the temperature dependence of the water proton resonance. As temperature increases, the hydrogen bonds between water molecules change such that the hydrogen proton is more strongly screened, leading to a lower resonant frequency................................................................. 47 3.5 Phase maps from MR images acquired before heating (reference) and during heating (current time frame). The difference between the phases is proportional to the temperature change............................................................... 50 3.6 Calibration experiment results relating the phase change to the temperature change. The phase change is linear with a slope of -0.015 ppm/°C................... 51 3.7 Motion related temperature errors. A) Magnitude image of the breast, imaged during free breathing without heating. B) Temperature plot of one voxel over time showing errors due to susceptibility changes caused by respiratory motion.................................................................................................................. 54 3.8 Phase and T1 information acquired simultaneously during a phantom cooling experiment using the DESPOT1 sequence. Both metrics show good linearity.. 58 3.9 Steps of the reference-less phase difference method. Instead of subtracting the phase from an earlier acquired image, a polynomial is fit to the background phase from a non-heated region of the current image. This fitted polynomial is used as the reference phase for subtraction......................................................... 59 ix 4.1 Thermal dose rate as a function of temperature change. The thermal dose rate doubles with every 1°C increase in temperature and at 20°C above body temperature, a lethal dose is delivered every second.......................................... 66 4.2 A one dimensional plot of the voxel sensitivity function (VSF)......................... 68 4.3 Simulated temperature distribution with VSF..................................................... 69 4.4 Simulated temperature map with 0.lmm isotropic spatial resolution. Maximum temperature time frame shown in a transverse orientation (parallel to the beam path). The maximum temperature rise was 23.3°C. The simulated power deposition matrix had FWHM dimensions of 7.8 x 2.3 x 2.3mm and the corresponding simulated temperature maps had FWHM dimensions of 10.7 x 3.3 x 3.3mm at the maximum temperature time frame............................ 72 4.5 The maximum measured temperature as a function of voxel dimensions for the case of 3-D imaging with isotropic voxels (A) and the case of 2-D imaging with varying slice thickness (B). For the 2-D case, 4 scenarios were considered: 1.0 x 1.0mm or 2.0 x 2.0mm in-plane resolution with the slices oriented perpendicular or parallel to the ultrasound beam.................................. 76 4.6 Grid placements effects A) Plot of maximum measured temperature as a function of grid location. Simulated temperatures reconstructed with 2.0 x 2.0 x 2.0 mm resolution and shifted at 0.06mm steps. B - D) Three representative transverse temperature maps with maximum temperatures of 20.3°, 18.9°, 17.1°, VD30 of 32mm3, 32mm3, 32mm3, and VD240 of 16mm3, 8mm3, 0 mm3, respectively................................................................................................ 77 4.7 Spatial resolution effects, experimental results. A) Coronal temperature maps from all 3 HIFU heating runs showing the maximum measured temperature. Each row is a different heating run acquired at a set spatial resolution and each column is a reconstruction of the original data using a different ZFI factor. B) Temperature-time plots for the hottest voxel of six different reconstructions..................................................................................................... 80 4.8 Grid placement effects. A) Plot of maximum measured temperature as a function of grid location. 2-D experimental temperatures acquired at 2.0 x 2.0 x 3.0 mm resolution and retrospectively shifted at 0.06mm steps in the x-direction only. B - D) Three representative coronal temperature maps with maximum temperatures of 26.3°, 24.2°, 21.8°, VD30 of 40mm3, 64mm3, 64mm3, and VD240 of 32mm3,40mm3, 32mm3, respectively............................. 81 x 5.1 Determining optimal values for α and number of iterations. A) The temperature RMSE of a single time frame from the first the ex vivo porcine training data set as a function of number of iterations. Convergence is almost complete after 100 iterations. B) The temperature RMSE for a number of different α/iteration combinations from the same data set................................... 101 5.2 Magnitude images and the corresponding temperature maps of time frame 16 of one ex vivo porcine tissue validation data set are shown. The full data are shown on top (A and D), TCR at R=4 is shown in the middle (B and E) and the difference images are shown on the bottom (C and F). The temperature difference map is scaled from -2°C to 2°C........................................................... 103 5.3 Correlation between error and temperature changes. A) Absolute value of temperature error plotted against the absolute value of dT/dt for all voxels in the heated region. B) Absolute value of temperature error plotted against temperature for all voxels in the heated region................................................... 104 5.4 Temperature versus time curves from the hottest voxels of four different ex vivo porcine tissue data sets. The TCR temperatures at R = 4 agree quite well with the full data temperature.............................................................................. 104 5.5 Correlation between error and temperature changes. A) Absolute value of temperature error plotted against the absolute value of dT/dt for all voxels in the heated region. B) Absolute value of temperature error plotted against temperature for all voxels in the heated region................................................... 107 5.6 Temperatures maps from in vivo data. A) Temperature map from in vivo canine data reconstructed with the full data. B) Same temperature map reconstructed with the TCR algorithm at R=3. C) Temperature difference map scaled from -2°C to 2°C........................................................................................ 109 5.7 Temperature plots from in vivo data. A) through C): Temperature versus time curves for a voxel in the focal zone from the three in vivo canine data sets. The TCR temperatures at R=3 agree quite well with the full data temperatures. D) through F): The mean and standard deviation of temperature error for an ROI over the heated region from the same three data sets. D) is the training data set, E) and F) are the validation data sets. The horizontal bars represent times when the heating is on. The nosier canine data had a larger standard deviation than the ex vivo porcine tissue data sets, even at R=3. Time frames acquired every 6.2 seconds................................................................................... 110 xi 5.8 Mean and standard deviation of temperature error over an ROI for an ex vivo porcine tissue validation data set. The TCR algorithm at R=4 is shown in A); the sliding window reconstruction in shown in B); and the low resolution reconstruction is shown in C). Time frames acquired every 6.2 seconds........... 112 5.9 Dose versus time plots, comparing the TCR algorithm to the full data. The plots in A) and B) are from ex vivo porcine tissue data sets and the plot in C) is from an in vivo canine data set. The dashed black lines represent the dose that would be calculated from the full data if the temperatures were systematically shifted by +/- 1°C. Time frames acquired every 6.2 seconds................................................................................................................ 117 6.1 Schematic showing the experimental set up. An MRI compatible 256-element phased array ultrasound transducer is housed in a water bath, coupled to ex vivo porcine tissue samples. In-house built surface coils are used for imaging.. 130 6.2 Model performance. A) Temperature plot of one voxel in the focal zone comparing fully sampled PRF temperatures with model only temperatures for 48W heating. B) Mean and STD of the model only temperature errors over a 24x7x3 ROI for all time frames........................................................................... 134 6.3 Comparing the MPF temperatures of run 1.3 versus the fully sampled PRF temperatures of run 1.3 (fully sampled PRF temperatures minus MPF temperatures; mean and STD over a 24x7x3 ROI). The temperature errors for the MPF algorithm at three different reduction factors are shown................................................................................................................... 136 6.4 The error metrics as a function of reduction factor. A) The error over the 5 hottest voxels in the hottest time frame (mean and STD over all 10 heating runs) at the 6 different reduction factors. B) Temperature RSME over a 24x7x3x25 ROI and all 10 heating runs at 6 different reduction factors............ 137 6.5 2-D MPF temperatures compared to 2-D fully sampled temperatures. A) Mean and STD of the hottest voxel in the focal zone over the first 4 runs of fully sampled PRF temperatures (experiment 2.A, black line) and 2-D MPF temperatures (experiment 2.B, white line). B) Mean and STD of the temperature differences between runs 2.A and 2.B over a 24x7x3 ROI...................................................................................................................... 139 6.6 3-D MPF temperatures compared to 2-D fully sampled temperatures. A) Mean and STD of the hottest voxel over the first four runs of fully sampled PRF temperatures (black line) and 3-D MPF temperatures (white line). B) Mean and STD of temperature differences over a 24x7x3 ROI for all four runs...................................................................................................................... 140 xii 6.7 2-D undersampled MPF temperatures with discrete motion. A) Temperature plots of three voxels in and near the focal zone, spaced 6mm apart. B) Phantom displacement relative to the starting point, as measured by the method of Mendes et al. C) Temperature maps showing the hot spot before moving, just after the first move, at peak heating, and just after the second move.................................................................................................................... 142 6.8 2-D undersampled MPF temperatures with continuous motion. A) Temperature plots of two voxels in the focal zone, spaced 12mm apart. B) Phantom displacement relative to the starting point, as measured by the method of Mendes et al. C) Temperature maps showing the hot spot during heating at four consecutive times over one half of a period................................ 143 6.9 Evaluating the effect of model accuracy on the MPF algorithm. MPF temperatures were created using erroneous values for Qrel, k, and ρ and were compared to fully sampled temperatures. A) Average difference over the 5 hottest voxels as a function of model error. B) Temperature RMSE over a 24x7x3x25 ROI as a function of model error...................................................... 144 xiii LIST OF TABLES Table 4.1 Thermal and acoustic parameters used for numerical simulation of temperature maps................................................................................................ 71 4.2 Effects of sampling grid location. Simulated temperature maps with 0.1 mm3 isotropic voxels had a maximum temperature rise of 23.3°C, VD30 = 30.2 mm3 and VD240 = 10.8 mm3. For the same data reconstructed at 1.0mm3 and 2.0mm3 isotropic resolution, the measured maximum temperatures and dosed volumes were heavily dependent on the sampling grid location........................ 78 4.3 Effects of sampling grid location on experimentally acquired data at different spatial resolutions. Grid shifting was done by applying a linear phase in k-space. The range of measured values are reported for each metric..................... 82 5.1 α values, number of iterations, and temperature RMSE for all parameter sets. The training data was used to determine the α values and number of iterations. Only the validation data sets were used to determine the temperature RMSE.............................................................................................. 102 5.2 Temperature RMSE values for TCR algorithm, sliding window reconstruction and low resolution reconstruction....................................................................... 112 6.1 Naming and description of all 39 HIFU heating experiments............................ 131 6.2 Summary of results for HIFU heating experiments. Temperature differences over the 5 hottest voxels in the hottest time frame and a temperature RMSE over a 24x7x3x25 ROI are displayed.................................................................. 140 6.3 Summary of results comparing MPF and full data temperature maps as a function of image SNR using k-space data from run 1.2.................................... 145 LIST OF ACRONYMS ADC Analog to Digital Converter CEM Cumulative Equivalent Minutes DESPOT1 Driven Equilibrium Single Pulse Observation of T1 EPI Echo Planar Imaging FFT Fast Fourier Transform FID Free Induction Decay GRE Gradient Echo FWHM Full Width at Half Maximum HIFU High Intensity Focused Ultrasound MPF Model Predictive Filtering MRgHIFU Magnetic Resonance guided High Intensity Focused Ultrasound MRI Magnetic Resonance Imaging NAA N-acetyl-aspartate NMR Nuclear Magnetic Resonance ppm parts per million PRF Proton Resonance Frequency RF Radio Frequency RIGR Reduced-encoding MR Imaging with Generalized-series Reconstruction RMSE Root Mean Square Error ROI Region of Interest SENSE Sensitivity Encoding SMASH Simultaneous Acquisition of Spatial Harmonics SNR Signal to Noise Ratio TCR Temporally Constrained Reconstruction TE Echo Time TR Repetition Time UNFOLD Unaliasing by Fourier-encoding the Overlaps Using the Temporal Dimension VSF Voxel Sensitivity Function ZFI Zero Filled Interpolation xvi ACKNOWLEDGEMENTS There have been numerous people who have helped me immeasurably throughout my time as a PhD student. I would not have been able to complete my degree without their assistance and support. I would like to express my gratitude to everyone and acknowledge a few individuals who have been especially important to my education and research during these years. Special thanks go to Professor Dennis Parker for being an outstanding advisor and mentor. Recognition is also due to my colleagues on the HIFU research project, in particular Allison Payne, Josh De Bever, and Urvi Vyas. I am grateful to Professors Ed DiBella and EK Jeong for their collaboration and help over the years. I want to thank the rest of my colleagues at UCAIR, my past professors, and the members of my committee. Finally, I am very appreciative of help and support given to me by my family and friends during my time as a graduate student. CHAPTER 1 PRINCIPLES OF MAGNETIC RESONANCE IMAGING 1.1. Introduction The purpose of this thesis is to present work done on improving magnetic resonance temperature imaging. Thermal therapies are increasingly being used for applications such as drug delivery or tumor ablation. Such therapies seek to selectively heat a targeted area of tissue until a temperature related end goal is reached. In order to carry out these procedures safely and effectively, the process must be continuously monitored. Magnetic resonance imaging (MRI) is ideally suited for the task due to its ability to detect irreversible tissue damage and measure temperature changes in a variety of tissue types. Proper monitoring requires that the MR temperature maps have spatial resolution adequate to accurately measure the induced temperature distribution, temporal resolution fast enough to track changes in regions where the temperature is increasing most rapidly, and sufficient volume coverage to image everywhere where thermal energy may be deposited. Currently, these needs are not met for all applications. The goal of this work is to develop MRI data acquisition and image reconstruction techniques to improve the spatial resolution, temporal resolution, and volume coverage of MR temperature imaging methods. 2 This thesis is divided into seven chapters. In the rest of Chapter 1, background and theory on MRI are presented. Both quantum mechanical and classical treatments are given for the magnetic resonance phenomenon and basic imaging techniques are presented. In Chapter 2, a brief introduction to MR-guided high intensity focused ultrasound (MRgHIFU) is given. HIFU is a thermal procedure that can be used for tumor ablation and is the tissue heating technique used in most of the experiments described in this thesis. Chapter 3 presents a thorough discussion of MR thermometry methods, including the physical basis of temperature-dependent MR parameters, current state-of-the- art techniques, and remaining challenges. In Chapter 4, original results are given quantifying the effects of MR spatial sampling properties on temperature measurements during HIFU heating. The minimum spatial resolution necessary for accurately measuring the HIFU induced temperature increases is determined. Achieving this spatial resolution using conventional MR thermometry would require excessively long scan times or minimal volume coverage. Chapters 5 and 6 present two novel MR temperature imaging methods for overcoming this challenge. Each undersamples the MR data to save acquisition time and then employs a novel image reconstruction technique to create artifact-free images from the undersampled data. The reconstruction technique presented in Chapter 5 is called Temporally Constrained Reconstruction. In this method, the artifact free image is created from the undersampled data by iteratively minimizing a cost function that contains two terms. The first term is a data fidelity term that penalizes image estimates that deviate from the originally acquired data and the second term is a temporal constraint term that penalizes image estimates that change abruptly in time. The implicit model used in this approach only assumes that whatever changes are taking place 3 happen slowly in time. In Chapter 6, a second novel reconstruction algorithm is presented where a more explicit model is used as part of the reconstruction process in a technique called model predictive filtering. This method combines temperature evolution predictions made from the Pennes bioheat equation with the undersampled data to create the temperature maps. Chapter 7 summarizes the work accomplished in this thesis, draws conclusions about the importance of the imaging methods presented, and describes potential future work to be done. 1.2. Nuclear Magnetic Resonance Background Magnetic resonance imaging (MRI) has been an incredibly successful imaging modality due to its ability to noninvasively acquire high resolution images of the internal structure of subjects without the use of ionizing radiation. The MRI signal is dependent on a number of different physical parameters that can be utilized to provide various contrasts amongst different soft tissue types. This flexibility allows MRI to be used for a range of purposes, including discrimination between normal and pathological tissues for diagnosis and tracking dynamic changes in tissue properties over time. Unlike other imaging modalities, MRI generates its signal source from within the image subject itself. In order to form an image, this signal must be created, detected, and localized. The creation of the signal is based on the phenomenon of nuclear magnetic resonance (NMR). While MRI is a fairly new field, being established in the 1970s and early 1980s, work on NMR principles dates back to the late 1930s. The NMR phenomenon occurs for certain types of nuclei when they are inserted into a magnetic field and exposed to electromagnetic radiation of the proper frequency. The nuclei 4 absorb the energy from the electromagnetic pulse, and a time-varying magnetic field is created whose frequency depends on the proton density of the sample, the strength of the magnetic field, the chemical environment of the nuclei and other factors. A signal can be detected by tuning a resonant coil to the frequency of this time-varying magnetic field. Isidor Rabi first demonstrated the NMR phenomenon in 1939 using molecular beams, an achievement for which he was awarded the Nobel Prize in Physics in 1944 [1]. A few years later, in 1946, independent groups working under Felix Bloch and Edward Purcell separately published methods for detecting the NMR phenomenon in solid paraffin and water, respectively [2-3]. The two groups' important contribution of extending the application to solids and liquids was recognized with a shared Nobel Prize in Physics in 1952. For the next 20 years, NMR research focused on determining the molecular and chemical properties of bulk materials. Then in the mid-1970s a number of investigators developed ideas for localizing the NMR signal, leading to the creation of two dimensional images. In 1973 Paul Lauterbur published an article in Nature demonstrating that a weaker spatially varying magnetic field added to main magnetic field would impart spatially dependent frequency differences to the nuclei that generate the total signal [4]. With knowledge of the gradient field, a given resonant frequency could be mapped back to a certain location from within the sample, and the strength of the signal from that location would indicate the relative concentration of the nuclei from that region. In the same year, Peter Mansfield showed that an MR signal acquired in the presence of a gradient field was related to the nuclei density through the Fourier transform [5]. Two years later, Kumar, Welti, and Ernst extended this idea to a 2-D 5 Fourier imaging method, which is still used today [6]. The transformation from NMR to MRI had begun and Lauterbur, Mansfield, and Ernst all received Nobel Prizes for their contributions. In the following sections, a more thorough treatment of MR signal generation, detection, and localization is presented. The theory is outlined based on the excellent works by Griffiths, Slichter, and Liang and Lauterbur [7-9]. 1.3. Signal Generation - Quantum View All elementary particles possess an intrinsic form of angular momentum known as spin. Unlike the orbital angular momentum of a particle, which can take on different discrete values, the spin of a particle is specific and fixed. Protons, neutrons, and electrons each have spin ½, and atomic nuclei that are made up of protons and neutrons can have 0, ½-integer, or integer spin quantum number. The magnetic resonance phenomenon occurs for all nuclei that have nonzero spin quantum number. MR imaging and spectroscopy are performed with a number of different spin ½ nuclei, including 1H, 3He, 13C, 19F, and 31P. Other nuclei such as 17O (spin 5/2) and 23Na (spin 3/2) are also used. MRI is most commonly performed using the spin ½ hydrogen nucleus due to its large gyromagnetic ratio and natural abundance in the human body. The work in this thesis is based solely on 1H MRI studies and therefore the theory presented will only deal with spin ½ systems. 1.3.1. Spin in a Magnetic Field Consider an isolated particle with spin s and the quantum operator for spin angular momentum: 6 S S i S j S k x y z = ˆ + ˆ + ˆ r [1.1] As with the orbital angular momentum operator, eigenvalues for S2 and only one of Sx, Sy, or Sz can be simultaneously specified [7]. The eigenvalues for S2 will be ћs(s+1) and the eigenvalues for Sz will be ћm, where s = 0, ½, 1, 3/2, …, m = -s, -s+1, … s, and ћ is Planck's constant divided by 2π. For a spin ½ nucleus, s = ½ and therefore there are only two possible eigenstates, spin up ( ↑ ) with m = ½ and spin down ( ↓ ) with m = -½, and the general state can be specified as a linear combination of the two eigenstates Ψ = ↑ + ↓ + − c c [1.2] where c+ and c- are complex constants determined by initial conditions. A nucleus will possess spin angular momentum, S r , and a magnetic moment, μ r , and relation between the two is given by S r r μ = γ [1.3] where γ is called the gyromagnetic ratio and is different for different nuclei. When such a nucleus is placed in a uniform magnetic field, B r , the interaction energy between the field and the magnetic moment is E B r r = −μ ⋅ [1.4] For a magnetic field of strength B0 aligned along the z-direction, this leads to a Hamiltonian of z B S 0 Η = −γ [1.5] 7 The eigenvalues of this time-independent Hamiltonian will be the same as those for the operator Sz, scaled by γB0. The energies of the 2s + 1 different eigenstates will therefore be: E B m m = −γh 0 [1.6] Summing the individual eigenstates and adding the time dependence will give the general solution: ( ) ( ) Σ− = Ψ = − s m s m mt c s,m exp iE t h [1.7] For a spin ½ system with only two states, Equation [1.7] reduces to: ( ) exp( 2) exp( 2) 0 0Ψ t = c ↑ − iγB t + c ↓ iγB t + − [1.8] To get a sense of the behavior of the nucleus in the magnetic field, the expectation value for the three components of the nucleus' magnetic moment can be computed according to: ( ) ( ) ( ) μ t μ (t)dτ x y z = ∫ Ψ x y z Ψ , * , [1.9] Using Equation [1.8] for the calculation, it can be shown that [9]: ( ) ( t) x 0 0 sin cos 2 θ φ ω γ μ = + h [1.10] ( ) ( t) y 0 0 sin sin 2 θ φ ω γ μ = + h [1.11] (θ ) γ μ cos 2 h = z [1.12] where ω0 = γB0 is called the Larmor frequency. The solution indicates that the magnetic moment precesses about the main field at a fixed angle, θ, with frequency ω0. 8 1.3.2. The Zeeman Effect and Net Magnetization Equation [1.6] shows that a nucleus in a magnetic field will have multiple energy levels separated by 0 0 ΔE = γhB = hω [1.13] This splitting of energy levels is known as the Zeeman effect. For spin ½ systems, in which there are only two energy levels, the spin up state corresponds with the magnetic moment being aligned parallel to the main field and has the lower energy. The energy splitting for a spin ½ system is illustrated in Figure 1.1. Transitions can occur between adjacent energy levels through stimulated or spontaneous emission. When an unmagnetized sample is placed in an external field, the sample becomes magnetized because slightly more spins enter the lower energy state. For an ensemble of spin ½ particles, the two states will be populated according to the energy difference established by the main field and the Boltzmann relationship [8]: Figure 1.1. The Zeeman effect. A spin ½ particle placed in an external field will have two energy levels. 9 ⎟⎠ ⎞ ⎜⎝ ⎛ Δ = ↓ ↑ kT E N N exp [1.14] where N↑ is the number of spins aligned parallel to the field, N↓ is the number of spins aligned anti-parallel to the field, k = 1.38 x 10-23 J/K is the Boltzmann constant, and T is the absolute temperature of the spin system. Equation [1.14] indicates that there is a very slight preference for the spins to be in the spin up, lower energy state. This can be seen by expanding the exponential term using a first order approximation kT B N N 1 0 γh ≈ + ↓ ↑ [1.15] Then using Ntot = N↑ + N↓, the population difference is obtained as: kT B N N Ntot 2 γh 0 − ≈ ↑ ↓ [1.16] Substituting in γ = 42.58 MHz/T (the gyromagnetic ratio for 1H), B0 = 3 T (the main field strength for a typical MRI scanner), and T = 310 K (human body temperature), the population difference is on the order of one in a million. Although the excess number of spins aligned with the magnetic field is very small, it still leads to an observable net magnetization due to the very large number of spins present in a sample. With each of the N↑ spins having magnetic moment μ = ½γћ, and each of the N↓ spins having magnetic moment μ = -½γћ, the net magnetization will be given by M (N N )kˆ 2 1 ↑ ↓ = h − r γ [1.17] Upon substitution of Equation [1.16], an approximation for the net magnetization is obtained as 10 kT M B Ntot z 4 0 γ 2h 2 ≈ [1.18] It can be seen that the net magnetization is oriented along the direction of the main field, is directly proportional to the main field strength and total number of spins, and inversely proportional to the temperature of the system. The establishment of a macroscopic magnetization allows the treatment to proceed from a classical point of view. 1.4. Signal Generation - Classical View In section 1.3, it was shown that quantum mechanics predicts that when a nucleus with a magnetic moment is placed in an external magnetic field the expectation values of the magnetic moment components perpendicular to the field will precess about that field. A similar results can be shown approaching the problem from the point of view of classical mechanics. According to classical mechanics, an object with a magnetic moment, μ r , placed in an external magnetic field, B r , will experience a torque that is equal to the rate of change of the object's angular momentum, J r B dt dJ r r r = μ × [1.19] Using J r r μ = γ and assuming that the magnetic field is aligned along the z-axis, equation [1.19] is converted to B k dt d ˆ 0 = γμ × μ r r [1.20] 11 Solving this differential equation for the individual components of μ r and substituting in ω0 = γB0 gives (t) ( ) ( t) ( ) ( t) x x 0 y 0 μ = μ 0 cos ω + μ 0 sin ω [1.21] (t) ( ) ( t) ( ) ( t) y x 0 y 0 μ = −μ 0 sin ω + μ 0 cos ω [1.22] ( ) (0) z z μ t = μ [1.23] where μx(0), μy(0), and μz(0) are the initial conditions. As in the quantum mechanical treatment, the classical approach predicts that the transverse components of the magnetic moment will precess about the main magnetic field at the Larmor frequency, with the longitudinal component determined by the initial conditions. In a system such as a tissue sample, the large number of magnetic dipoles can be summed together to obtain the net magnetization: M i j k n z n n y n n x n ˆ ˆ ˆ , , , ⎟⎠ ⎞ ⎜⎝ ⎛ + ⎟⎠ ⎞ ⎜⎝ ⎛ + ⎟⎠ ⎞ ⎜⎝ ⎛ = Σμ Σμ Σμ r [1.24] The x- and y- components of the individual magnetic dipoles will not be in phase as they rotate and will average out to zero when summed together. However, the z-component of the magnetic dipoles will add together to form a net magnetization, 0 z M , that is parallel to the main field with strength of the magnetization being given by Equation [1.18]. With the initial magnetization aligned along the z-direction and therefore no initial transverse components, Equations [1.21] - [1.23] indicate that the magnetization vector will remain stationary in time. In order to generate a signal, the system must be perturbed out of equilibrium to create a component of transverse magnetization. This is 12 accomplished by supplying external energy in the form of an electromagnetic wave with its accompanying oscillating magnetic field, referred to as the 1 B field. B (t) 2Be (t)cos( t )iˆ 1 1 0 = ω +φ r [1.25] where Be (t) 1 is an envelope function that can take different forms depending on the application. In order to have the desired effect on the magnetization, the 1 B field must be perpendicular to the main field and rotating at the same frequency at which the magnetization is precessing about the main field. Therefore 0 ω in Equation [1.25] must be the Larmor frequency. Note that this is consistent with the quantum point of view, where radiation with energy equal to the energy difference between the two spin states is needed to induce a transition. The B1 pulse is often referred to as an RF pulse because the electromagnetic radiation supplied is in the radio frequency range for most MRI applications. For example, the Larmor frequency for 1H imaging at 3 Telsa is 127.7 MHz. In order to describe how a magnetization vector, M r , in an external magnetic field responds to the application of a B1 field, an extension of Equation [1.19] was developed, known as the Bloch equation. In its most general form, the Bloch equation is: ( ) 1 0 2 ˆ ˆ ˆ T M M k T M i M j M B dt dM x y z − z − + = × − r r r γ [1.26] where T2 is the spin-spin relaxation time constant and T1 is the spin-lattice relaxation time constant. T2 is the time that it takes for the magnetization in the x-y plane to dephase and reach 37% of its original value, while T1 is the time that it takes for the z-magnetization to regrow to 63% of the thermal equilibrium value. 13 To simplify the treatment of B1 fields, it is convenient to transform the Bloch equation into a rotating frame of reference. This is done by defining unit vectors i (wt)i (wt) j r ˆ = cos ˆ − sin ˆ j (wt)i (wt) j r ˆ = sin ˆ + cos ˆ [1.27] k k r ˆ = ˆ and expressing the relevant vectors in terms of these new rotating unit vectors. The Bloch equation in the rotating frame then becomes ( ) 1 0 2 ˆ ˆ ˆ T M M k T M i M j M B dt dM xr r yr r zr z r r eff r − − + = × − r r r γ [1.28] where r xr r yr r zr r M = M iˆ + M ˆj + M kˆ r [1.29] r xr r yr r zr r B = B iˆ + B ˆj + B kˆ r [1.30] γ ω r r r = + eff r B B [1.31] During the application of a B1 pulse, eff B r will be made up of the main B0 field in the z-direction and the B1 component in the transverse direction. ( ) γ ω r r = + + eff r r B B t iˆ B kˆ 1 0 [1.32] If the frequency of the rotating frame is matched to the Larmor frequency, kˆ 0 ω ω = − r , then it can be seen that the z-component of eff B r is cancelled out, leaving only the time dependent B1 field to act on a magnetization that is stationary with respect to the rotating frame. The application of a B1 pulse is usually on the order of a few milliseconds, while 14 the T2 and T1 time constants are on the order of tens to hundreds of milliseconds. The second and third terms from the right hand side of Equation [1.28] can therefore be dropped when considering the effects of a B1 pulse on the magnetization. With this assumption and the substitution of Equation [1.32] with kˆ 0 ω ω = − r , the Bloch equation during an RF pulse equation becomes ( ) r r r M B t i dt dM ˆ 1 = × r r γ [1.33] A specific solution to Equation [1.33] is only possible once the form of B1(t) has been specified. An exact solution can be obtained if B1(t) is a rect function, with the form: ( ) ⎪⎩ ⎪⎨ ⎧ > ≤ = 0, 1 2 , 1 2 1 1 t B t B t [1.34] Then the solution is M (t) = 0 xr M (t) M ( B t) yr z 1 = 0 sin −γ [1.35] M (t) M ( B t) zr z 1 = 0 cos −γ The solution indicates that the magnetization, initially pointing the z-direction, is rotating in the yr-zr plane about the applied B1 field at frequency 1 1 ω = −γB . For a pulse length of τ, the magnetization will rotate through an angle of α ω τ 1 = . Of course, in the laboratory frame the magnetization would also be precessing about the z-axis as is it being tipped into the transverse plane. At the conclusion of an applied B1 pulse, the magnetization will be some mix of longitudinal, Mzr(0), and transverse components, Mxyr(0), where Mxyr = Mxr + iMyr.. The 15 behavior of these vectors can still be described by Equation [1.28], the Bloch equation in the rotating frame, where Beff is now zero. This allows the Bloch equation to be broken down into two ordinary differential equations: 1 0 T M M dt dM zr zr z − = − r [1.36] 2 T M dt dM xyr = − xyr r [1.37] The solution in rotating frame is given by ( ) ( ( )) ( ) ( ) 1 1 M t M 0 1 exp t T M 0 exp t T zr z zr = − − + − [1.38] ( ) ( ) ( ) 2 M t M 0 exp t T xyr xyr = − [1.39] And the solution is laboratory frame is ( ) ( ( )) ( ) ( ) 1 1 M t M 0 1 exp t T M 0 exp t T z z z = − − + − [1.40] M (t) M ( ) ( t T ) ( i t) xy xy 2 0 = 0 exp − exp − ω [1.41] The solutions indicate that as the magnetization is precessing in the laboratory frame at the Larmor frequency, the transverse magnetization is decaying with time constant T2, while the longitudinal magnetization is regrowing to its initial magnitude, Mz 0, with time constant T1. Figure 1.2 shows a plot of the transverse and longitudinal components of the magnetization after a 90° RF pulse. It is the transverse magnetization that gives rise to the MR signal and must be sampled before it decays to zero. 16 Figure 1.2. Relaxation behavior for a spin system with T2 = 50 ms and T1 = 600 ms that has been tipped into the transverse plane with a 90° RF pulse. The transverse component decays to zero while the longitudinal component regrows to its original value. 1.5. Signal Detection Thus far it has been shown that a sample placed in a strong external magnetic field will become magnetized along the direction of the main field and that application of an additional magnetic field orthogonal to the main field and oscillating at the proper frequency will rotate the magnetization into the transverse plane whereupon it precesses about the main field at the Larmor frequency. This rotating magnetic field will create a detectable signal based on the principles of Faraday's law. According to Faraday's law, an EMF will be induced in a closed circuit that is equal to the negative of the time rate of change of the magnetic flux, Φ, through that closed circuit. 17 ( ) ( ) dt V t d t Φ = − [1.42] (t) B (r ) M(r t)dτ Φ = ∫ p ⋅ , r r r r [1.43] where Bp is the magnitude of the magnetic field that would be created by a unit current in the closed circuit at location p. The MR signal is detected by placing a conducting loop coil that resonates at the Larmor frequency near the sample such that the precessing transverse magnetization creates a changing flux through the coil, thus inducing a voltage. The signal voltage is calculated by inserting Equation [1.43] into Equation [1.42] and assuming that the z-component of the magnetization is changing slowly enough to be negligible. The signal voltage then becomes ( ) ( ) ( ) ( ) ( ) ∫ ⎥⎦ ⎤ ⎢⎣ ⎡ ∂ ∂ + ∂ ∂ = − dτ t M r t B r t M r t V t B r y p y x p x , , , , r r r r [1.44] The final form of the signal expression is obtained by using Equation [1.41] for the values of Mx and My, demodulating the signal using a reference signal of ( t) 0 2cos ω , low-pass filtering the signal, and assuming that Bp is homogenous over the region of interest: S(t) M (r ) ( iω (r )t) ( t T )dτ = ∫ xy − − 2 ,0 exp exp r r [1.45] where ω is now the deviation from the Larmor frequency and therefore represents the precessional frequency of the magnetization in the rotating frame. The most basic type of signal measured in an MR experiment is known as a free induction decay (FID) signal. An FID signal is one that follows from RF excitation with no other manipulation done on the magnetization. For a sample made up of a spin isochromat in a homogeneous magnetic field, the FID signal oscillates in time and decays 18 in magnitude as Equation [1.45] indicates. Such an ideal situation is rarely found in an MR experiment and in practice most samples are composed of spins with varying resonant frequencies due to their differing chemical environments and the fact that main magnetic field is not exactly homogeneous. To account for magnetic field inhomogeneities, ΔB, a new decay time constant, T2 *, is introduced: B T T = +γΔ 2 * 2 1 1 [1.46] The T2 * time constant is shorter than T2 due to the variations in the magnetic field causing individual spins to precess at different frequencies and therefore dephasing the net magnetization. Figure 1.3 gives examples of an FID signal with T2 * decay for a single Figure 1.3. FID signals. A) FID signal from a single spin isochromat. B) FID signal from a system with two different spin isochromats. A) B) 19 isochromat and an FID signal from two different spin isochromats with slightly different Larmor frequencies. The loss of signal due to T2 * decay cannot be avoided during an FID. However, some of the decayed signal can be recovered through the use of a signal detection technique known as spin echo acquisition. In this method, two RF pulses are played out with some delay time, τ, in between them in order to refocus the dephasing spins. For optimal signal generation, 90° and 180° RF pulses are used, although any combination will produce the same effect on a smaller scale. The refocusing effect of the 90° - τ - 180° RF pulse sequence is described in Figure 1.4. Figure 1.4. Schematic showing the behavior of the magnetization during the spin echo signal acquisition process. 20 A 90° RF pulse along the negative y-axis will tip the magnetization onto the x-axis (Figure 1.4A). The spins begin to precess, albeit at slightly different frequencies due to their different local magnetic fields, and therefore dephase over the delay time, τ. The 180° RF pulse flips the spins about the x-axis in the transverse plane (Figure 1.4C). Now that the order of the dephasing spins has been reversed, their differing frequencies will cause them to begin to rephrase, with total alignment achieved after the same delay time, τ (Figure 1.4D). The magnetization is refocused, although slightly diminished in magnitude due to irreversible T2 effects. A timing diagram for the spin echo signal is shown in Figure 1.5. Figure 1.5. Timing diagram for the spin echo acquisition, showing the signal decaying at a rate determined by T2 * and being refocused to an amplitude determined by T2. 21 Another common technique for creating a signal echo is the gradient echo acquisition. The timing diagram for this method is shown in Figure 1.6. After the RF pulse has created some amount of transverse magnetization, a magnetic field gradient is applied for some time, τ. The gradient field is aligned with the main field, but its magnitude varies linearly in one or more directions: B (G x G y G z)k G x y z = + + ˆ r [1.47] Figure 1.6. Timing diagram for the gradient echo showing the signal decaying at a rate faster than T2 * and being refocused to an amplitude determined by T2 *. 22 The application of this spatially varying field will create additional small differences in the local fields that determine the precessional frequency of the individual spins. This will cause the magnetization to dephase even faster than the T2 * rate. After the delay time, the gradient field is reversed and the signal echo will be formed at 2τ after the RF pulse. Note that the peak of the echo is determined by T2 * decay, not T2 decay as with the spin echo. 1.6. Signal Localization The theory presented thus far has demonstrated how the magnetic resonance phenomenon is used to generate a signal from within a sample and then how that signal is detected. These methods form the basis of both NMR and MRI experiments. In this section, the MRI theory diverges. In order to form an image, the detection scheme must be able to localize where the signal originated from in the sample. Up until now, no spatial information has been included in the MR signal. Signal localization is accomplished through the use of imaging magnetic field gradients. Just as with the gradient fields used in the gradient echo acquisition technique, the spatially varying magnetic fields impart different resonant frequencies to the spins in the sample based on their location. If the spatial dependence of the gradients is known, then the different spins can be localized and an image formed. 1.6.1. Frequency Encoding To frequency encode the MR signal, a linear gradient is applied in one direction, for example the x-direction, and the total magnetic field has the form 23 B (B G x)k x ˆ 0 = + r [1.48] which creates a spatially dependent resonant frequency of the form (x) G x x ω = ω +γ 0 [1.49] Consider a one dimensional object with a distribution of spins that is described by ρ(x). Using Equation [1.45] and neglecting the transverse relaxation term for the moment, the signal coming from the small region over dx will be dS(t) (x)dx ( i( G x)t) x = ρ − ω +γ 0 exp [1.50] The total signal will be a summation of all infinitesimal regions, which after demodulation of the carrier frequency ω0 has the form S(t) = ∫ (x) (− i G xt)dx x ρ exp γ [1.51] This signal is now spatially encoded along the direction of the applied gradient. However, this type of encoding alone does not suffice for 3-D objects, as each spin within the same y-z plane will have the same Larmor frequency and there will be no way to differentiate the locations of the different spins within those planes. To do this, another type of encoding, known as phase encoding, is employed. 1.6.2. Phase Encoding Phase encoding is based on the same principle as frequency encoding with a slightly different implementation. To do the phase encoding, a linear magnetic field gradient is applied in a direction orthogonal to the frequency encoding direction, say the y-direction, but only for a brief time period, PE τ . Once this gradient is turned off, the signal from a small 1-D region will look like 24 dS(t) (y)dy ( i G y ) ( i t) y PE 0 = ρ exp − γ τ exp − ω [1.52] Now, the different spins along the y-direction will have accumulated a particular phase offset based on their location, y. After again removing the carrier frequency, the total signal equation for combined frequency and phase encoding will look like S(t) = ∫ (x y) (− i G xt) (− i G y )dxdy x y PE ρ , exp γ exp γ τ [1.53] When combined with frequency encoding, every spin in a 2-D plane can be localized due to its unique frequency and phase offset combination. 1.6.3. 2-D and 3-D Imaging The task still remains of localizing the signal in the third dimension, the z-direction. This is accomplished in one of two different ways, with the different techniques referred to as 2-D and 3-D imaging. Each method starts by selectively exciting only the spins that fall within a predetermined region along the z-direction. For 2-D imaging, this region is a few millimeters thick and the excitation is known as slice selection. For 3-D imaging, the region is larger and the excitation is referred to as slab selection. In either case, the problem is how to design an RF excitation that only tips the desired spins into the transverse plane. As demonstrated in section 1.4, a B1 pulse will affect any spin that is precessing at the same oscillation frequency as the pulse. The solution, therefore, is to apply a magnetic field gradient in the z-direction during RF excitation in order to have the spins' frequencies vary linearly along the z-axis: (z) G z z ω = ω +γ 0 [1.54] 25 Exciting spins over a spatial width of Δz, centered at z0, requires exciting only the spins within the frequency range Δω, centered at ω0, where G z z Δω = γ Δ [1.55] For small tip angles, it can be shown that the desired frequency selection function, ρ(ω), is related to the RF pulse waveform, B1(t), by the Fourier transform (see Liang and Lauterbur for a full treatment of the derivation). For uniform excitation over a finite region, the desired frequency selection will be given by ( ) ⎟⎠ ⎞ ⎜⎝ ⎛ Δ − = Π ω ω ω ρ ω 0 [1.56] where Π represents the rect function. The Fourier transform of a rect function is a sinc function, giving the following form for the RF pulse: B (t) = ΔwΣ(Δωt)exp(− iωt) 1 [1.57] where Σ represents the sinc function. For 2-D imaging, a sinc RF pulse played out in the presence of a strong z-gradient will selectively excite only a thin 2-D plane of spins within the sample. The spins within the 2-D slice can then be spatially encoded using frequency and phase encoding as described above. To image over a volume, multiple 2-D slices are used. For 3-D imaging, a weaker z-gradient is used to excite a thicker slab of spins. In order to get spatial localization within this excited slab, the signal is phase encoded in the z-direction in addition to frequency encoding in the x-direction and phase encoding in the y-direction. 26 1.6.4. A k-space Interpretation When acquiring data to form an entire image, it useful to think of the data acquisition process in terms of a k-space formulation which makes an explicit connection between the spatial encoding and the Fourier transform. Upon examination, it can be seen that the spatially encoded signal equation given in Equation [1.53] has the form of a Fourier transform. This is made more obvious by introducing the variables: k γG t 2π x x = [1.58] γ τ 2π y y PE k = G [1.59] Then Equation [1.53] becomes S(k k )= ∫ (x y) (− i (k x + k y))dxdy x y x y , ρ , exp 2π [1.60] and the Fourier relationship between the encoded signal and the spin density function becomes explicit. With this notation, k-space is defined by the variables kx and ky. In order to obtain an image, enough data must be acquired to sufficiently cover k-space. A timing diagram for covering 2-D k-space in a Cartesian fashion using gradient echo imaging is shown in Figure 1.7, where the data are acquired over N excitation cycles. The nth RF pulse is played out in a sinc waveform while a z-gradient is simultaneously applied for slice selection. A rephrasing z-gradient is applied to counteract the dephasing that occurs during excitation. Next, the nth phase encoding gradient is applied in the y-direction with strength Gy,n and the first lobe of a frequency encoding gradient is applied in the x-direction with strength -Gx. The effect of these two pre-phasing gradients is to move in k-space from the origin to γ τ 2π x x PE k = − G [1.61] 27 γ τ 2π y y,n PE k = G [1.62] Next, the phase encoding gradient is turned off and the direction of the frequency encoding gradient is reversed. Sampling commences using an analog to digital converter (ADC), and data is acquired over one line in k-space as the gradient echo signal is formed. During this data readout phase, the ky location remains constant because no y-gradient is on, while the kx location is moving from left to right due to the constant x-gradient. This process is repeated many times with different strengths for the phase encoding gradient such that the readout lines cover the entire k-space. A) B) Figure 1.7. Pulse sequence timing diagram. A) Timing diagram for gradient echo sequence showing the RF pulses, all gradient waveforms, and ADC signal acquisition. B) Cartesian k-space path determined by x- and y-gradients. 28 1.7. Image Reconstruction With the data fully sampled, the final step in the MRI experiment is to reconstruct the image. Consider a one dimensional object, I(x), and its Fourier transform samples, S(kn): S(k ) = ∫ I (x) (− i k x)dx n n exp 2π [1.63] The reconstruction problem concerns how to obtain an accurate representation of I(x) based on the limited number of data points contained in S(kn). For uniform sampling of k-space, kn can be replaced with nΔk, with n = …, -2, -1, 0, 1, 2, … Equation [1.63] becomes: S(n) = ∫ I (x)exp(− i2πnΔkx)dx [1.64] To extract an expression for I(x), one must make use of the fact that a finite object replicated ad infinitum at a period of 1/Δk can be expressed as a Fourier series according to the formula: Σ Σ ( ) ( ) ∞ =−∞ ∞ =−∞ Δ = ⎟⎠ ⎞ ⎜⎝ ⎛ Δ − Δ n n S n i n kx k I x n k 1 exp 2π [1.65] For an object limited in space to |x| < Wx/2, the periodic extension will not have any overlap if Δk < 1/Wx. In this case, one period of the object can be selected such that ( ) ( ) ( ) Σ∞ =−∞ = Δ Δ n I x k S n exp i2πn kx [1.66] This relation implies that infinite sampling is necessary to accurately reconstruct an image of the object, I(x), which is clearly not feasible for any real MRI experiment. In practice, only N samples are acquired and Equation [1.66] is rewritten as 29 ( ) Σ ( ) ( ) − =− = Δ Δ / 2 1 / 2 exp 2 N n N I x k S n i πn kx [1.67] This truncation of sampling data leads to some inaccuracies in the reconstruction due to the Gibbs ringing phenomenon, but it is a trade off that must be made. I(x) is a continuous function, which is also not feasible to create in a real setting where digitization is required, and therefore the final step is to formulate the reconstructed image as N discrete values obtained from the N sampled data points. Using the Nyquist sampling criteria to set the pixel size as N k x Δ Δ = 1 [1.68] Equation [1.67] can be rewritten in discrete form as: ( ) Σ ( ) ( ) − =− = Δ / 2 1 / 2 exp 2 N n N I m k S n i πnM N [1.69] This is the final discrete imaging equation for one dimension. It is easily extended to two or three dimensions. Figure 1.8 shows an example of a fully sampled k-space data set and the reconstructed image that it forms. Figure 1.8. MR k-space and image. A) Fully sampled k-space data set from a gradient echo sequence for a 128 x 128 image. B) Magnitude image created from k-space data A) B) CHAPTER 2 MR GUIDED HIGH INTENSITY FOCUSED ULTRASOUND High intensity focused ultrasound (HIFU) can heat regions of tissue inside the body selectively and noninvasively without the use of ionizing radiation. The low impact nature of the technique means that HIFU treatments have the potential to have a significant impact on how certain types of cancers are treated. They could be used to treat cancerous tumors without the drawbacks of invasive surgery, they would not have the complications associated with radiation, and they could potentially reach tumors that were previously considered to be inoperable. Already, HIFU has been used to treat both benign and malignant tumors in a variety of sites, including uterus, liver, kidney, pancreas, bone, breast, prostate, and brain [10-26]. Because of its noninvasive nature, safety and efficacy considerations require a method to monitor and control the heating during a HIFU ablation procedure. Many investigators have turned to MRI for this task, referring to the combined process as MR-guided HIFU (MRgHIFU) [27-34]. This chapter will present the background and basics of MRgHIFU treatments and illustrate the necessity for accurate temperature mapping to properly control the process. 31 2.1. Ultrasound Physics Ultrasound is an acoustic pressure wave with a frequency above 20 kHz. For HIFU applications, the frequencies range from several hundred kHz to several MHz. Ultrasound transducers work by applying an alternating voltage across a piezoelectric material. The applied voltage deforms the material in a periodic fashion, leading to the creation of pressure waves. As the ultrasound wave travels from the transducer to the target and beyond, its path will be affected by the physical properties of the conducting media, including the speed of sound (c), impedance (Z), and absorption (μ). At all medium interfaces some percentage of the incoming ultrasound will be transmitted and the rest reflected, with the percentage being determined by the acoustic impedances of the two materials. Equation [2.1] gives the formula for the fraction of energy that will be reflected when ultrasound passes between media with impedances Z1 and Z2. 2 2 1 2 1 ⎟ ⎟⎠ ⎞ ⎜ ⎜⎝ ⎛ + − = Z Z R Z Z [2.1] Tissue/air interfaces have a very large acoustic impedance mismatch and therefore ultrasound needs a coupling medium in order to travel from the transducer to the subject. For a large impedance mismatch, such as at tissue/air or tissue/bone interfaces, the majority of the ultrasound will be reflected. The ultrasound transmitted through a medium interface will be deflected from its original path due to refraction effects that depend on the speeds of sound in the two materials. The formula for the refraction angle is Snell's Law, the same equation that governs the refraction properties of light: 32 2 2 1 sin 1 sin c c θ θ = [2.2] where θ1 and θ2 are the angle of incidence and the angle of refraction, and c1 and c2 are the respective speeds of sound. A schematic of an ultrasound wave incident on an interface is displayed in Figure 2.1, showing the incident wave, VI, the reflected wave, VR, and the transmitted, refracted wave, VT. As the ultrasound wave propagates it will also lose energy through absorption. Different tissue types have different absorption properties, with muscle being approximately twice as absorptive as fat. Water is not highly absorbing, and has an acoustic impedance close to that of tissue, and is therefore often used as the coupling medium between the transducer and the subject. All of these factors must be taken into account when determining how the ultrasound will travel through the subject and where the ultrasound energy will be deposited. Figure 2.1. An ultrasound wave, VI, incident on an interface where the two different media has different impedances, Z, speeds of sound, c, and absorptions, μ. Some of the ultrasound wave will be reflected, and the transmitted portion will be refracted. 33 At low power intensities, ultrasound can travel harmlessly through tissue. Diagnostic ultrasound emits plane waves of acoustic radiation and forms an image based on the timing and power of the reflected waves that come back to the transducer. The absorption is low enough that no significant tissue heating occurs. HIFU uses the same ultrasound waves as diagnostic imaging, but at a higher intensity and focuses them such that they interfere constructively within a small volume. The resulting large deposition of energy into this area will cause heating of the tissue, which can lead to necrosis and cell death. To do the focusing, the transducer can be a single element transducer made up of one large curved piezoelectric piece, or a phased array transducer made up of many smaller elements. For the phased array transducers, the phases and amplitudes of the individual elements can be specified separately, allowing electronic steering and shaping of the focus. Temperature rises are induced by HIFU as the tissue absorbs the ultrasound and the acoustic energy is converted into thermal energy. The oscillating pressure waves of the ultrasound couple to particles in the tissue, causing them to vibrate and rotate. This microscopic mechanical motion is what causes the heating of tissue. It has been shown that the deposition of thermal energy causes tissue necrosis, but the physical process of how the heating induces cell death is not well understood. One metric for estimating the amount of damage caused by heating is thermal dose. It is an empirical metric that accounts for the cumulative effects of heating over time and is defined as [35]: Σ= = +Δ − Δ N n n D R T Tn t 1 ( 43) 43 0 [2.3] 34 where D43 is the thermal dose that would be equivalent to continuous heat application at 43°C (Cumulative Equivalent Minutes or CEM), T0 is the baseline temperature and is set to 37°C in this work, ΔTn is the measured temperature change at time n (°C), Δtn is the time step between measurements in minutes, and R is an empirically determined constant set to 4 for ΔTn less than 6°C and to 2 for ΔTn greater than 6°C. The amount of thermal dose necessary to cause cell death is still a topic of experimental research, but most investigators agree that reaching 240 CEM will ensure necrosis [36]. 2.2. HIFU Challenges A typical MRgHIFU tumor ablation procedure is carried out by applying numerous sonications with a cooling time in between them. Each sonication heats a small volume within the tumor and is delivered until an end-point temperature rise or accumulated thermal dose is reached, with the MR scanner used for measuring temperature changes. After the cooling time, a new location is sonicated. This process is repeated until the entire targeted area in and around the tumor is treated to satisfaction. While the basic steps of the procedure are simple enough, there are a number of complicating factors that make MRgHIFU treatments difficult to carry out in practice. Figure 2.2 shows a schematic of a HIFU sonication with the numbers indicating areas where problems may arise. 1) At the center of the focal zone the temperatures are rising most rapidly and the spatial temperature gradients are the largest. It is in this area where real time temperature monitoring is most critical to control the treatment process. 2) Even though it is before the focal zone, some energy will be deposited in the near field of the ultrasound. The acoustic window is wider at this point, and ultrasound will pass 35 Figure 2.2. Schematic of the set up for a HIFU tumor ablation. Ultrasound energy is used to selectively heat the tumor to the point of necrosis. A number of potential challenges exist, including 1) rapid temperature rises at the focus, 2) near field heating, 3) beam distortion due to fat, 4) heating at tissue/bone interfaces, and 5) motion of the target organ. through the same tissue areas for most sonications. If too little cooling time is prescribed between sonications, thermal build up could occur and cause unwanted damage to healthy tissue. 3) The presence of fatty tissue in the path of the ultrasound beam will have two effects. First, MR temperature measurements are more difficult in fat and therefore it will be harder to monitor heating in that tissue. Second, the different ultrasound tissue properties of fat will distort the beam path, causing the focus to blur or change location from where it was intended to be. 4) As noted above, tissue/bone interfaces strongly reflect ultrasound. This could potentially induce heating at the bone surface, even if the 36 interface is far from the focal zone. 5) Many organs of interest for MRgHIFU treatments are in constant motion due to respiration. Not only does this create a moving target for the ultrasound, it also complicates the MR temperature imaging. 2.3. Need for Temperature Monitoring Many of the challenges noted above require accurate and rapid temperature mapping to be overcome. The temperature mapping in the focal zone must be accurate to ensure that the end-point temperature or thermal dose is properly reached. As can be seen from Equation [2.3], a one degree error in temperature leads to a factor of two error in the calculated thermal dose. The temperature mapping of the focal zone must also have high spatial resolution in order to accurately measure the underlying temperature distribution. The full width at half maximum (FWHM) dimensions for the power deposition of a typical HIFU transducer are approximately 2 x 2 x 8 mm, leading to temperature rises that vary spatially over a small volume. If the temperature mapping voxels are too large, averaging effects will lead to an underestimation of the true temperature. The temperature maps must also be acquired rapidly. HIFU heating can induce temperature rises as fast as 2 °C/sec in the focal zone, requiring image acquisition rates of one per second or greater in order to properly control the energy deposition. Additionally, the image acquisition must be rapid enough to "freeze" any motion that may be occurring. Finally, the temperature mapping must cover a large volume. As indicated in Figure 2.2, undesirable heating can occur in the near or far fields of the ultrasound path. These regions must be monitored to ensure that healthy tissue is not inadvertently damaged. These temperature imaging requirements of high spatial 37 resolution, high temporal resolution, and large volume coverage constitute a challenging list of demands to meet. MRI is a promising modality for meeting these needs. MR temperature imaging can measure temperature changes in a number of different soft tissue types with an accuracy of ±1°C [20]. However, the imaging parameters of spatial resolution, temporal resolution, and volume coverage are self-conflicting in MRI and one can only be improved at the expense of another. Currently, the image acquisition and reconstruction process is efficient enough to meet the requirements for one of these factors, but not for all three simultaneously. CHAPTER 3 MRI THERMOMETRY This chapter presents the basics of MRI thermometry, outlining the different temperature measurement techniques, current challenges, and most recent advancements in the field. A number of different parameters measurable by MRI have a temperature dependence, making MR temperature mapping feasible. The three most commonly measured parameters are the longitudinal relaxation time (T1) [37-40], the diffusion coefficient (D) [41-42], and the proton resonant frequency (PRF) [43-45]. The temperature related properties of the three parameters vary considerably due to their different physical basis and the different methods utilized to measure them. For additional treatments of the material that follows, see the excellent review papers by Quesson et al [46] and Rieke et al. [47]. 3.1. Temperature Dependence of the T1 Relaxation Time of Water Protons The mechanism for the longitudinal relaxation time, T1, is the interaction of the water protons with the thermal motion of their surrounding molecules. These randomly tumbling molecules create a randomly varying local magnetic field at the site of the water proton. This local field will have a frequency spectrum density, J(ω), whose width is inversely related to the time scale of the fluctuations. The transitions necessary for 39 relaxation to thermal equilibrium can only occur if there are magnetic field fluctuations at the Larmor frequency, ω0. Frequency spectrum densities that have larger values at the Larmor frequency will result in shorter T1 values. An example of a randomly varying local magnetic field and three different examples of frequency spectrum density plots are shown in Figure 3.1. The rate of motion of these randomly moving molecules is described using their correlation time, τc. More rapid motion results in a lower correlation time. It can be shown that the T1 relaxation time is related to the correlation time according to the formula [9]: ⎟ ⎟ ⎠ ⎞ ⎜ ⎜ ⎝ ⎛ + + + ⎟ ⎟⎠ ⎞ ⎜ ⎜⎝ ⎛ = 2 2 2 2 2 3 2 1 0 0 1 4 4 10 1 1 3 c c c c T r ω τ τ ω τ hγ τ [3.1] The correlation time, τc, is inversely proportional to temperature and for 3T experiments on water protons at 20°C, ω0·τc<<1. Under these conditions, Equation [3.1] reduces to T c ∝τ 1 1 [3.2] And therefore, T1 is directly proportional to temperature. T ∝ T 1 [3.3] A number of different investigators have shown experimentally that T1 does vary linearly with temperature over ranges of 15 °C - 55 °C. The measured temperature dependencies where all around 1%/°C, with reported values of 1.4%/°C in bovine muscle [38], 1 - 2%/°C in liver [48], and 1.0%/°C in fat [49]. Based on this experimental evidence, it is possible to estimate the relationship between T1 and temperature with a linear equation: 40 Figure 3.1. Local field effects T1. A) Randomly varying local magnetic field at the site of a water proton. B) Frequency spectrum density for fields varying at different rates. T (T ) T (T ) m T ref = + Δ 1 1 [3.4] where Tref is the reference temperature, m = dT1/dT, and ΔT is temperature change. T1 values at body temperature will differ for different tissues, and as the reported results indicate, the temperature dependencies will vary as well. Therefore, the slope, m, must be determined empirically for each different tissue type. An experiment was performed to demonstrate the temperature dependence of T1. An agar phantom was heated to 45 °C 41 and allowed to cool in the MR scanner. Images were acquired during cooling for T1 measurements, and an alcohol thermometer was used to independently measure the temperature of the agar. A plot of the measured T1 against temperature is shown in Figure 3.2. A number of different strategies exist for measuring T1 directly. Inversion recovery techniques use a 180° RF pulse to invert the magnetization and then sample the signal after the z-magnetization has been allowed to recover for a time, TI. This is done for a number of different TI values, and the acquired signals will have a dependence that Figure 3.2. Plot of T1 values against temperature. The temperature dependence of T1 shows good linearity over a small range of temperatures around body temperature. (°C) 42 can be derived from Equation [1.37] which describes the z-magnetization recovery to equilibrium: ( ) ( ( )) 1 S t = A 1− 2exp − TI T [3.5] The value for T1 can be extracted through any parameter fitting method, using the known values for TI and the measured values for S(t). The imaging time for such a method is far too long to be of use for HIFU applications and therefore investigators have proposed techniques for more rapid measurement of T1. A Look-Locker segmented EPI sequence was implemented to measure T1 values in a few seconds [50]. Another method, referred to as driven equilibrium single pulse observation of T1 (DESPOT1), measures T1 values by acquiring gradient echo signals at multiple different flip angles [51-52]. The spoiled gradient echo signal is a function of the net magnetization, M0, the flip angle, α, the sequence repetition time, TR, and T1: 1 1 0 1 cos sin 1 E S M E α α − − = [3.6] ( ) 1 1 E = exp −TR T [3.7] Equation [3.6] can be linearized as follows: ( ) 1 0 1 1 sin tan S = E S + M − E α α [3.8] When two or more acquisitions are made at different flip angles, the signals can be plotted according to Equation [3.8] and the fitted slope will give an estimate of T1 based on Equation [3.7]. As an alternative to directly measuring the T1 value, changes in a T1-dependent magnitude signal can be measured and correlated to temperature changes. If the T1 43 temperature dependence is assumed to be linear, then the signal equation of Equations [3.6] and [3.7] can be rewritten with the expression for T1(T) of Equation [3.4] inserted. Imaging parameters should be chosen to maximize the relative sensitivity of the signal magnitude changes with respect to temperature, dS/SdT. This expression is given by [46]: ( ) ref( )( ) ref T T E E T ATR E SdT dS 1 1 1 cos 1 cos 1 1 2 1 1 − − − − = − α α [3.9] The second term on right in Equation [3.9] comes from the fact that the equilibrium magnetization, M0, is temperature dependent itself, as can be seen by referring back to Equation [1.18]. This approach of measuring T1-dependent signal changes has the advantage that a standard gradient echo or spin echo sequence can be used, which can have adequate SNR and temporal resolution for a limited amount of volume coverage. Even if the problem of temporal resolution is overcome, a number of challenges remain for using T1 changes to map temperature changes. As noted above, the temperature dependence will vary from tissue to tissue, and empirical calibration will be necessary for each tissue type in the region of interest. Additionally, the T1 relaxation process is based on the tumbling rates of neighboring water molecules. When HIFU heating induces cell death or tissue coagulation, it is very possible that these mechanisms could change nonlinearly. Indeed, different investigators have measured nonlinear T1 changes in tissue during heating [53]. When using techniques that rely on measuring changes in the signal magnitude, care must be taken to account for all other factors that could affect the signal strength. The temperature dependence of the equilibrium magnetization has already been mentioned, but temperature related phase changes could also play a role (see section 3.3 for a description of how image phase changes with 44 temperature). If a large temperature gradient exists across an imaging voxel, the temperature-dependent phases will disperse and cause a drop in signal strength. The use of small voxel sizes and/or a spin-echo acquisition will mitigate these effects. Finally, the presence of two different tissue types within the same voxel will cause complications due to their differing T1 temperature dependencies. If one of the tissue types is fat, suppression techniques should be used. 3.2. Temperature Dependence of the Diffusion Coefficient of Water Water molecules at body temperature have thermal energy that keeps them in constant motion on a microscopic scale. The movements of an ensemble of water molecules due to this thermal energy are described by the theory of thermal Brownian motion. The randomly interacting particles will spread at a rate determined by their diffusion coefficient, D. Similar to T1, the dependence of D on temperature can be approximated by an energy activation equation [42]: D(T ) ( E kT ) a,D ∝ exp − [3.10] where Ea,D is the activation energy of the diffusion process, k is Boltzmann's constant and T is the absolute temperature. The activation energy for diffusion has been reported to be significantly higher than the activation energy for T1. While this gives D(T) a better sensitivity with respect to temperature changes, the linearity of D as a function of T is not as good over ranges around body temperature. The same cooling experiment as described in 3.1 was performed to demonstrate the temperature dependence of the diffusion coefficient. Figure 3.3 shows a plot of the apparent diffusion coefficient measured using MR against the independently measured temperature. 45 Figure 3.3. Apparent diffusion coefficient (ADC) plotted against temperature. The diffusion coefficient also shows reasonable linearity in the range of temperatures around body temperature. Although the molecular displacements due to Brownian motion are very small, they can still be detected by MR using a pair of strong, opposing gradients played out in succession after RF excitation. The hydrogen spins will acquire a certain phase due to the first gradient, and will be fully refocused by the second gradient if they are stationary. Spins that are moving in the direction of the applied gradients due to diffusion will not be fully refocused. This phase dispersion will result in a signal loss of e-bD, with b = −γ 2G2δ 2 (Δ −δ 3) [3.11] where G is gradient strength, δ is the duration of the gradient, and Δ is the time between the starts of the two gradients. In this way the diffusion coefficient can be measured and temperature information extracted. (°C) 46 Using the diffusion coefficient for measuring temperature changes poses a number of challenges. First, the temperature dependence of the diffusion constant will be tissue-specific and therefore needs to be calibrated for each new tissue type. Diffusion in fat is particularly slow and therefore the fat signal should be suppressed when necessary to avoid partial volume complications. Because the diffusion technique relies on detecting microscopic motion, the methods are very sensitive to macroscopic motion. Additionally, structures present in tissue, such as muscle fiber, force the water molecules to diffuse anisotropically. This directional dependence implies that either the full diffusion tensor or its trace should be measured, adding scan time to the process. As with the T1 measurement, coagulation of tissue is problematic for the measurement in that it leads to structural changes that could have major effects on the diffusion pathways. This will cause the temperature dependence to be nonlinear at temperatures in the range necessary for tissue necrosis. 3.3. Temperature Dependence of the Proton Resonance Frequency As described in Chapter 1, the resonant frequency of a hydrogen proton is determined by the local magnetic field. A number of factors can create deviations from the main field strength at the local level, however for the purpose of temperature measurement, the important factor is the chemical shift. The chemical shift is a change in a nucleus' resonant frequency due to its chemical environment. The external magnetic field interacts with the electrons in an atom to set up a current density. This current density creates a magnetic field at the site of the nucleus in the opposing direction of the main field, lowering the local magnetic field that the nucleus experiences, and therefore 47 also lowering its resonant frequency. For water protons the chemical shift is temperature dependent. The 1H nucleus in an H2O molecule, which is hydrogen bonded to another H2O molecule, is screened by an electron cloud. The hydrogen bonds between H2O molecules affect the extent to which the electron cloud screens the 1H nucleus. Increases in temperature can bend, stretch, or break these hydrogen bonds, causing the electron cloud to screen its nucleus more strongly. The additional screening lowers the local magnetic field, and therefore the resonant frequency will decrease with increasing temperature [44, 54]. A schematic of this process is shown in Figure 3.4. Figure 3.4. A schematic depicting the temperature dependence of the water proton resonance. As temperature increases, the hydrogen bonds between water molecules change such that the hydrogen proton is more strongly screened, leading to a lower resonant frequency. 48 It must be noted that a number of other elements and molecules have been shown to have a temperature dependent chemical shift. For example, studies performed on solid lead nitrate demonstrated the temperature dependence of the chemical shift of 207Pb [55], and the temperature dependence of the resonant frequency of 119Sn is commonly used is various compounds as an NMR "chemical shift thermometer" [56]. The physical mechanism behind these temperature dependencies is not well understood for these cases, but it is clearly not due to hydrogen bonds being broken. Some other cause must exist that may also be present in the case of the temperature dependence of the 1H chemical shift. The local magnetic field strength as a function of the main field, the temperature dependent chemical shift, σ(T), and other nontemperature related contributions, ΔB0, can be written as: ( ) ( ( )) 0 0 B T 1 T B B loc = +σ +Δ [3.12] The temperature dependence of the chemical shift has been shown to be linear over a wide range of temperatures around body temperature, giving: σ (T ) = αT [3.13] where α is the chemical shift coefficient and has been measured to be -1.03 ± 0.02 * 10-8 /°C [44]. The temperature dependent resonant frequency is detected in MRI using the phase of the acquired image: ( ) ( ) [( ( ) ) ] 0 0 T B T T T 1 T B B loc E E Φ = γ = γ +σ + Δ [3.14] where TE is the echo time of the sequence. The nontemperature-related changes in the local field that contribute to the phase must be removed and therefore a temperature 49 change is calculated based on the change in phase between images acquired at different times. The formula is derived from Equations [3.13] and [3.14] [43]: ( ) ( ) E ref B T T T T 0 αγ Φ − Φ Δ = [3.15] When taking the difference between the current phase image and the reference phase image, a complex subtraction is preferable to direct subtraction because of phase wrapping issues. The complex phase subtraction between two images, I1 and I2, is carried out as: ( ) ( ) ( ) ( ) ( ) ( ) ( ) ( ) 1 2 1 2 1 2 1 1 2 1 2 Re Re Im Im tan Re Im Re Im I I I I I I I I + − Φ − Φ = − [3.16] where Re and Im represent the real imaginary parts of the image. Figure 3.5 shows an example of an MR phase image and the phase difference between two images acquired at different times. The hot spot from HIFU heating can be seen in the phase difference image, although this phase difference has not yet been converted to a temperature difference. A gradient echo sequence must be used to acquire the phase information as a spin echo sequence will refocus the temperature-dependent phase contribution. When optimizing the parameters of the gradient echo sequence, the signal to noise ratio of the image phase, SNRΦ, should be considered. SNRΦ is defined as (Φ) Φ = Φ σ SNR [3.17] where σ represents the standard deviation. The phase accumulation will increase with TE according to Equation [3.14] and it can be shown that σ(Φ) is inversely proportion to the 50 Figure 3.5. Phase maps from MR images acquired before heating (reference) and during heating (current time frame). The difference between the phases is proportional to the temperature change. signal magnitude SNR, which decreases exponentially with TE, giving [46]: ( * ) 2 SNR T exp T T E E ∝ − Φ [3.18] To maximize this expression, TE should be chosen equal to T2 *. For most tissues of interest, T2 * is approximately 20 - 50 ms. Choosing a TE this long would greatly decrease the temporal resolution of the imaging, and therefore some intermediate value should be chosen that balances scan time and SNR issues. The chemical shift coefficient, α, appears to be independent of tissue type, provided the tissue is an aqueous-based one. A number of different investigators have calibrated the phase change as a function of temperature change in a variety of tissue types and most have measured α to be within -0.8 * 10-8 /°C to -1.2 * 10-8 /°C [36]. The temperature dependence of the PRF shift does not seem to be effected by tissue coagulation. The linearity range has been measured to extend up to 80 °C, well past the point of tissue necrosis. 51 While the different measurements of α do not correlate with different tissue types, this spread in the calibration results is still over 30%. It is therefore still advised to run a calibration experiment to determine the precise value of α for a given experimental set up. Figure 3.6 shows an example of such a calibration experiment where an agar phantom was heated uniformly and allowed to cool in the scanner while its temperature was being monitored with an alcohol thermometer. The slope of the line fit to the phase-temperature curve was -1.15 * 10-8 /°C. Figure 3.6. Calibration experiment results relating the phase change to the temperature change. The phase change is linear with a slope of -0.015 ppm/°C. m = -0.015 ppm/°C 52 3.4. Challenges with the PRF Method The PRF technique has many advantages over the other methods of measuring temperature changes with MRI and that is why it has become by far the most widely used technique. However, there are still a number of limitations and challenges involved when using the PRF method. The PRF method can measure only relative temperature changes, not an absolute temperature. This can be a problem when calculating accumulated thermal dose. As seen in Equation [2.3] from Chapter 2, the calculation of thermal dose uses a reference temperature which is assumed to be 37 °C. While this is a good approximation for body temperature, it may not be accurate across different organs or subjects. For example, locations in the breast may be a few degrees below regular body temperature. In the thermal dose calculation, each degree of error in the absolute temperature will lead to a factor of two error in the dose calculation. The PRF method does not work in adipose tissues. As noted above, the physical basis of the temperature dependent chemical shift of water protons is the hydrogen bonds between H2O molecules that give rise to the temperature dependent electron screening of the nucleus. The molecules of the various types of fat that make up adipose tissue do not have such bonds. Their 1H protons have a chemical shift due to the electrons surrounding them, but it is not temperature dependent. In regions of the body that only contain fat, temperature measurements with the PRF technique are not possible. PRF temperature measurements from a voxel containing both fat and water will have some error because the fat still contributes signal, but its phase does not change with temperature. Therefore, 53 fat suppression should be employed when using the PRF method to avoid partial volume effects. As indicated by Equation [3.14], the phase of an image will be affected by anything that affects the local magnetic field. This includes a number of phenomena that are not related to the temperature dependent chemical shift. By taking the difference of two images, all of these spurious effects are subtracted out, and any remaining phase differences are assumed to be due to temperature changes. However, it is possible that one or more of these field-changing phenomena could change in between scans, making the local field different for the two different time frames. In this case, the change in the local field will create a phase difference that will be erroneously interpreted by the PRF method as being due to a temperature change. For example, the strength of the main magnetic slowly drifts over time, changing the resonant frequency on the order of a few Hertz per minute. Over the course of heating experiment that lasts several minutes, this will induce temperature errors of a few degrees C. Another nontemperature related phase change could come from interscan motion. Motion, whether it is inside or outside the imaging field of view, will change the bulk susceptibility within the field of view, which in turn affects the local magnetic field. Sometimes this motion is due to inadvertent subject movement, but often it is due to unavoidable physiological motion like breathing or heart beats. When imaging the breast, for example, motion of the chest wall from breathing can cause temperature errors up to 5 °C - 10 °C. Figure 3.7 shows an example of these temperature errors. Figure 3.7A displays a magnitude image of a breast that was acquired with the breast held still, during free breathing, and without any external heating. Ten 54 Figure 3.7. Motion related temperature errors. A) Magnitude image of the breast, imaged during free breathing without heating. B) Temperature plot of one voxel over time showing errors due to susceptibility changes caused by respiratory motion. images were acquired every 8.3 seconds and the phase information was used to create temperature maps. Figure 3.7B plots the measured temperatures from one voxel within the glandular tissue over the course of 75 seconds, showing the periodic nature of the errors due to breathing. In addition to the changes in bulk susceptibility that motion causes, it is also problematic for the PRF method because of registration issues. If inter-scan motion occurs for tissue inside the field of view, the images will be misregistered and the phase difference will be calculated between mismatched voxels, leading to temperature measurement errors. Finally, the spatial and temporal resolutions of the PRF temperature maps are not adequate for properly monitoring the temperatures changes induced by HIFU heating over the entire region of interest. For all MR imaging sequences, trade-offs must be made amongst SNR, spatial resolution, temporal resolution, and coverage. In the standard RF-spoiled gradient echo sequence used for the PRF method, the user can 55 modify the imaging matrix size, the field of view, the sequence repetition time (TR), TE, and whether to do 2-D or 3-D imaging (amongst other things). A typical protocol might have the following parameters: 2-D imaging, five slices acquired in an interleaved fashion, 128x128x5 imaging matrix, TR = 65 ms, TE = 8 ms, 2.0 x 2.0 mm in-plane resolution, and 3.0 mm thick slices with 1.0 mm gap between slices. This sequence protocol would provide adequate SNR, but the spatial resolution is suboptimal, there would not be coverage in the near or far field of the ultrasound, and the scan time of 8.3 seconds is too long to capture the rapidly rising temperatures in the focal zone. The spatial resolution could be improved by choosing 1.0 x 1.0 mm in-plane resolution, but only at the expense of worse SNR, less coverage or longer scan times. Similarly, the scan time could be shortened by reducing the TR, but this would require sacrificing SNR, spatial resolution, or the number of slices acquired. Similar trade-offs exist when trying to improve any one of the imaging parameters. To accurately capture the spatial and temporal features of a HIFU-induced temperature distribution, the MR images would ideally be acquired at 1.0 x 1.0 x 3.0 mm spatial resolution, 1.0 second temporal resolution, and cover a region of approximately 256 x 128 x 72 mm. Using a conventional 3-D gradient echo sequence, it would take approximately 60 seconds to get this spatial resolution and volume coverage. 3.5. Overcoming MR Temperature Imaging Limitations Each of the challenges discussed above are well known and have been worked on by many different investigators for a number of years. While none are fully solved, a 56 number of innovations have made significant improvements to the MR temperature imaging landscape. Two approaches have been presented for acquiring absolute temperature measurements. The first approach involves using a temperature sensitive contrast agent, where a gadolinium or manganese-based compound is enclosed by a phospholipid membrane [57-58]. The membrane is not permeable to water below a critical phase transition temperature, Tc, which is accurately known. Once the critical temperature is reached, water permeates the membrane and interacts with the contrast agent, causing a large decrease in the T1 relaxation time, which is detectable in the MR signal. Absolute temperature can therefore be measured, but only at Tc. The second approach uses a spectroscopic method where the water signal is referenced against another, temperature independent signal such as lipid or N-acetyl-aspartate (NAA) [59-60]. Each of these reference signals has a chemical shift that is different from water. The water signal and the reference signal will each be equally affected by nontemperature related changes in the local magnetic field, and therefore the frequency separation between the two signals will be temperature dependent. If the frequency separation between water and the reference signal is known for a specific temperature, then absolute temperature imaging can be done as the water peak changes relatively to the reference peak with temperature. One drawback is that each species must be present within a voxel in order to do temperature mapping. This is not always the case, and even when it is, the strength of the reference signal may be weak enough that it requires large voxels or long imaging times to be detectable. 57 The problem of how to best measure temperature changes in voxels that contain fat is also an area of on going research. The most promising approaches to measuring temperature changes in voxels that only contain adipose tissue have utilized T1. As noted above, methyl and methylene molecules do not have a temperature dependent chemical shift which rules out the PRF method, and diffusion in fatty tissues is very slow, making measurements of the diffusion coefficient difficult. Despite its noted challenges, this leaves T1 as the most promising candidate for measuring temperature changes in fatty tissue. Hynynen et al. implemented an approach using the T1-dependent magnitude signal of a fast spin echo sequence [49]. HIFU heating experiments were carried out in fatty tissue of rabbits and the sequence parameters were optimized to maximize the contrast-to- noise ratio for changes in T1. The authors demonstrated that the signal intensity changed with temperature changes at a rate of 0.97 %/°C and that magnitude signal changes were clearly visible during sonications. The DESPOT1 technique was implemented for simultaneous measurement of T1 and phase changes. As described in section 3.1, the T1 information can be extracted after two or more repetitions of a gradient echo sequence run at different flip angles, and because a gradient echo sequence is used, phase information can be gathered after every scan. To the test the method, a vial of canola oil and a vial of doped water (0.03 g/l000 ml magnesium sulfate), heated outside of the scanner, were imaged as they cooled inside the scanner. Three different flip angles were used, allowing the PRF temperature maps to be acquired every 4.5 s and T1 data to be collected every 13.5 s. Room temperature reference vials of oil and doped water were included to correct for field drift and fiberoptic temperature probes were used to monitor the actual temperatures of all vials. Results are shown in Figure 3.8. Both data sets show 58 Figure 3.8. Phase and T1 information acquired simultaneously during a phantom cooling experiment using the DESPOT1 sequence. Both metrics show good linearity. good linearity, indicating that phase and T1 temperature information can be acquired simultaneously using this method. This would allow the PRF technique to be applied for temperature measurements in voxels containing aqueous tissue, while the T1 information could be used to track temperature changes in voxels containing fat. The problem of nontemperature related field changes as it applies to the PRF method has been successfully approached in two different ways. The first method applies only to periodic motion. In this technique, many images are taken before the HIFU treatment begins in order to capture the tissue of interest during all different stages of the motion cycle [61-62]. These images are stored to be used as a look-up library of reference images. When the HIFU treatment begins, the current image is acquired and the most closely corresponding image from the look up library is used as the reference for the phase difference calculation. Because the two images will be from the same part of the motion cycle, the voxels should be closely registered and the bulk susceptibility of the 59 two images should also almost equivalent, leaving only phase changes that are due to temperature changes. This approach does not correct for inadvertent subject motion during the treatment because this type of motion cannot planned for and imaged in advance. The second method is known as reference-less PRF [63-64]. This technique does not rely on a previously acquired image to use as a reference in the phase subtraction calculation, but rather obtains this information from nonheated regions within the current time frame. The various steps are shown in Figure 3.9. The phase of the current image (Figure 3.9A) is unwrapped (Figure 3.9B) and a region surrounding the HIFU hotspot, where no heating has occurred, is identified (Figure 3.9C). A 2-D, Nth-order polynomial is fit to this region (Figure 3.9D). The extension of this polynomial into the heated region is used as the estimate for the background phase and subtracted from the actual phase of the current image, with the phase difference (Figure 3.9E) converted into a temperature Figure 3.9. Steps of the reference-less phase difference method. Instead of subtracting the phase from an earlier acquired image, a polynomial is fit to the background phase from a nonheated region of the current image. This fitted polynomial is used as the reference phase for subtraction. 60 difference. The method can correct for bulk susceptibility changes that are periodic or unexpected, provided that there is a non-heated region of aqueous tissue surrounding the temperature rise with sufficient SNR. In addition to nontemperature related phase changes, motion causes many other problems for MR thermometry. Intrascan motion creates image artifacts that typically present as ghosting or blurring in the phase-encode direction. The best way to avoid intrascan motion artifacts is to image with a scan time that is short compared to the time scale of the motion. Respiratory motion, for example, causes displacements of a few centimeters over a period of 2 - 6 seconds. It has been found that imaging durations of 400 ms or less are fast enough to adequately "freeze" this motion and avoid the unwanted image artifacts [65]. Even if images can be acquired fast enough to freeze the motion during scanning, registration problems remain due to interscan motion. The motion needs to be detected, estimated, and compensated for. One way of detecting motion is to acquire navigator echoes, where only a few central lines of k-space are acquired and projected into image space to obtain a profile of the object [66]. This method can detect bulk, rigid body motion, but cannot measure more complex motion or deformations. Another detection method relies on tracking easily identifiable landmarks within the imaging field of view. For example, larger vessels within the liver can be used. Motion estimation and compensation is a robust field that extends far beyond MRI applications. For the purposes of MR thermometry, and MRgHIFU in particular, the technique used should be compatible with real time applications. This requires an algorithm that does not depend 61 on user inputs or intervention and that is computationally efficient. For a good review of the methods, see the review paper from de Senneville et al. [65]. The final area where major advancements in MR thermometry have been made is in accelerating the image acquisition time. In general, two different approaches have been taken: modifying the image pulse sequence to more efficiently acquire the k-space data, and creating reconstruction algorithms that produce artifact free images from undersampled k-space data. The time savings gained from more efficient sampling or from sampling less data can be used to improve any combination of temporal resolution, spatial resolution, volume coverage, and SNR. Two noteworthy changes to the standard gradient echo sequence have been presented for improved temporal resolution in MR temperature imaging. The first change involves shifting the echoes such that spins excited by the nth RF pulse are read out during the (n+1)th acquisition period [67-68]. This is accomplished by applying an additional gradient of strength M right after the RF pulse and another gradient of strength -2M right before the next RF pulse. The purpose is to create a situation where TE > TR, so that long TE values can be used without lengthening TR. This allows the desirability of a TE that is closer to the optimal value of T2 * without sacrificing temporal resolution due to a long TR. The second change is to add an echo-planar readout to the sequence which acquires two or more k-space lines after every RF pulse. Both multishot and single-shot versions of echo-planar imaging (EPI) have been implemented [69-70]. The single shot version can acquire data as fast as 200 ms per slice; however, distortions at tissue interfaces can be a problem due to the long echo train length. 62 Many different techniques have been presented for improving temporal resolution by undersampling the k-space data and then using various schemes to reconstruct the images without aliasing artifacts. Parallel imaging methods take advantage of the different sensitivities of multiple receiver coils. Using m coils and SENSE or SMASH reconstructions [71-72], one can achieve a reduction factor of R (R ≤ m), but at the cost of having the SNR reduced by at least a factor of R . The reduction factors achieved by these parallel imaging techniques are also limited by the number of coils in use. Other methods to reconstruct images from undersampled measurements use an implicit or explicit a priori model in the time direction. The most basic of these is the sliding window reconstruction, which assumes that the signal is not varying rapidly in time. In the sliding window reconstruction, a subset of the k-space data is acquired for the current time frame, and the missing lines are filled in with the k-space data from the previous time frame. As the process repeats, the full k-space for time frame (n) will be made up of the most recently acquired data, plus data from time frames (n-1), (n-2), (n-3), etc. Other undersampled reconstruction methods include UNFOLD, k-t BLAST, keyhole, and RIGR [73-76]. The reduction factors realized by these techniques range from 2 to 8. While these techniques represent a promising step forward, the acceleration factors achieved are still not quite enough to obtain the desired spatial resolution, temporal resolution, and volume coverage necessary for proper monitoring and control of MRgHIFU procedures. CHAPTER 4 SPATIAL SAMPLING CONSIDERATIONS FOR ACCURATE MR TEMPERATURE MEASUREMENT This chapter is based on a paper titled, "The Effects of Spatial Sampling Choices on MR Temperature Measurements", authored by Nick Todd, Urvi Vyas, Josh de Bever, Allison Payne, and Dennis L. Parker. The paper was submitted to the journal Magnetic Resonance in Medicine in December 2009. 4.1. Introduction As discussed in Chapters 2 and 3, MR temperature measurements have been shown to play a crucial role in the planning, monitoring, control, and assessment of thermal therapies for noninvasive treatment of tumors [24, 29, 33, 77-78]. Because these treatments are noninvasive, it is paramount for clinicians to be able to predict the heat deposition pattern, monitor the treatment in real time, determine when a region of tissue has been treated to satisfaction, and retrospectively assess the amount of damage done to the treated area. MR imaging can play an important role in each of these stages due to its ability to measure temperature changes through the temperature dependence of the proton resonance frequency. 64 Investigators have used MR temperature measurements to obtain estimates for tissue thermal and acoustic properties which, in turn, can be used to assist in the patient treatment planning process [79-81]. MR temperature imaging is also widely used during thermal therapies to monitor in real time where the thermal energy is being deposited in the tissue, ensuring that the targeted area is receiving treatment and that normal tissue is not being damaged [77, 81-84]. As part of this monitoring, many investigators use either a maximum temperature or an accumulated thermal dose as the criterion for determining when a region of targeted tissue has been sufficiently necrosed. Total accumulated thermal dose can be calculated retrospectively after a treatment is complete and used as part of the histological analysis to compare predicted tissue necrosis with actual necrosis [36]. For all of these reasons, it is vitally important that MR temperature measurements be accurate, reliable, and repeatable. MR temperature maps are necessarily a discrete representation of a physical quantity that is continuously varying in both space and time. The MR sampling scheme may affect the discrete temperature values measured. Due to averaging effects, different choices for the sampling grid location, voxel size, and scan time can lead to different measurements of the same underlying temperature distribution. These differences could have a large impact on how a thermal therapy procedure is planned, carried out, and assessed. This is especially true for high intensity focused ultrasound (HIFU) treatments, where the spatially nonuniform deposition of energy into a small volume creates large temperature gradients and a focal spot with dimensions on the order of an MR voxel. In this chapter, the extent to which sampling choices affect MR temperature measurements are investigated and quantified. While temporal averaging does affect 65 temperature accuracy, this chapter only considers the effects due to variations in the spatial sampling scheme, specifically the spatial resolution and image sample grid location. Simulation results are presented quantifying the effects that variations in the sampling scheme have on measuring the maximum temperature and total accumulated thermal dose. Experimental results are presented to validate the simulation results and demonstrate that postprocessing techniques can be used to reduce temperature measurement errors by manipulating the sampling grid properties. 4.2. Methods 4.2.1. Theory Thermal dose is a method of estimating the accumulated damage done to tissue cells through heating. The accumulated thermal dose over N time frames is calculated according to [35]: Σ= = +Δ − Δ N n n D R T Tn t 1 ( 43) 43 0 [4.1] where D43 is the thermal dose that would be equivalent to continuous heat application at 43°C (Cumulative Equivalent Minutes or CEM), T0 is the baseline temperature and is set to 37°C in this work, ΔTn is the measured temperature change at time frame n (°C), Δtn is the time step between frames in minutes, and R is an empirically determined constant set to 4 for ΔTn less than 6°C and to 2 for ΔTn greater than 6°C. Although tissue necrosis thresholds vary for different tissue types, most investigators use 240 CEM as the threshold for complete tissue necrosis [36]. The nonlinear relationship between dose and 66 temperature means that thermal dose can accumulate very rapidly at high temperatures and that the calculation is extremely sensitive to errors in temperature measurement. Figure 4.1 plots thermal dose rate as a function of temperature rise above body temperature (37°C). At 16.8°C above baseline, 30 CEM of dose is delivered every second; and at 19.8°C above baseline, a lethal dose of 240 CEM is delivered every second. For every degree error in the temperature measurement, the calculated thermal dose rate will be scaled by a factor of two. For example, if an MR image measures a temperature rise of 17°C in a voxel when the actual rise was 20°C, then the thermal dose accumulated during that time frame will be calculated to be eight times less the actual value, an error that could result in a mistaken judgment of nonnecrosis. Figure 4.1. Thermal dose rate as a function of temperature change. The thermal dose rate doubles with every 1°C increase in temperature and at 20°C above body temperature, a lethal dose is delivered every second. 67 For HIFU heating applications with a small focal spot, the temperature increase will vary spatially over the volume of a voxel. The discrete temperature measured within a voxel is therefore an average of the continuous temperature over the entire voxel. In MR imaging, the signal sensitivity is not constant throughout the voxel, but rather is determined by a sinc-like voxel sensitivity function (VSF) [85]. In one dimension, the VSF is given by: ( ) ⎟ ⎟⎠ ⎞ ⎜ ⎜⎝ ⎛ Δ ⎟⎠ ⎞ ⎜⎝ ⎛ = ⋅ Δ N x x x x N q x x X π π sin sin 1 [4.2] where Nx is the number of voxels and Δx is the voxel dimension. The 1-D VSF is plotted in Figure 4.2 for the voxel of interest and its 5 neighbors. The three dimensional VSF would be a product of q(x), q(y) and q(z). The signal sensitivity is maximal at the center of the voxel volume, but drops off quickly away from the center. For a three dimensional isotropic voxel, the sensitivity is 64% of the maximum at the sides of the voxel, 41% at the edges, and only 26% at the corners. This nonuniform sensitivity has implications for how the underlying temperature distribution is measured and discretized. Figure 4.3 shows a one dimensional example of the averaging effects due to the VSF and what happens when the voxel sampling grid is shifted in space. In the top plot the 2mm voxel is centered over the peak of a temperature distribution that has a maximum temperature rise of 23.2°C. The voxel sensitivity is shown on the same plot, scaled by a factor of 10 for visibility. In MR imaging, the temperature information is contained in the phase of the signal, but the magnitude is 68 Figure 4.2. A one dimensional plot of the voxel sensitivity function (VSF). Signal sensitivity falls off rapidly away from the center of the voxel. For a 3-D voxel, the sensitivity is 64% of the maximum at the sides of the voxels 41% at the edges, and 26% in the corners. scaled by the VSF. When the phase and scaled magnitude information is combined over this one dimensional voxel the maximum temperature is measured as 22.4°C. The bottom plot of Figure 4.3 shows the same situation, but with the sampling grid shifted to the left by 0.5mm. Now the peak of the temperature distribution is away from the center of the voxel and weighted less because of the fall off in the VSF. In this case, the measured temperature comes out to be 21.5 °C, almost two degrees less than the actual peak temperature. The effect will be more pronounced in three dimensions as the peak temperature could be located in the corner of a voxel where the VSF is even less sensitive. 69 Figure 4.3. Simulated temperature distribution with VSF (scaled by a factor of 10 for clarity). A) When centered in one voxel, the maximum temperature is measured to be 22.4 °C. B) When the sampling grid is shifted, the peak temperature is weighted less by the VSF and the maximum temperature is measured to be 21.5 °C. An MR image is the discrete Fourier transform of discrete samples of the continuous spatial frequency space (k-space) of the object, which allows the image grid positions to be chosen retrospectively. In practice, manipulation of the sampling grid can be done in two different ways: linear phase shifting and zero-filled interpolation (ZFI). Each technique only affects the positioning of the image points but has no effect on the shape or size of the VSF and therefore has no effect on the actual image spatial resolution. 70 The application of a linear phase in one direction in the k-space domain will result in a shift of the sampling grid in the same direction in the image domain. For 2-D MR data, the sampling grid can be shifted in-plane, while for 3-D MR data the sampling grid can be shifted in all three directions. Arbitrarily small slopes can be chosen for the applied linear phase, resulting in subvoxel shifts of the sampling grid. The technique does not alter the size of the image matrix and therefore the only cost is a small penalty in computation time. Zero-filling the k-space data to a larger matrix size will create a denser grid of overlapping VSF's, reducing the distance between sampling points in the image domain. While ZFI does not increase the spatial resolution of the original data, the more tightly packed VSF's will reduce the amount of area with low signal sensitivity and allow for better capture of image details that are smaller than the VSF dimensions. In theory, ZFI can achieve arbitrary placement of the sampling grid, but this would require creating an arbitrarily large image matrix. 4.2.2. Simulations Ultrasound beam propagation software was used to model the power deposition of a 256-element randomly spaced phased array ultrasound transducer (Imasonics, Inc.) at 0.1 x 0.1 x 0.1 mm (0.1 mm isotropic) spatial resolution [86]. Using this power deposition matrix and the Pennes bioheat equation [87], simulated temperature maps with 0.1 mm isotropic spatial resolution and 1 second temporal resolution were created for a 15-second, 100W step input of power followed by a 5-second cooling period. All thermal and acoustic parameters used in the simulation process were assumed to be constant and are summarized in Table 4.1. The full width at half maximum (FWHM) dimensions for 71 Table 4.1. Thermal and acoustic parameters used for numerical simulation of temperature maps. Property Value k (W/m*K) 0.47 Ct & Cb (J/kg*°C) 4186 ρ (kg/m3) 1000 W 0 c (m/s) 1500 (tissue & water) α (np/cm*MHz) 0.04835 the ultrasound power deposition were 7.8 x 2.3 x 2.3 mm, while the FWHM dimensions for the temperature distribution at the hottest time frame, 10.7 x 3.3 x 3.3 mm, were slightly larger due to thermal diffusion. Complex MR image data was created from the simulated temperature maps by converting the temperature into phase [43] and assuming constant image magnitude over the region of interest. These complex images were then transformed into k-space using a 3-D FFT. An image of the central slice of the 0.1 mm isotropic temperature maps at the hottest time frame is shown in Figure 4.4. Modifications to the spatial sampling scheme of the simulated temperatures were done by manipulating the simulated k-space data in two ways. First, the actual spatial resolution of the temperature maps was changed to create 2-D and 3-D voxels with varying voxel dimensions. The 2-D voxels had f