T2 Decay Curve Acquisition and Analysis in MRI Noise Considerations, Short T2, and Bx Field Encoding by Craig K. Jones M.Sc, The University of Western Ontario, 1997 B.Sc, Simon Fraser University, 1992 A THESIS SUBMITTED IN PARTIAL FULFILMENT OF THE REQUIREMENTS FOR THE DEGREE OF DOCTOR OF PHILOSOPHY in The Faculty of Graduate Studies (Department of Physics and Astronomy) We accept this thesis as conforming to the required standard THE UNIVERSITY OF BRITISH COLUMBIA July 31, 2003 © Craig K. Jones, 2003 In presenting this thesis in partial fulfilment of the requirements for an advanced degree at the University of British Columbia, I agree that the Library shall make it, freely available for reference and study. I further agree that permission for extensive copying of this thesis for scholarly purposes may be granted by the head of my department or by his or her representatives. It is understood that copying or publication of this thesis for financial gain shall not be allowed without my written permission. Department of Physics and Astronomy The University Of British Columbia Vancouver, Canada Abstract The myelin sheath is a lipid bilayer wrapped around axons in the human brain. The myelin allows faster conduction and uses less energy than along non-myelinated axons. Diseases, such as multiple sclerosis, break down the myelin sheath resulting in cognitive and physical disability. Water trapped between the myelin bilayers has a T2 ~ 15 ms compared to intra/extra-cellular water with a T2 ~ 80 ms and cerebrospinal fluid which has a T2 ~ 2000 ms. A multi-echo MRI pulse sequence is used to acquire a T2 decay curve which is fit using a non-negative least-squares (NNLS) curve fitting algorithm to calculate the signal as a function of the T2 (T2 distribution). The myelin water fraction (MWF) is calculated as the signal with a T2 < 50 ms divided by the total signal in the T2 distribution. Quantification of the MWF in vivo, is slow, prone to errors due to artifacts in the T2 decay curve, and requires many data acquisition averages to attain the necessary high signal-to-noise ratio. Each of these areas were addressed in this thesis. First, I compared the accuracy and precision of MWF estimates from a set of simulated and in vivo multi-echo data with and without noise reduction filters applied. The MWF estimated from anisotropically filtered multi-echo data had the smallest variability, over homogeneous regions, compared to MWF estimates from other filtered data. Second, I created a new technique to optimize coefficients that when linearly combined with multi-echo data, result in fast, accurate estimates of the MWF. Simulations and in vivo estimates of the MWF were as accurate as those from the NNLS algorithm and 20,000 times faster. Finally, the standard multi-echo sequence was optimized to remove artifacts due to inaccurate refocusing pulses. I created a novel T2 curve fitting algorithm, based on NNLS, to accurately estimate the T2 distribution and refocusing pulse flip angle from a T2 decay curve. Simulations and data acquired on phantoms showed the new technique was as accurate as the NNLS algorithm in quantifying the T2 distribution parameters. In vivo myelin water maps were as good as those estimated from the optimized multi-echo pulse sequence. ii Table of Contents Abstract ii Table of Contents iiList of Tables viList of Figures ix List of Algorithms xi Acknowledgements xiGlossary xiv 1 Introduction 1 1.1 Magnetic Resonance Imaging 1 1.2 T2 and Ti Relaxation !' 2 1.2.1 Spectral Density 2 1.2.2 Bloch Equations 3 1.3 Central Nervous System 4 1.4; Tissue Relaxation 6-1.5 T2 Decay Curve 7 ' 1.5.1 Mono-Exponential T2 Relaxation 7 .1.5.2 Multi-Exponential T2 Relaxation . . 8 1.6 T2 Decay Curve Acquisition 9 1.6.1 K-Space . . .1.6.2 Optimized Multi-Echo Pulse Sequence 10 1.6.3 Standard Multi-Echo Pulse Sequence 1 1.6.4 Noise in Decay Curve 14 1.6.5 Artifacts Influencing the T2 Decay Curve 16 1.7 T2 Decay Curve Analysis 8 1.7.1 Decay Curve Inversion (Decay Curve to T2 Distribution) 11.7.2 T2 Decay Curve Requirements 24 1.8 Phantoms 25 1.9 Overview of Thesis . . 26 in Table of Contents 2 Myelin Water Fraction Noise Reduction 27 2.1 Introduction 22.2 SNR Effect on Myelin Water Fraction 8 2.3 Description of Spatial Filters 29 2.3.1 Anisotropic Diffusion Filter 30 2.3.2 Susan Filter 31 2.3.3 Median Filter2.3.4 Wavelet Filter 2 2.4 Spatial Filters on Simulated Data . . . 32 2.4.1 Methods2.4.2 Results 36 2.4.3 Conclusions 40 2.5 Spatial Filters Applied to Zn Vivo Data , .40 2.5.1 Data Acquisition and Analysis 42.5.2 Results 1 2.5.3 Conclusions 44v 2.6 Spatial Filtering vs Averaging .44 . 2.6.1 Methods v.. • • • ••'44. 2.6.2 Results 46 2.6.3 Discussion :. . . . 50 2.6.4 Conclusions 53 2.7 Conclusions3 Myelin Water Fraction from Linear Combination 55 3.1 Introduction3.1.1 Linear Combination for T2 Selectivity 6 3.1.2 Properties of a Linear Combination . . 57 3.1.3 Overview 53.2 Calculation of Coefficients 8 3.2.1 Coefficients for Short T2 59 3.2.2 Coefficients for All T2 60 3.2.3 Description of the Algorithm '. . . . 61 3.2.4 Results 62 3.2.5 Discussion 3 3.3 Linear Combination Simulations 64 3.3.1 Model3.3.2 Results and Discussion3.3.3 Conclusions 65 3.4 Selectivity of Mono-Exponential Phantoms 66 3.4.1 Methods3.4.2 Results and Discussion 63.4.3 Conclusions 8 3.5 In Vivo Myelin Water Fraction: NNLS vs. T2 Filter 63.5.1 Methods 69 3.5.2 Scan Results3.5.3 Conclusions 72 iv Table of Contents 3.6 Calculation of Coefficients by Matrix Inversion 72 3.6.1 Calculation of Coefficients 73.6.2 Results and Discussion 4 3.6.3 Conclusions 76 3.7 Calculation of Coefficients by Backus-Gilbert 73.7.1 Calculation of Coefficients : 76 3.7.2 Results and Discussion 77 3.7.3 Conclusions 80 3.8 Optimal Echo Times .3.8.1 Filter Measures3.8.2 Equally Spaced Echoes 82 3.8.3 Non-Equidistant . . .: 85 3.8.4 Discussion. . . 8 3.8.5 Four Echo Non-equidistant Sampling . 90 3.9 Conclusions • • • • .-92 4 Non-180 Refocusing Pulses 94 4.1 Introduction . .4.1.1 Motivation4.1.2 Previous Work 5 4.1.3 Implementation of Algorithm 98 4.1.4 Overview of Work 102 4.2 Mono-Exponential T2 Decay Curve Fit 104.2.1 Mono-Exponential Fit Algorithm4.2.2 Monoexponential Fit - Noise Simulation 103 4.2.3 Phantom Verification 105 4.3 Multi-Exponential T2 Decay Curve Fit : ,111 4.3.1 Multi-Exponential Fit Algorithm . . . .'Til 4.3.2 Bi-Exponential Parameters From Incorrect a 115 4.3.3 NNLS° Multi-Exponential Simulations 116 4.3.4 Phantoms 117 4.3.5 In Vivo Multi-Exponential Fits . . • 120 . 4.4 Multi-Slice, Multi-Echo : 126 4.4.1 Introduction 124.4.2 Methods4.4.3 Results and Discussion '. . . . 127 4.4.4 Conclusions 135 4.5 Non 90° Excitation Pulse - - 136 4.5.1 Decay Curve Calculation 134.5.2 Simulation Methods : 136 4.5.3 Simulation 1: Inaccurate Excitation and Accurate Refocusing, Assume both Accurate 137 4.5.4 Simulation 2: Inaccurate Excitation and Inaccurate Refocusing, Assume Accurate Excitation 138 4.5.5 Simulation 3: Inaccurate Excitation and Inaccurate Refocusing 139 4.5.6 Discussion 140 v Table of Contents 4.5.7 Conclusions 141 4.6 Non-Constant Train of Flip Angles 144.6.1 Simulations 3 4.6.2 Scanning 144 4.6.3 Conclusion 9 4.7 Conclusions5 Conclusions 150 5.1 Filtering5.2 Linear Combination 155.3 B!/T2 2 5.4 Preferred T2 Quantification Protocol 153 Bibliography ......... 154 A Related Publications • 162 List of Tables 1.1 Ti and T2 of mono-exponential nickel/agarose phantoms 25 2.1 The mean global error in myelin water fraction over the 10 noisy data sets for each filter 37 2.2 The mean error in myelin water fraction over the 10 noisy data sets for each filter for the area of 15% myelin water 32.3 The mean error in myelin water fraction over the 10 noisy data sets for each filter for the area of 30% myelin water ; 38 2.4 In vivo myelin water fraction as a function of filter ,41 2.5 The number of times each filter had the smallest myelin water fraction ....... . 42 2.6 Number of holes per ROI as a function of iterations of the anisotropic diffusion filter 48 3.1 The coefficients for the short T2 and all T2 linear combination method 62 3.2 Mean and standard deviation of the signal in the short T2 phantom image ...... 67 3.3 Coefficients for the short, medium and long T2 selection based on matrix inversion. 74 3.4 Coefficients calculated by the Backus-Gilbert method 73.5 Coefficients for short and all T2 filter based on four unevenly spaced echoes. .... 90 4.1 Mono-exponential noise simulation as a function of a v 104 4.2 Mono-exponential noise simulation as a function of a and SNR scaled by a function of the refocusing pulse . . . 104.3 Parameters of phantom derived from mono-exponential fit ,. 106 4.4 Bi-exponential noise simulation to assess accuracy and precision of parameter esti-.... mation from non-1800 refocusing pulse train ;. 116 4.5 Noise simulation to estimate variability of multi-exponential decay curve parameters. 117 4.6 T2 time calculated from the MESEQ sequence and the MESEC (with the flip angle calculated) .119 4.7 Refocusing pulse flip angle calculated from multi-echo decay curve •'. 119 4.8 Estimated Ti from curve fits of five points (TR = 100,300,800,2000,3000 ms) from a pulse saturation experiment .114.9 Percent myelin water estimated from scans with (xnom = 150° and anom = 180°. . 121 4.10 Refocusing pulse flip angle calculated over the regions from the scans with <xnom — 150° and anom = 180° .123 4.11 Geometric mean T2 of the medium component from select regions throughout the image 124 4.12 Decay curve parameters estimated from in vivo scan as a function of pulse sequence. 133 4.13 The excitation pulse was varied from 90° to 50° and the parameters were calculated assuming only the refocusing pulse was incorrect 138 vii List of Tables 4.14 The excitation and refocusing pulses were varied with ae = 90°, 85°,... ,50° and ar = 2ae 139 4.15 The excitation and refocusing pulses were varied with ae = 90°, 85°,..., 50° and ar = 2ae 140 4.16 Simulation to assess decay curve parameter estimation from non-constant refocus ing pulse train (16 echoes) 144 4.17 The T2 and refocusing pulse flip angle (for the first echo) calculated from each train of refocusing pulses 7 4.18 The resulting x2 misfit given the scaling of the RF pulses was 1.0,1.0,1.0, S90 from s<4) 148 viii List of Figures 1.1 Myelinated axon '•• 5 1.2 Electron-micrograph of a myelinated cell 6 1.3 Simulated mono-exponential decay curves 8 1.4 Standard MRI image and corresponding k-space ' 10 1.5 Optimized single-slice, multi-echo pulse sequence (MESE0) 11 1.6 Conventional multiple echo 2D spin-echo pulse sequence (MESEC) .......... 12 1.7 Conventional multiple echo 3D spin-echo pulse sequence (MESEC) . . . ... ... . 13 .: 1.8 Rice distribution as a function of the mean . .14 1.9 Decay and T2 distribution of a bi-exponential white matter model . . . 19 1.10 Least squares and regularized T2 distributions of a white matter model ...... 21 2.1 Histograms of myelin water fraction from simulated decay curves 29 2.2 Simulated myelin water fraction image annotated with six regions 33 2.3 Simulated myelin water fraction image (truth) from which the 32-echo data set was constructed 35 2.4 Myelin water maps calculated from spatially filtered, simulated multi-echo data. . 39 2.5 Myelin water maps calculated from spatially filtered MRI data 43 2.6 Myelin water maps from unfiltered and filtered (anisotropic diffusion filter) multi-echo data .42.7 Histograms of unfiltered and filtered myelin water fractions . 51 . 3.1 Constraints on the short T2 filter. 60 3.2 Constraints on the "all" T2 filter 1 3.3 T2 filters for short T2 components and all T2 components 63 3.4 Myelin water fraction noise simulation to compare linear combination and NNLS. . 65 3.5 Linear combination short T2 selectivity of phantoms 7 3.6 Mean signal from phantoms in selected image versus T2 filter 68 3.7 In vivo myelin water maps calculated by linear combination and NNLS 70 3.8 Myelin water fractions calculated per region from the linear combination method and NNLS 71 3.9 Bland-Altman plots for white matter regions and gray matter regions 73.10 T2 filters based on coefficients calculated by matrix inversion 73 3.11 Brain images T2 filtered from coefficients calculated by matrix inversion 75 3.12 Coefficients calculated by Backus-Gilbert method 78 3.13 T2 filters based on coefficients calculated by Backus-Gilbert method 79 3.14 The measures of the effectiveness of a short T2 filter 81 3.15 The coefficients as a function of the number of echoes and TE 83 3.16 Width of T2 filter for equally spaced echo times 84 ix List of Figures 3.17 The coefficients as a function of the number of echoes and TE 86 3.18 Width of T2 filter for unequally spaced echo times 87 3.19 The T2 filter, exp (-TE/T2) as a function of TE 9 3.20 The short T2 filter as a function of a consecutive subset of echoes based on the thirty-two coefficients 83.21 Short T2 filter calculated from four echoes compared to filter calculated from 32 echoes 91 3.22 Myelin water map calculated by linear combination from four echoes and from 32 echoes4.1 Phase coherence diagram 96 4.2 Decay curves, collected with ap — 180°, fit to obtain T2 and a 107 4.3 Decay curves, collected with ap = 150°, fit to obtain T2 and a 108 4.4 Decay curves, collected with ap = 110°, fit to obtain T2 and a 109 4.5 Decay curves, collected with ap = 90°, fit to obtain T2 and a 110 4.6 The x2 misfit plotted as a function of refocusing pulse flip angle ' 112 4.7 The rank and inverse of the condition number of AQ based on a 32 echo decay . . 114 4.8 The rank and inverse of the condition number of Aa based on a 200 echo decay . . 115 4.9 The myelin water maps calculated from the pulse sequence with a refocusing pulse flip angle prescribed as 180° and 150° , 121 4.10 Refocusing pulse flip angle maps calculated with a refocusing pulse flip angle pre scribed as 180° and 150° 122 4.11 The arithmetic mean T2 maps calculated from the pulse sequence with a refocusing pulse flip angle prescribed as 180° and 150° . . 123 4.12 The x2 maps calculated from the pulse sequence with a refocusing pulse flip angle prescribed as 180° and 150° 125 4.13 The myelin water maps calculated from the pulse sequence with a refocusing pulse flip angle prescribed as 180° and 150° .126 4.14 The TE = 10 ms image from the four average, 32 echo MESEG sequence and six slices of the six slice slab MESEC acquisition 129 4.15, Myelin water maps calculated from a 6-slice slab MESEC . . . .... 130 4.16 Difference between the myelin water map calculated from the MESEC sequence and the myelin water map calculated from the MESEC sequence 131 4.17 In vivo decay curve from MESEC and MESE0 from the genu ofthe corpus callosum (WM) 134 4.18 Maps ofthe refocusing pulse flip angle calculated from a 6-slice slab MESEC .... 135 4.19 Spin-echo sequence with three 180° followed by 13 90° . . . 142 4.20 Spin-echo sequence with an exponential decay of the refocusing pulse amplitudes. . 143 4.21 Decay curve fit of decay curves collected with non-constant refocusing pulse train, phantom with T2 = 85 ms 145 4.22 Decay curve fit of decay curves collected with non-constant refocusing pulse train, phantom with T2 = 23 ms 6 4.23 Measured and fitted data from non-constant train of refocusing pulses 148 x List of Algorithms 1.1 Small norm NNLS regularization 23 4.1 Signal calculation from a train of refocusing pulses . 100 4.2. Ti relaxation 101 4.3 T2 relaxation . . .: . . 101 4.4 Phase evolution of magnetization . 10xi Acknowledgements When I first moved to Vancouver in 1997, I began working for Dr. Don Paty, Dr. David Li, and Andrew Riddehough at the MS/MRI Research Group. After a little over a year, it became clear that to pursue further research I would have to do a PhD. My first conversations with Dr. Paty and Dr. Li, were overwhelmingly supportive. Dr. Paty's comment was, "I was wondering when you would come and talk to me about this." With a lot of work from Andrew, an arrangement was worked out such that I worked in the MS/MRI Research group very part-time while I did my PhD. The financial, academic and other support (Ken Bigelow and Vilmos Soti, in particular), during the past four years, has allowed me to complete this PhD. Without the MS/MRI Research Group, this thesis would not exist. Thanks! Anne, my wife, was the one who pushed me, in the beginning, to do a PhD. I will never forget the evening we sat at Starbucks, talking over the direction of my career. It became more evident the only way to pursue what I wanted was to do a PhD and the only thing holding me back was the desire to support my family. Anne challenged me on that desire and made it clear that we would make it through just fine. In hindsight, she was correct. We have been able tc live our lives, for the most part as if I had a (real paying) job, and, in fact, have probably had a better lifestyle. The first year of my PhD was spent with my supervisor, San. I was doing courses and trying to find a project that to which I could take hold. Then, after some discussion with Alex, and some choices by other grad students, some interesting possibilities came along that led to the work in this thesis. It is interesting to see that the direction of the thesis, in the first year, is very different than the content of the work herein. I would like to thank Alex for the support, interesting projects and allowing me to take on the project. I would like to thank San for the many discussions we had on the different topics. The meetings with San were always extremely interesting, his enthusiasm for research and science always rubbed off on me. His knowledge of xn Acknowledgements physics and MRI as well as his natural ability to describe complex ideas, using analogies, greatly enhanced my understanding. On more than one occasion my patient wife was waiting for me as I got lost in conversation with San ("talking to San again?") Other people with influence in the thesis were my committee members: Dr. Anna Celler, Dr. Carl Michal, and Dr. Ken Poskitt. Each of them brought knowledge from their respective areas, which enhanced the quality of the science in this work. My committee meetings were very positive times and put the work in perspective. . There are many friends who have been a great influence and who have enabled such an interesting time for Anne and me in Vancouver. For myself, in particular, Ken Nixon, who I met with Wednesday mornings for several years. The one question Ken always asked me, "Do you have a chapter done?". He is a great friend and someone who helped put the work in perspective. .... Finally, the one person who has had a significant impact in my life is Dr. Ken Whittal. Ken has been a great friend, and more importantly, a mentor for me. He was very patient, with my. many questions, and always ready and available. The word mentor is defined as "...[a] wise and trusted counselor and teacher..." and Ken was... xiii Glossary Abbreviation Meaning BG Backus-Gilbert CSF Cerebrospinal Fluid FOV Field of View (cm) FWHM Full-Width-at-Half-Max GM Gray Matter 'ME" " • - •• Multi-echo MESE Multi-echo Spin-Echo Pulse Sequence MRI Magnetic Resonance Imaging MS Multiple Sclerosis MT Magnetization Transfer MWF Myelin Water Fraction NNLS Non-Negative Least Squares RF Radiofrequency ROI Region-of-Interest RMS Root-Mean-Square SAR Specific Absorption Rate SI Signal Intensity SNR Signal to Noise Ratio T Echo spacing (ms) TE Echo Time (ms) Ti Spin-Lattice Relaxation Time (ms) TR Repetition Time (ms) T2 Spin-Spin Relaxation Time (ms) xiv Glossary Glossary WM White Matter xv Chapter 1 Introduction In magnetic resonance imaging, information about the local water environment is quantified by Ti and T2 proton relaxation times and the density of the protons, p. The measurement process and analysis is complex, but enables different information to be obtained from the local'water environment. The thesis will describe new techniques for better quantification and analysis of T2, as well as the encoding of extra information using "artifacts" that arise from data collection. The introduction has an overview of magnetic resonance imaging (Section 1.1) and the physics of relaxation (Section 1.2). Then, the central nervous system is descrbied (Section 1.3) and the motivation to use quantitative T2 to assess tissue and pathology (Section 1.4). The T2 decay curve is discussed (Section 1.5) and then the data collection (Section 1.6) and data analysis (Section 1.7). A set of phantoms was created to validate the acquisition and curve fitting of the decay curves (Section 1.8). 1.1 Magnetic Resonance Imaging Magnetic resonance imaging (MRI) uses a strong main magnetic field, Bo, to align the protons in the human body parallel and anti-parallel to the main field. A set of smaller magnetic field gradients, along X, Y, and Z, are used to spatially encode the protons and a radiofrequency (RF) coil is used excite the sample and to receive RF signal from the protons. The spatially encoded RF signal is reconstructed into an image. The signal intensity of each voxel in the image is a function of the number of protons in the voxel (p), the spin-spin relaxation time (T2) and the spin-lattice relaxation time (Ti). 1 CHAPTER 1. INTRODUCTION 1.2. T2 AND Ti RELAXATION 1.2 T2 and Ti Relaxation The T2 relaxation time calculated from a multi-echo decay curve acquired using the MRI scanner is related to the local water environment, the acquisition method, and the main magnetic field (Bo). The local water environment affects the interaction of the water molecules (protons). The effect of the interaction is modeled as a spectral density function and is used to define the Ti and T2 relaxation times. The relaxation times relate, classically, to the magnetization through a set of phenomenological equations derived by Bloch [1]. The local water environment in human tissue, on the scale of an MRI voxel, typically contains several distinct water compartments. Each of the water compartments may have a different T2, therefore, a modified set of magnetization' equations must be analyzed to calculate the T2 of each compartment. 1.2.1 Spectral Density The Ti and T2 relaxation times are related to the Larmor frequency of the proton, wo, and the correlation time, Tc, ofthe spins through the spectral density function [2]: Jin) M = ^iry~T^ (1.1) where is the distance between spins, A;(0) = 24/15, k(l) = 4/15, and k(2) = 16/15. The correlation time is a characteristic time scale of molecular motion. The relaxation rates, inverse: of the relaxation times, are then defined as [2]: 1/T, = k\jW(cj0) + J{2)(2uJo)\ (1.2) 1/T2 = k \jM(0) + 5-jU("o) + \jW(2u0) (1.3) where (w) is defined in Equation 1.1 and k = (^f 74/i2§J (I + 1) [2]. 2 CHAPTER 1. INTRODUCTION 1.2. T2 AND Tj RELAXATION 1.2.2 Bloch Equations The magnetization, as a function of time, was defined in a set of phenomenological equations by Felix Bloch [1]: c? . . r,,/N -i M,.(/.) „ MM) „ Mo MM) „ ,- x -M(t)=7[M(t)xB]-^-^iy+ Tl * (L4) where M(i) is the magnetization vector; B is the applied magnetic field; T2 is the spin-spin relaxation time; Ti is the spin-lattice relaxation time; 7 = 42.577 MHz/T is the gyromagnetic ratio of a proton; Mo is the magnetization vector at thermal equilibrium; and x, y, and z are unit vectors along X, Y, and Z. Equation 1.4 can be written as a function of the magnetization in each orthogonal direction: jtMM) = 7M,50-^ (1.5) jMy{t) = lMxB0-^ (1.6) ±MM (,7The solution to the coupled differential equations in Equations 1.5,1.6,1.7 is: MM) =[Mx(0)cos(w0<) + Afy(0)sin(woO]exR(-t/T2). . (1.8). My(t) = [-Mx(0)sia(uQt) + My{Q)cos(u0t)}exp(-t/T2) (1.9) MM) = M0 + [Mz(0)-M0]exp(-i/T1). (1.10where OJQ = jBo and Mo is the equilibrium magnetization. By defining the transverse magnetiza tion, Mxy(t), as Mxy(t) = Mx(t) + iMy(t), special solutions to Equations 1.5, 1.6, 1.7, where the initial magnetization is aligned in the transverse plane (e.g., after a 90° excitation pulse) is given by: Mxy{t) = M0 e^2 e1Wot (1.11) Mz(t) = M0(l-e-'/Tl). (1.12) CHAPTER 1. INTRODUCTION 1.3. CENTRAL NERVOUS SYSTEM where Mxy(0) = 1 and M2(0) = 0. The frequency of precession OJQ is typically demodulated out of the transverse magnetization to obtain: Mxy(t) = M0e-*/T2 (1.13) Mz(t) = M0 (l -e-*/Ti) . (1.14) The measured signal intensity, y(ti) is a function of the magnitude of the transverse magne tization y(t) oc K\Mxy(t)\ = K Mn e-*/T2 The measured signal from the RF coil is digitized and filtered to obtain the signal intensity. The primary interest, for this thesis, is T2 relaxation, therefore, the acquisition and analysis of T2 decay curves will be discussed further. 1.3 Central Nervous System Myelin [3, 4] is a fatty sheath, composed of a bi-layer of lipids, which surround axons in the human brain. The biochemical composition is approximately 70% lipid and 30% protein (myelin basic protein, proteolipid proteins) [5]. In the central nervous system (brain, cerebellum, and spinal cord) the myelin sheath are processes from oligodendroglial cells, arid in the peripheral nervous system, from Schwann cells. Each oligodendrocyte can myelinate as many as 50 different axons [6]. Each oligodendrocyte process wraps around the axon multiple times. The sheath is interrupted, along the axon, by the Nodes of Ranvier, approximately 1 mm apart [6], and, therefore, the electrical conduction jumps from one Node of Ranvier to the next at approximately TOO m/s. This jumping (saltatory conduction) results in faster conductivity, relative to an unmyelinated axon, and requires less energy. To maintain the conduction velocity of a myelinated axon 20 microns in diameter, an unmyelinated axon would have to be several millimeters in diameter [3]. The thickness of the each myelin bilayer is approximately 0.1 /um [7] and each radial line in the myelin sheath is approximately 300 - 400 A apart [7]. There are many diseases (e.g., multiple sclerosis, adrenoleukodystrophy [9, 10], and attention deficit-hyperactivity disorder [11]) which affect the myelin, and therefore, the water trapped be tween the myelin bilayers. To assess the disease progression or drug effect, it is important to analyze the progression of demyelination in vivo. MRI is well suited for this and therefore, will continue to be an important factor in assessment, research and treatment of disease. 4 CHAPTER 1. INTRODUCTION 1.3. CENTRAL NERVOUS SYSTEM Figure 1.1: Myelinated axon [8] showing the multiple layers of the myelin sheath around the nerve and the nodes of Ranvier. The integrity of the myelin bilayer can be measured by measuring the relaxation time, T2, which is a function of the local water environment. The T2 decay curve acquisition, from a stan dard MRI pulse sequence, results in signal which decays as an exponential function of T2. Simple systems decay as a function of a single T2 (mono-exponential decay), as there is a single water compartment. Complex tissues (e.g., brain and muscle) do not have a single water compartment at the resolution of MRI (ss 1 mm in-plane resolution) and therefore, the signal decays as a sum of single T2 exponentials (multi-exponential decay). 5 CHAPTER 1. INTRODUCTION 1.4. TISSUE RELAXATION Figure 1.2: Electron-micrograph of a myelinated cell 1.4 Tissue Relaxation In 1978, Vasilescu et al. [12] published quantitative results of multi-exponential T2 decay curve analysis of frog sciatic nerve. They found three distinct T2 water compartments based on a "peeling-off" mathematical procedure applied to multi-echo data. A slowly relaxing compartment, T2 ~ 300 — 500 ms, attributed to water in the intra-cellular space. An intermediate relaxing component, T2 ~ 80 ms ascribed to axoplasmic water. And a fast relaxing component, T2 ~ 20 ms, ascribed to water closely associated with proteins and phospholipids. This paper was the first one to describe multi-exponential behavior in myelinated nerve. A little later, in 1986, Kroeker et al. [13] suggested the problem of analyzing biological re laxation times to allow the data determine the appropriate model. They appropriately assumed the data can be described by an integral of weighted exponentials and used the CONTIN curve fitting program (written by Provencher) to solve integral equations for multi-exponential data. It wasn't until 1991, when Menon et al. [14] analyzed myelinated nerves from cat brains, the short T2 peak (T2 ~ 15 ms) was attributed to water trapped between the myelin bilayers. Menon fit T2 decay curves, using NNLS, acquired from myelinated cat brain nerves and found CHAPTER 1. INTRODUCTION 1.5. T2 DECAY CURVE four distinct components: component near T2 = 1 ms attributed to phospholipid protons in the tissue, component near T2 = 12.7 ms attributed to water in the myelin layers, component near T2 = 89 ms attributed to intra and extra-cellular water, and a component not attributed to anything at T2 ~ 340 ms. The short component, T2 = 12.7 ms, contributed 6.8% to the total water in the T2 distribution. Menon also noted the myelin bound water would be in slow exchange with the intra- and extra-cellular water due to myelin lipid bilayers acting as a diffusion barrier. The first in vivo corroboration of a short T2 component was by MacKay et al. [15, 16]. In vivo T2 relaxation measurements [17, 18, 19] distinguish three water compartments in the brain: water trapped between the myelin bilayers (10 < T2 < 50 ms), intra/extra axonal water (T2 ss 80 ms) and water associated with cerebrospinal fluid (CSF) (T2 > 2000 ms). The fraction of the water trapped between the myelin bilayers, relative to the total, was approximately 1:6 based on x-ray diffraction measurements [7] among other techniques [17]. Since these first few studies, there have been many other studies to further describe and assess the short T2 component. Beaulieu et al. [20] compared T2 distributions from myelinated and non-myelinated nerves from the garfish and found the short T2 component only in the myelinated axon. More recently, many studies have compared the short T2 component to other MRI phe nomenon such as magnetization transfer imaging [21, 22, 23] and diffusion imaging [24]. Each water compartment is "visible" to MRI [17, 18, 19, 25], therefore, MRI is a unique tool to probe the water environments and assess normal and pathological tissue in vivo. ' 1.5 T2 Decay Curve 1.5.1 Mono-Exponential T2 Relaxation The T2 decay curve is the signal intensity measured from the radiofrequency coils as a function of time t, for this thesis the signal intensity will be denoted y{t). The signal intensity, y(t), is defined as the magnitude ofthe transverse magnetization y(t) = \Mxy(t)\ as defined in Equation 1.13. A vector of signal intensities will be denoted y and a vector of echo times denoted t. For standard multi-echo pulse sequences, t is the time ofthe spin echo and is typically measured in milliseconds. Figure 1.9 shows three simulated decay curves of T2 = 20 ms (circles), T2 = 80 ms (crosses), and T2 = 2000 ms (diamonds). For each curve, y = pexp (—t/T2), the proton density p = 1000 and the time of the echoes was t = 10,20,... , 320 ms. 7 CHAPTER 1. INTRODUCTION 1.5. T2 DECAY CURVE 1000r 150 TE (ms) Figure 1.3: Simulated decay curve with signal intensity S(t) = pexp (—£/T2) with p = 1000, t = 10,20,... , 320 ms and T2 = 20 ms (circles), T2 = 80 ms (crosses), and T2 = 2000 ms (diamonds). TR was considered infinite. 1.5.2 Multi-Exponential T2 Relaxation Many tissues in the human body have more than one water compartment visible to MRI (Sec tion 1.4). Multiple water compartments within one imaging element give rise to a multi-exponential relaxation process given by the equation: rTmax y{U) = / s(T)exp(-U/T)dT (1.15) where y(t{) is the signal intensity at echo time ti, s(T) is the amplitude of the component at T2 = T. Equation 1.15 is a Fredholm equation of the first kind. A Fredholm equation [26] is defined as: g(x) = h(x)f(x)+ f K(x,y)f(y)dy (a < x < b) (1.16) J a and is of the first kind if h(x) = 0. Fredholm equations can be solved [26] by linear and non-linear fitting methods. Non-linear methods require a starting model and have inherent difficulties at converging to the global minimum [27]. Linear inversion techniques [28] assume a large number of T and solve for the amplitudes s(T), most of which are 0. 8 CHAPTER 1. INTRODUCTION 1.6. T2 DECAY CURVE ACQUISITION To solve integral equations on the computer, Equation 1.15 is discretized: M s(T) = Y,Sj8(T-Tj) 3=1 where M is the number of T2 components and 6 (x) is the Dirac delta function defined as: y 0, otherwise. Substituting Equation 1.17 into Equation 1.15 gives: y(U) = ^2sjexp{-U/Tj) i = l...N: (1.19) 3=1 Equation 1.19 is a matrix equation which can be solved by linear inversion techniques (discussed further in Section 1.7). 1.6 T2 Decay Curve Acquisition Decay curve data, y, is acquired with a multi-echo, spin-echo pulse sequence. Two types of spin: echo pulse sequences were used in this work: an optimized, single slice spin-echo and a conventional spin-echo, slice selective pulse sequence. A standard MRI pulse sequence is described by gradients along three orthogonal axes Gx, Gy, and Gz; the RF pulse train; and the acquired data. The acquisition ofthe data and the relationship to the image is discussed (Section 1.6.1). Each pulse sequence will be discussed (Sections 1.6.2 and 1.6.3), as well as noise considerations (Section 1.6.4) and artifacts (Section 1.6.5). 1.6.1 K-Space Standard MRI images are created by calculating the magnitude of the data encoded in the Fourier transform space - called k-space [29]. Mathematically, the relationship between the data in image space, M (x, y), and the data in k-space, S (kx, ky), is related by: (1.18) S(kx,ky)= / / M'(x,y)exp(i x kx)exp(i y ky)dy dx (1.20) 9 CHAPTER 1. INTRODUCTION 1.6. T2 DECAY CURVE ACQUISITION Figure 1.4: Standard MRI image of a brain (left) and the magnitude ofthe corresponding k-space image (right). where S (kx, ky) is the signal at kx, ky and M (x, y) is the image space signal at x, y. The variables kx and ky describe the "distance" moved in k-space and are defined as: fc* = 7 / Gx(t')dt' (1.21) Jo ky = i I Gy{t')dt' (1.22) Jo where 7 is the gyromagnetic ratio, Gx(t') and Gy(t') are the gradient shapes, t is the duration of the gradient. Each point in the MRI image is a weighted sum over all points in k-space. The k-space formalism allows the description of all possible pulse sequence data acquisition schemes (e.g., Cartesian, spiral, concentric circles). The theory is easily extended for 3 dimensions. 1.6.2 Optimized Multi-Echo Pulse Sequence An optimized, multiple echo spin-echo pulse sequence (Figure 1.5) was used to acquire multi-echo images from which multiple component T2 measurements were calculated. The implementation of the pulse sequence used for this thesis had a train of 32 composite refocusing pulses (902-180,,-10 CHAPTER 1. INTRODUCTION 1.6. T2 DECAY CURVE ACQUISITION 90x) where the duration of the 90° hard pulse was 300 fis and the duration of the 180° hard pulse was 600 ps. Composite pulses were necessary to reduce the effect of inhomogeneous Bi [30]. The excitation pulse is a sine pulse of duration 3.2 ms. Alternating, descending gradient crushers along Z [31] further reduce the effect of stimulated echoes and signal from out of slice. Decay curve artifacts are discussed in Section 1.6.5. This pulse sequence will be annotated as MESE„ in this thesis. RF G DACQ / v 1 1 1 1 i I 1 A A A A -A A- A A '• A / A ... V V V t-\z v~ 10 15 20 25 Time (ms) 30 35 40 45 50 Figure 1.5: Optimized single-slice, multi-echo pulse sequence (MESE0). The excitation pulse is a slice-selective, sine pulse and the refocusing pulses are composite pulses 90°-180°-90°. Large alternating-descending gradient crushers on G2 were used to remove signal from stimulated echoes and signal out-of-slice. 1.6.3 Standard Multi-Echo Pulse Sequence The standard multi-echo pulse sequence (Figure 1.6) consists of a selective excitation pulse fol lowed by a train of selective refocusing pulses. This pulse sequence was used in both 2D and 3D mode. In 2D mode, both the excitation and refocusing pulses were slice selective. In 3D mode, shown in Figure 1.7, the excitation pulse was slab selective the refocusing pulses were slice selective, and phase encodes along Gy and Gz were used to encode the data acquisition in the 3D k-space. A slab selective pulse selectively excites spins within a thickness of Ns where N is the 11 CHAPTER 1. INTRODUCTION 1.6. T2 DECAY CURVE ACQUISITION number of slices and s is the thickness of each slice. In the 3D implementation used in this thesis, the phase encodes along Gz were added to the constant gradient crushers along G2. The large, constant gradient crushers along Gz were applied to reduce signal from out-of-slice spins, but did not reduce the stimulated echo artifact. This pulse sequence will be denoted by MESEC in this thesis. , J v y v y v J v y v n (1 -4 -Jl -Jl (1 n (1 n • 1 nf1-D a--n II (L n u n n II' u j V.. .j -5 0 5 10 15 20 25 30 35 40 45 50 Time (ms) Figure 1.6: Conventional multiple echo 2D spin-echo pulse sequence (MESEC). The excitation pulse and refocusing pulses are slice selective. 12 CHAPTER 1. INTRODUCTION 1.6. T2 DECAY CURVE ACQUISITION RF — -4 -A- A DACQ Slab Selective Balanced Z Phase Encode: Balanced Z Phase Encode: Balanced Z Phase Encode: 10 15 20 25 Time (ms) 30 35 40 45 50 Figure 1.7: Conventional multiple echo 3D spin-echo pulse sequence (MESEC). The excitation pulse is slab selective (gray bar centered around 0 ms) and signal is phase encoded along Gy and Gz. The slice phase encodings are added to the large crushers in a way analogous to the phase encoding already on G^. 13 CHAPTER 1. INTRODUCTION 1.6. T2 DECAY CURVE ACQUISITION 1.6.4 Noise in Decay Curve Standard imaging techniques acquire complex k-space data, then an image is calculated as the magnitude of the Fourier transform of the k-space data. The noise in the magnitude images [32] has a Rician distribution, based on the assumption that the noise on the real and imaginary channels is Gaussian. The probability density function for a Rice distribution is defined [33] as: -PRice(z) = (1.23) where a, is the "standard deviation" of the function, /i is the mean, x is the parameter, and I0 (x) is a modified Bessel function of the first kind. The two interesting implications of the Rice distribution is what happens as a function of /i/<r. °4 : Figure 1.8: Rice distribution as a function of fi/a. When /i/er is zero, for example in the black area of an MRI image, then the noise is Rayleigh distributed. When ///cr is large, for example in the object of an MRI image, then the noise approximates a Gaussian distribution. As /i/cr —>• 0, the Rician distribution becomes a Rayleigh distribution. The Bessel function 14 CHAPTER 1. INTRODUCTION 1.6. T2 DECAY CURVE ACQUISITION 2O(A«) = 1 when JJL = 0 and exp (-^f2) = exp ^—when |t/ = 0. Therefore, -PRice(z) = -^exp '2^2 = ^fcayleigh (^) (1.24) (1.25) As p/a —> oo, the Rician distribution becomes similar to a Gaussian distribution. The asymp totic expansion [34] of the Bessel function IQ(X) when x is large is: I6{x) V2 7TX J_ 1•9 1-9-25 + fe + 2!(8x)2 + 3!(8x)3 + (1.26) then for sufficiently large x, = 77!= • Therefore, considering PRice(£) around a; = fi: x P(x) — exp x exp (J,2 + x2 2a2 x2 + ^i2 a exp 2<72 7 v^Ix 2a2 exp (S) (1.27) (1.28) (1.29) (1.30) then, because x ^ fi the scaling factor in front of the exponential can be simplified leaving: P(x) = •. exp V2ira2 — -^Gaussian (^0 2a2 (1.31) (1.32) which is the probability function for Gaussian distribution. All simulations, in the following chapters, used Rician noise. The Rician noise was created as ye{U) = yj[y(ti) + ei]2 + e?,, where y is the true signal, and e\ and e2 are random numbers from a Gaussian distribution with zero mean and standard deviation a. The standard deviation, cr, for the Gaussian distribution was defined as y(ii)/SNR. 15 CHAPTER 1. INTRODUCTION 1.6. T2 DECAY CURVE ACQUISITION 1.6.5 Artifacts Influencing the T2 Decay Curve Artifacts in the T2 decay curve are non-random changes to the decay curve during data acquisition. The artifacts could result from mis-tuned hardware or limitations of the underlying assumptions. Four effects that arise during acquisition of the T2 decay curve are: 1) signal from out-of-slice, 2) magnetization transfer, 3) stimulated echoes, and 4) T2 dependence on the echo spacing r. Out-of-Slice The signal measured from a voxel includes signal from protons outside the desired slice. This is because the slice profile of an RF pulse is non-rectangular. The goal of shaped RF pulses is to create a slice profile that is close to rectangular, but there is always some ripple within the slice (pass band) and outside the slice (stop band). The goal of hard RF pulses is the opposite, the desire is to flip all protons across a sample. Therefore, whether shaped RF pulses or hard RF pulses, there is a need to suppress signal from outside the desired slice profile. The standard method to suppress signal from outside the desired slice profile is to include a large pair of gradients on Gz around each refocusing pulse [31, 35]. These gradients add a large phase to the spins outside the slice, so that the signals do not get refocused during data acquisition. Magnetization Transfer Effects The pulse sequence must be carefully designed so as to not introduce artifacts during acquisition which may affect the quantification of the myelin water fraction. The standard MRI acquisition for quantification of the myelin water is a single slice pulse sequence described in Section 1.6.2. A multi-slice acquisition would be a desirable feature, but would need to be implemented care fully. For slice-selective imaging, a magnetic field gradient is applied along Gz which spatially encodes the spins with a frequency. It was shown previously [36] that an off-resonance RF pulse (magnetization transfer [MT] pulse) will reduce the signal from water trapped between the myelin bilayers due to direct saturation. Therefore, the acquisition of multi-echo, multi-slice data may reduce the signal from the myelin water as an RF pulse on-resonance for a given slice will be an off-resonance pulse for other slices in the volume. There are three possible ways around this: 1) acquire single slice images only, 2) acquire data in which a slab is excited and multiple slices are 16 CHAPTER 1. INTRODUCTION 1.6. T2 DECAY CURVE ACQUISITION phase encoded within the slab, or 3) it may be possible to acquire multiple slices if the slices are distributed in time and space so the magnetization transfer effect is minimized. Stimulated Echoes Stimulated echoes are overlayed on top of the primary echo and are due to imperfectly refocused magnetization. Stimulated echoes are a major artifact in T2 decay curves [31]. In the presence of an RF pulse [37], dephasing magnetization will be rotated such that a portion will continue dephasing, a portion will begin to rephase and a portion will be stored along the longitudinal axis. The magnetization along the longitudinal axis does not accrue phase, therefore, when the magnetization along the longitudinal axis is rotated back into the transverse plane, it, will not have the same phase as the magnetization that was in the transverse plane. The lack of phase coherence results in magnetization refocusing at a different echo time. The stimulated echo artifact is described further in Chapter 4 T2 Dependence on Inter-Echo Spacing The standard method to measure T2 in vivo is to use a multi-echo spin-echo pulse sequence where the constant echo spacing is T (TE = 2r). The T2 estimated from a multi-echo spin-echo pulse sequence was shown [38] to be dependent on r due to iron, in the form of ferritin, in the brain (e.g., basal ganglia). Iron in the brain has strong magnetic properties leading to microscopic magnetic field inhomogeneities. Jensen et al. modeled the relaxation rate R2 (R2 = 1/T2) as R2 = R2a + AR2. The relaxation rate R2a was defined independent of microscopic field inhomogeneities, and AR2 was derived to be a function which was approximated as being proportional to re3/2 where x = ADT/R2, D is the diffusion time, and R = 3.5 [im was the radius of the iron which were modeled as spheres [38]. As AR2 is a function of T, the measured R2 would also be a function of r. R2 is several orders of magnitude smaller than the measurement error for the typical r of 5 ms. Therefore, comparing T2 from different techniques, requires knowledge of r and the underlying tissue composition. 17 CHAPTER 1. INTRODUCTION 1.7. T2 DECAY CURVE ANALYSIS 1.7 T2 Decay Curve Analysis Decay curve data acquired by either technique described in Section 1.6 must be fit to determine the underlying T2. 1.7.1 Decay Curve Inversion (Decay Curve to T2 Distribution) Decay curves acquired from voxels which contain multiple water compartments are composed of a linear combination of decay curves each corresponding to one water compartment. Figure 1.9 (top) shows three decay curves: the short T2 decay curve of yshort(U = 200exp•(—t/20), the medium T2 decay curve of ymediumW = 800 exp (—1/80), and the total bi-exponential decay curve ybiexp(i) = 200 exp {-t/20) + 800exp (-t/80). Echo times at t = 10,20,..., 320 ms. The goal of any data inversion technique is to calculate the single or multiple T2 times and amplitudes, s, for a given decay curve y. The bottom plot of Figure 1.9 shows the amplitudes, Sj, and T2 times of the T2 decay curves in the top plot of Figure 1.9. In this example, s\ = 200 at T2)i = 20 ms and s2 = 800 at T2 2 = 80 ms. The plot of s versus T2 is called the T2 distribution. Each data inversion technique may assume positions for the spikes in the T2 distribution, upper and lower bound for the positions, or non-negativity of the amplitudes. Most inversion techniques minimize a misfit based on the measured data. 1.7.1.1 Misfit Measures The misfit is a scalar measure of the difference between the predicted data y and the measured data y. Two common definitions for the misfit between the predicted and measured data are Li, which measures ||y - y||i = Y,iLi \Vi ~ Vii and L2 which measures ||y - y||2 = 2~2iLi \Vi ~ Vi?• Another misfit measure, x2, accounts for the difference between the predicted and measured data, as well as the noise in the measurement, is defined as: where N is the number of data points, yi is a measured datum, is a reconstructed datum, and ai is the standard deviation of the ith datum. The x2 misfit is a logical misfit to minimize in curve fitting algorithms, as it accounts for noise in the measurements. (1.33) 18 CHAPTER 1. INTRODUCTION 1.7. T2 DECAY CURVE ANALYSIS 900 TE (ms) 800r ••. 700 k 600H & 500 -CO c CD ° 400 -o o £ 300 -200h 100";;:I j j • j i i ! : Q\ 1 1 1 i i i i i i J 20 40 60 80 100 120 140 160 180 200 T2 (ms) Figure 1.9: Decay curves (top) and T2 distribution (bottom) of a bi-exponential white matter model. The short T2 peak at 20 ms in the T2 distribution (bottom) corresponds to the fast decaying exponential (exp (—i/20), dashed-dotted line in the top plot). The medium T2 peak at 80 ms in the T2 distribution (bottom) corresponds to the slower decaying exponential (exp (—1/80), dash line in the top plot). The two peaks in the T2 distribution (bottom) correspond to the sum of the two decaying exponentials and is the solid line in the top plot (200exp (-i/20) + 800exp (-t/80)). The echo times are t = 10, 20,... , 320 ms. 19 CHAPTER 1. INTRODUCTION 1.7. T2 DECAY CURVE ANALYSIS The expected value E(x2) = N and the standard deviation ax2 = y/2N. Solutions with x2 ^ N reproduce the noise and solutions with x2 3> N misrepresent the measurement. Intuitively, if .the only difference between the predicted and measured data is noise, then Equation 1.33 is EiLi = ££i i = N. 1.7.1.2 Non-Negative Least Squares The primary curve fit algorithm used in this thesis was a non-negative least squares algorithm. The T2 distribution is calculated from a decay curve by applying the non-negative least squares (NNLS) algorithm [39] to the decay curve [40]. The objective function NNLS minimizes is: min ||v4s - y|| (1-34) s where A\j = exp (—£j/T2 j), yi is the datum at ti, and s is the set of amplitudes over T2. Therefore, As is the predicted data and y is the measured data. To incorporate the standard deviation, ai, at the data point yi, each datum yi is divided by Gi as is the corresponding row Ai. When NNLS is used to solve this modified system, it minimizes the x2 misfit in Equation 1.33. NNLS is guaranteed to minimize Equation 1.34 [39], or the x2; if the standard deviation is divided through both sides As = y. It is interesting to note that each column of A is a mono-exponential decay curve with T2 = T2j at echo times t. Therefore, NNLS minimizes the misfit of a linear combination of mono-exponential decay curves to create multi-exponential decay curve with the minimum x2 misfit. Regularization The standard NNLS technique can be modified to apply different types of smoothing to the T2 distribution. This is accomplished by including a second term in the objective function (Equa tion .1.34). The modified objective function to minimize is: min ps-y||2 + M||#s-f||2 (1.35) s where /J. controls the contribution of the second term. If /j, = 0, then the minimization yields the least squares NNLS solution. 20 CHAPTER 1. INTRODUCTION 1.7. T2 DECAY CURVE ANALYSIS The matrix H and vector f could take on many forms depending on the type of regularization desired. To minimize the changes in the distribution, H would be defined as: H = 1 -1 0 0 1 -1 0 0 : 0 0 -1 (1.36) and f — 0. Or the matrix could be designed to minimize the amplitudes as in: H = 10 0 0 1 0 0 0 : 0 0 ••• ... o ... o 1 0 0 1 (1.37) and f = 0. 800 700 600 ^500 </> c CD D c o o 400 CL 300 200 100 20 40 60 80 100 T2 (ms) 120 140 160 180 200 Figure 1.10: T2 distribution of a two compartment model. The spike distribution (dashed-dot lines) represents the least squares solution to a bi-exponential distribution. The continuous distribution represents a regularized solution to the bi-exponential dis tribution. The area of the spike at 20 ms is the same as the area of the distribution centered at 20 ms (same for 80 ms spike and distribution). 21 CHAPTER 1. INTRODUCTION 1.7. T2 DECAY CURVE ANALYSIS Implementation The baseline effect due to the non-zero mean ofthe Rician distribution at low signal (Section 1.6.4) must be accounted for within the NNLS curve fitting [41]. To account for the baseline offset during the curve fitting, an extra column of ones is added to the A matrix and a variable B is included as one more row in the s vector, then: an a12 1 Si r -. 1 2/1 a.2i a2N S2 1 V2 032 <Z3iV — SN 1 B VM AMI C-M2 • • O-MN (1.38) The set of T2 was defined to be 120 logarithmically spaced T2 = 10.00,10.52,11.06,... , 4000.0 ms. Therefore N, in Equation 1.38 was 120 and M was the number of data collected from the pulse sequence (discussed in Section 1.6). A regularized NNLS solution was used and consisted of finding the minimum x2 solution, Xmim then regularizing the T2 distribution until a solution was found with 1.02 Xmin — X2 < 1-025 Xmin-The regularization matrix H = u-I, as shown in Equation 1.37, was used for the work in this thesis. The algorithm to determine the regularized x2 is described in Algorithm 1.1. 1.7.1.2.1 T2 Distribution Bins The T2 distribution calculated using NNLS from a tissue sample is expected to have broad distri butions due to the underlying heterogeneity of the tissue [42] and the regularization method [42, 43]. To quantify the T2 relaxation time or amplitude, the T2 distribution is typically binned. A bin is defined with a start T2 and end T2. There are several measures defined for a bin. The most important measure, in particular for the short T2, is the fraction of signal in the bin relative to the total signal: fr = (1.39) 22 CHAPTER 1. INTRODUCTION 1.7. T2 DECAY CURVE ANALYSIS Algorithm 1.1 Small norm NNLS regularization. Regularized NNLS 1: Define flow = 1.02 {X2low = fiowXmin} 2: Define fhigh = 1.025 {x2high = fhighX2min\ {Calculate the least squares solution} 3: Solve s = nnls (A, y) 4: Calculate Xmin fr°m A y>s 5: Define a vector, /u, logarithmically spaced between 10~5 and 10_ {Find m such that x\ JS between x20W and Xhigh } 6: while x2 < X20W OR x2 > xligh do 7: for all m in \x do 8: Solve s = nnls (A, y, Hi) 9: Calculate xj from A, y, s 10: end for 11: Do linear fit of x2 vs fj, 12: From linear fit, calculate fi at midpoint between xfow and x\ig 13; Solve s = nnls (A, y, Hi) 14: Calculate x2 from A, y, s 15: if X?<XL then 16: Set H to be linear spaced between Hi and M(1) 17: else if Xi > xligh then 18: Set H to be linear spaced between /i(end) and Hi 19: else 20: Do Nothing 21: end if 22: end while 23 CHAPTER 1. INTRODUCTION 1.7. T2 DECAY CURVE ANALYSIS The arithmetic mean T2 of a bin is defined as: T = ET-Tend T=Tstart s{T) T , s(T) (1.40) ET—Tenti T=T„tart and the energy for the bin is £ start For this thesis, two compartments were defined. The short T2 compartment (related to water trapped between the myelin bilayers, see Section 1.4) was defined as a bin with 10 < T2 < 50 ms. The medium T2 bin was defined as 50 < T2 < 120 ms and is related to the inter and intra-cellular water. 1.7.2 T2 Decay Curve Requirements 1.7.2.1 SNR Requirements It has been shown [44] that a signal-to-noise ratio (SNR) of 700 is required to separate T2 compo nents of one order of magnitude apart in a T2 distribution. The standard method to attain high SNR, in the MESEG sequence, is to acquire four averages for each line acquired in k-space [18]. The SNR is proportional to VN where N is the number of averages - therefore, an image acquired from four averages has twice the SNR of an image acquired from one average. The total scan time for a sequence with four averages, 128 phase encode lines and TR = 3s is approximately 26 minutes. This is a long time for a volunteer or patient to attempt to remain motionless in the scanner for one sequence. Therefore, other methods of attaining the required: SNR must be 1.7.2.2 T2 Decay Curve Sampling The T2 decay curve must be sampled sufficiently to accurately estimate the T2 components and amplitudes. Whittall et al. [45] analyzed 95 in vivo white matter decay curves and concluded that 91 of them (96%) were multi-exponential based on an F-test of the misfit. Eighty-nine of the 114 gray matter decay curves (78%) and 15 of the 30 multiple sclerosis lesion decay curves (50%) were mono-exponential. A white matter model was created and mono-exponential fits of simulated four echo data (TE = 30, 60, 90,120 ms) missed the short T2 component, underestimated the T2, and overestimated the proton density. Multi-exponential fits of the same simulated data sampled at assessed. 24 CHAPTER 1. INTRODUCTION 1.8. PHANTOMS 32 echoes (TE = 10, 20,... , 320 ms) accurately reproduced the two T2 components. Whittall et al. compared mono-exponential T2 solutions from simulated decay curves sampled with different sets of echo times. They found the T2 and proton density calculated from two different sampling schemes to be statistically different even though the underlying T2 and proton density were the same. Therefore, in vivo white matter has a multi-exponential decay and few echoes (four or less) are insufficient to calculate a robust set of T2 and amplitudes. 1.8 Phantoms Decay curve acquisition and analysis were validated using a set of mono-exponential phantoms. The set of six mono-exponential phantoms were created with Ti and T2 times (see Table 1.1) similar to the brain [46]. Nickel-chloride (NiCl2) was chosen, as the paramagnetic ion Ni is an effective Ti modulator and is relatively independent of temperature and field strength [46, 47] and remains stable over long periods of time. The nickel/agarose was mixed [48] and put into glass tubes that were vacuum sealed. The Ti were measured from a saturation recovery experiment with TR = 100, 300, 800, 2000, 3000 ms. Other parameters were TE = 11.0 ms, FOV = 16 cm, 2 averages. For each phantom, the mean signal intensity over the phantom was fit to the equation y(TRj) = p[l — exp (—TRj/Ti)]. The T2 measurements were done with a 48 echo MESE0 pulse sequence with TE = 10,20,... , 480 ms. Other parameters were TR = 3000 ms, 2 averages, and FOV = 16 cm. Phantom T2 (ms) Ti (ms) mM NiCl2 % Agarose 1 237 712 2 0.2 2 91 712 2 1 3 98 i858 0.5 1 4 24 333 5 4 5 84 342 5 1 6 26 690 2 4 Table 1.1: Ti and T2 of mono-exponential nickel/agarose phantoms. 25 CHAPTER 1. INTRODUCTION 1.9. OVERVIEW OF THESIS 1.9 Overview of Thesis Noise and artifacts which result from decay curve acquisition can confound the analysis and interpretation of T2 decay curves. The thesis focused on three aspects of the noise and artifacts as they relate to T2 decay curve analysis. First, the curve fitting algorithms used to fit T2 decay curves are strongly affected by the SNR of the acquired data. I analyzed noise reduction filters applied to the multi-echo data, to determine if post-processing was as good or better than multiple acquisitions (averaging). I compared the myelin water fractions and myelin water maps calculated from simulated and in vivo multi-echo data. For the simulated and in vivo data, I applied six noise reduction filters and compared the variability in the myelin water fraction. ; Second, the standard method to calculate the myelin water fraction from a T2 decay curve is to use a curve fitting algorithm. The standard algorithm is slow and problematic for low SNR multi-echo data. The short T2 is of particular interest as it relates to the myelin water fraction. I developed a method to linearly combine the multi-echo data, which acts as a filter in the T2 space. The myelin water fraction calculated from the linear combination method was compared to the myelin water fraction calculated from the standard curve fit technique. Third, it was previously shown that stimulated echoes encoded over the primary echoes will lead to inaccurate T2 estimation. Chapter 4 discusses a method to estimate the T2 and refocusing pulse flip angle from a decay curve that has stimulated echoes. The accuracy and precision of the estimation of T2 and refocusing pulse flip angle were assessed using noise simulations. The phantoms were scanned with the MESE0 sequence to validate the theory. Myelin water fractions and maps estimated from the MESEC decay curves were compared to those estimated from the MESE0 decay curves. Several extensions are discussed. 26 Chapter 2 Myelin Water Fraction Noise Reduction 2.1 Introduction Myelin is a fatty sheath that surrounds neurons and allows for faster propagation of electrical impulses using less energy than that required by non-myelinated neurons. Over 1,000,000 people in the US are affected by diseases that either breakdown the myelin (e.g., multiple sclerosis [MS]) or impair its initial growth (e.g., adrenoleukodystrophy); therefore, it is important to have an in vivo technique able to measure myelin. Magnetic resonance imaging inherently measures properties of water molecules associated with tissues in the brain and body. There is evidence [17, 19, 25] of three water compartments in the brain: water trapped between the myelin bilayers (10 < T2 < .50 ms), intra/extra axonal water (T2 ~ 80 ms) and water associated with cerebrospinal fluid (CSF) (T2 > 2000 ms). The multiple water compartments in normal white matter result in a multi-exponential MR signal decay as a function of echo time (TE) and the shortest T2 component may be capable of providing in vivo information about myelin. Previous work [49] showed that voxel-by-voxel decay curves from a spin-echo pulse sequence optimized to collect images at multiple echo times was able to sample the white matter decay curve. The relative proportion of water, as a function of T2 (T2 distribution), was calculated voxel-by-voxel from the T2 decay curves by a regularized non-negative least-squares (NNLS [28]) algorithm. A high signal-to-noise ratio (SNR) is required to distinguish components in the T2 distribution [44] and therefore, image processing techniques may be beneficial to reduce the variability due to random noise. The anisotropic diffusion filter was created by Perona [50] and introduced to MRI by Gerig [51]. It is a locally adaptive smoothing filter that incorporates the local image intensity gradient 27 CHAPTER 2. MWF NOISE REDUCTION 2.2. SNR AND MWF strength into the smoothing based on an input parameter, K, which is a measure of the strength of an edge in the image. Therefore, the algorithm smooths local areas of homogeneity with out blurring edges. Modifications to the method were presented by Gerig [51] to filter 3D MRI volumes and MRI volumes that contain multiple channels (e.g., multi-echo data). In MRI, the anisotropic diffusion filter has been applied to work in quantification of MS lesion volumes [52], intra-cranial boundary detection and RF correction [53], vessel visualization [54] and detection of brain contour [55]. My hypothesis is that myelin water fractions (MWFs) calculated from spatially filtered spin-echo data should yield better qualitative and quantitative results than MWFs calculated from unfiltered MRI data. A set of spatial noise reduction filters, with the property that edges are preserved, are described in Section 2.3. The filters were evaluated on a simulated image to determine the filter that reduces the noise the most (Section 2.4). I collected three sets of scans with 1, 2 and 4 averages on five normal volunteers. The myelin water fractions were calculated from filtered and unfiltered multi-echo data and the results are shown in Section 2.5. The myelin water fraction variability was compared from MRI data with different numbers of averages and with and without filtering and is shown in Section 2.6. 2.2 SNR Effect on Myelin Water Fraction The signal-to-noise ratio affects the variability and robustness of the myelin water fraction (Sec tion 1.7.2:1). The effect of SNR on myelin water fraction will be shown by a noise simulation study. The variability of the myelin water fractions, as a function of SNR, was calculated from simulated decay curves. A 32-echo decay curve was created with TE =10,20,..., 320 ms, a short water compartment with p s = 200 and T2S = 20 ms, and a medium water compartment with T2m = 80 ms and p m = 800. One thousand realizations of noise were added in quadrature for each of SNR = 50, 100, 200, and 400. The myelin water fraction, for each noisy decay curve, was calculated with NNLS. Histograms of the myelin water fraction were calculated as a function of SNR. The mean and standard deviation of the calculated myelin water fractions were plotted as a histogram (Figure 2.1). The mean and standard deviation of each ofthe histograms was: 0.1997 28 CHAPTER 2. MWF NOISE REDUCTION 2.3. FILTER DESCRIPTIONS ±0.1290 (SNR 50), 0.2126 ± 0.0828 (SNR 100), 0.2144 ± 0.0490 (SNR 200), and 0.2078 ± 0.0293 (SNR 400). The spread in the histogram was very large for SNR < 100 and the shape of the histogram was non-Gaussian. When the SNR was greater than 100, the estimated myelin water fractions were closer to the true value of 0.2 and the histogram was closer to a Gaussian. 250r 200 150 SNR 50 - -• SNR 100 -- SNR 200 SNR400 100 0.2 0.3 0.4 Myelin Water Fraction 0.6 Figure 2.1: Histograms of the myelin water fraction from simulated decay curves with SNR of 50, 100, 200, and 400. The true myelin water fraction was 0.2. Therefore, the signal-to-noise ratio must be high enough to obtain accurate and precise myelin water fraction measurements. There are two methods to obtain high SNR images:; averaging during data acquisition and noise reduction image processing. Data acquisition averaging is a standard technique in MRI to obtain higher SNR images, but the cost is a longer scan time. The scan time of the MESEG pulse sequence with four averages is 26 minutes and is considered long relative to standard clinical imaging (average scan time of 5 minutes). Therefore, it would be beneficial if post-processing of the MRI data could replace the data acquisition averaging. The rest of the chapter describes results from simulated and in vivo work to assess this possibility. 2.3 Description of Spatial Filters The post-processing of MRI data using spatial filters decreases the effective noise in the image without the need of a longer scan time. There are a large variety of spatial filters that decrease 29 CHAPTER 2. MWF NOISE REDUCTION 2.3. FILTER DESCRIPTIONS the noise (noise reduction filters) each with a different set of properties and assumptions. One class of filters has the property that the edges within an image are retained, or enhanced, as opposed to blurred. This property would be desirable in a noise reduction filter as small objects and edges would remain visible in the post-processed image. There are several spatial filters with this property: anisotropic diffusion filter, SUSAN filter, median filter and wavelet filter. Each of these filters is described in more detail below. 2.3.1 Anisotropic Diffusion Filter The anisotropic diffusion filter was created by Perona [50] and introduced to MRI by Gerig [51]. It is a non-linear, locally adaptive smoothing filter that incorporates local gradient strength into the smoothing based on an input parameter, K, which is a measure of the strength of an edge in the image. Therefore, the algorithm smooths local areas of homogeneity without blurring edges. The anisotropic diffusion filter was applied as follows. For each pixel, the signal difference was calculated in each of the four directions (E, W, N, S): DE = Ie(x + l,y)-Ie(x,y) (2.1) Dw = Ie(x, y) - Ie(x - 1, y) DN = Ie(x, y + 1) - Ie(x, y) Ds = Ie(x,y) - Ie{x,y - 1) where Ie(x, y) was the signal intensity at spatial location (x,y) in the eth echo image. Then the signal intensity at (x,y) was updated: I'{x,y) = I{x,y) + L\{cEDE-cwDw+cNDN-csDs) (2.2) where the integration constant A was defined as 1/5 for 2D data and 1/7 for 3D data [51] and cj, (where d is one of the four directions, E, W, N or S) was defined as: Cd(x,y) = exp < -i2^ N (2.3) 30 CHAPTER 2. MWF NOISE REDUCTION 2.3. FILTER DESCRIPTIONS where N was the number of echoes. The channel diffusion filter [51] was used for all filtering and the kernel defined as: exp < - >/V/iC?)2 + V/2(^)2 + ... + V7M(^)2' 2 (2.4) where Ji,... ,IM are the M echo images and of is the vector of spatial co-ordinates and V was the gradient operator [51]. For this work two spatial dimensions (2D) were used and the data was filtered across echoes using Equation 2.4. 2.3.2 Susan Filter The Susan filter [56] (Smallest Univalue Segment Assimilating Nucleus) is a non-linear spatial filter similar to the anisotropic diffusion filter. The Susan filter preserves image structure by smoothing neighbors which form part of the same region as the central pixel. The primary difference between the Susan filter and the anisotropic diffusion filter is that the Susan filter uses the Univalue Segment Assimilating Nucleus (USAN) weighting to determine similarity of signal intensities. One other important difference is the center pixel is not included in the summation over the neighborhood. This difference allows for a reduction of impulse noise - spikes in signal intensity. The primary equation governing the smoothing is: > J(r V) - Ui,J)^0,0)^X + ^y+^e^{ 2^ \ Z.(ij)^(o,o) exP \ ^ ——: . where r = \fi2 + j2, a controls the scale of the smoothing, t is the brightness threshold. 2.3.3 Median Filter The median filter is a non-linear filter based on finding the median of a set of neighbouring pixels. It is implemented as a convolution across the image and replaces each voxel by the median of the 4-connected neighboring voxels and the central voxel. The median of these voxels is calculated by sorting the voxel signal intensities in numerical order and replacing the voxel under consideration by the middle value of the sorted list. For all the work in this thesis, the median filter was applied in two dimensions to each of the 32 echo images in turn. 31 CHAPTER 2. MWF NOISE REDUCTION 2.4. SPATIAL FILTERS ON SIMULATED DATA 2.3.4 Wavelet Filter The wavelet packet filter [57, 58, 59] is a non-linear filter based on the wavelet decomposition. Each image was transformed, using the 12-point coiflet wavelet [60] into a complete tree of log2(A') levels where iV = 256. The best basis set was found by a minimum entropy criterion to find the most compact basis set which represented the original signal. The best basis set was then thresholded to divide the structure from the noise in the decomposed signal. All wavelet packet filtering was kindly processed by Dr. John Wood at the University of Southern California. 2.4 Spatial Filters on Simulated Data I created a simulated set of images with different areas of myelin water fraction. Noise was added to the set of images and each filter was applied, in turn, to the noisy images, and the estimated myelin water fraction was compared to truth. Simulated images enable the assessment of each of the filters and not in the presence of other scanner problems. As well, processing of simulated images allows results to be compared back to truth. 2.4.1 Methods A simulated 32-echo MRI data set was created which modeled important aspects of white matter in an MRI image: areas of constant myelin water fraction, areas with a low gradient in myelin water fraction, and edges from myelin water fraction to no myelin water fraction. The myelin water fraction was calculated from the filtered and unfiltered data and was compared to truth. 2.4.1.1 Simulated Image A simulated myelin water fraction image was created, with different types of edges and gradients. The image was 256 rows by 256 columns, where r is the row number and ranges from 1 at the top of the image to 256 at the bottom of the image, c is the column number and ranges from 1 at the left side of the image to 256 at the right side of the image. The pixels in the image are the true myelin water fraction, /, which ranges from 0 < / < 1. Six regions, shown in Figure 2.2, were used to define the image: 32 CHAPTER 2. MWF NOISE REDUCTION 2.4. SPATIAL FILTERS ON SIMULATED DATA A B C D 1 - F Figure 2.2: Simulated myelin water fraction image annotated with six regions. 33 CHAPTER 2. MWF NOISE REDUCTION 2.4. SPATIAL FILTERS ON SIMULATED DATA Region A: Large region of zero myelin water fraction, f(r, c) = 0 for 1 < r < 128 and 1 < c < 128. Region B: Each row was defined: f(r,c) = i h{r) 129 < c < 139 h(r)(U9 - c)/(149 - 140) 140 < c < 149 (2.6) 0 150 < c < 256 and h(r) = (r - 1)/117 for 11 < r < 128. For rows 1 < r < 10, h(r) was set to 0. Region C: Each row, 129 < r < 246, was defined as: { 0 1 < c < s(r) ~ ~ (2-7) 0.15[c - s(r)]/[117 - s(r)] s{r) < c < 118 where s(r) — r — 129 + 1. Region D: Constant region of /(r,c) = 0.15 for 129 < r < 256 and 119 < c < 128. Region E: Constant region of f(r,c) = 0.30 for 129 < r < 256 and 129 < c < 139. Region F: Each row, 129 < r < 246, was defined as: f(r,c) = { 0.30[s(r) - c]/[s(r) - 139] 139 < c < s(r) 0 s(r) < c < 256 (2.8) where s{r) = (246 - r)/117 * (256 - 140) + 140. The regions and boundaries between the regions were defined to simulate regions and boundaries between tissues in the brain. The upper two quadrants have a sharp edge between 0 myelin water fraction and a positive myelin water fraction to simulate the boundary between white matter and ventricles. The myelin water fraction gradients in Regions C and F simulated regions of changing myelin water fraction in the brain due to changes in tissue. The resulting simulated myelin water fraction image is shown in Figure 2.3. 34 CHAPTER 2. MWF NOISE REDUCTION 2.4. SPATIAL FILTERS ON SIMULATED DATA Figure 2.3: Simulated myelin water fraction image (truth) from which the 32-echo data set was constructed. A 32-echo dataset was created based on the simulated image. Each echo was defined as: y(r, c, t) = f{r, c) exp (-t/20) + [1 - f(r, c)] exp (-t/80) (2.9) where / was the fraction at location row r and column c defined in Figure 2.3 and t = 10, 20, ..., 320 ms. The 32-echo dataset, created from the simulated image, had 10 realizations of SNR =100:1 quadrature noise (Section 1.6.4) added to it. 2.4.1.2 Filtering The filters tested in these simulations were: anisotropic diffusion filter of 1 iteration (anisol), anisotropic diffusion filter of 3 iterations (aniso3), anisotropic diffusion filter of 5 iterations (aniso5), median filter (median), susan filter (susan), a median filter of the myelin water map calculated from the unfiltered data (postmedian) and a wavelet packet filter with thresholds of 46 (wavelet46), 48 (wavelet48), 50 (wavelet50), and 52 (wavelet52). Each filter was applied 35 CHAPTER 2. MWF NOISE REDUCTION 2.4. SPATIAL FILTERS ON SIMULATED DATA to each of the 10 noisy 32-echo dataset. 2.4.1.3 NNLS The small norm, non-negative least squares algorithm (Section 1.7.1.2) was applied point-by-point to each filtered dataset. 2.4.1.4 Analysis Two measures were used to assess the myelin water fraction calculated from the filtered data. First, a global error was calculated between the true myelin water fraction and the reconstructed. Second, an error was calculated for the region of 0.15 and 0.30 myelin water fraction. The myelin water fraction calculated from the noisy unfiltered and filtered images was com pared to the true myelin water fraction. The measure was defined as: where f{r,c) is the true myelin water fraction at row r and column c, f(r,c) is the measured myelin water fraction, and N is the number of pixels summed over. The global error was defined in the error measure. For each of the local measures, r, c was defined to range over the constant myelin water fractions. The mean and standard deviation of the errors were calculated over the ten realizations of noise. The mean and standard deviation were calculated for each filter and compared: . 2.4.2 Results Global Myelin Water Fraction The mean and standard deviation of the error is shown in Table 2.1. The filter that produced the lowest global myelin water fraction error was the anisotropic diffusion filter (5 iterations) followed closely by the susan filter. There was very little difference between the estimated myelin water fraction and the true myelin water fraction in the areas of the simulated image where the true myelin water fraction was 0. In the areas of the true myelin water fraction where the myelin (2.10) such that it did not come within 5 pixels of the border, so the boundary does not create a bias 36 CHAPTER 2. MWF NOISE REDUCTION 2.4. SPATIAL FILTERS ON SIMULATED DATA Filter Global Error Measure none 1.964 (0.005) anisol 1.208 (0.009) aniso3 0.758 (0.012) aniso5 0.637 (0.011) median 1.175 (0.014) susan 0.741 (0.013) postmedian 1.340 (0.010) wavelet46 0.723 (0.016) wavelet48 2.504 (5.809) wavelet50 0.818 (0.084) wavelet52 1.048 (0.008) .Table 2.1: The mean global error in myelin water fraction over the 10 noisy data sets for each filter. The standard deviation is shown in brackets. water fraction was greater than 0, the difference between the true fraction and the estimated fraction was relatively constant, that is, no bias with increasing true myelin water fraction. The wavelet48 filter was the only filter to increase the global RMS error and this is likely due to the threshold chosen. Error in 15% Myelin Water Fraction For this section, the error was calculated within the section in the lower half of the image toward the mid-line that was constant at 15% myelin water. The results are shown in Table 2.2. Filter Error Measure none 2.623 (0.081) anisol 1.773 (0.081) aniso3 1.038 (0.061) aniso5 0.813 (0.056) median 1.584 (0.075) susan 0.928 (0.054) postmedian 1.719 (0.086) wavelet46 1.021 (0.052) wavelet48 2.340 (4.329) wavelet50 1.208 (0.162) wavelet52 1.570 (0.077) Table 2.2: The mean error in myelin water fraction over the 10 noisy data sets for each filter for the area of 15% myelin water. The standard deviation is shown in brackets. 37 CHAPTER 2. MWF NOISE REDUCTION 2.4. SPATIAL FILTERS ON SIMULATED DATA The filter that produced the lowest global myelin water fraction error was the anisotropic diffusion filter (5 iterations) followed closely by the susan filter. The median filters did not perform well. Error in 30% Myelin Water Fraction The error was calculated within the section in the lower half of the image toward the mid-line that was constant at 30% myelin water. The results are shown in Table 2.3. Filter Error Measure none 2 891 (0 065) anisol 1 927 (0 079) aniso3 1 247 (0 113) aniso5 0 957 (0 123) median 1 733 (0 087) susan 1 124 (0 091) postmedian 1 326 (0 082) wavelet46 1 270 (0 079) wavelet48 4 030 (8 903) wavelet50 1 499 (0 147) wavelet52 1 819 (0 080) Table 2.3: The mean error in myelin water fraction over the 10 noisy data sets for each filter for the area of 30% myelin water. The standard deviation is shown in brackets. The filter that produced the lowest global myelin water fraction error was the anisotropic diffusion filter (5 iterations) followed closely by the susan filter. The median filters did not perform well. Myelin Water Maps The simulated myelin water maps were calculated from each of the filtered data sets. One example set of data is shown in Figure 2.4. 38 CHAPTER 2. MWF NOISE REDUCTION 2.4. SPATIAL FILTERS ON SIMULATED DATA Figure 2.4: Myelin water maps of the true fraction (top left), unfiltered (top middle), susan (top right), anisol (second row, left), aniso3 (second row, middle), aniso5 (second row, right), median (third row, left), postmedian (third row, middle), wavelet46 (third row, right), wavelet48 (fourth row, left), wavelet50 (fourth row, middle), wavelet52 (fourth row, right). 39 CHAPTER 2. MWF NOISE REDUCTION 2.5. SPATIAL FILTERS ON IN VIVO DATA 2.4.3 Conclusions The simulations provided an excellent environment to assess the filters as a function of noise and compare the estimated myelin water fraction back to the true myelin water fraction. Each of the . filters reduced the RMS error except for the wavelet48 filter. The anisotropic diffusion filter of five iterations was the filter that produced the smallest error measures. 2.5 Spatial Filters Applied to In Vivo Data In vivo MRI images of the brain are more complex, than simulated images, as they contain flow artifacts, truncation artifacts, and head motion. A study similar to the simulation study was done based on a set of MRI scans of five volunteers. The myelin water fraction was calculated from a set of regions for each volunteer and compared as a function of the filters. 2.5.1 Data Acquisition and Analysis An optimized multi-echo (ME), spin-echo pulse sequence (descending, alternating crushers [31] and composite 180° pulses) was used to collect 32 echoes at TE = 10,20,... , 320 ms from a single axial slice through the genu and splenium of the corpus callosum. Five normal volunteers had three sets of ME scans: 1) a scan of 4-averages (28 minutes), 2) a scan of 2 averages (13 minutes) and 3) a scan of one average (6 minutes 30 seconds). An N-average image is acquired by collecting each line of k-space N times (see Section 1.6.1). The volunteers remained in the magnet in the same position for all three scans. Other MRI parameters were FOV = 24 cm, TR = 3,000 ms, slice thickness = 5 mm, 256 x 128 frequency/phase encoding, and pixel size of 0.94 x 0.94 mm. All data was collected on a 1.5T GE Signa Horizon system with EchoSpeed gradients (General Electric Medical Systems, Milwaukee WI). Regions were drawn on the TE =80 ms image in white matter, gray matter and CSF areas. The myelin water fraction was calculated voxel-by-voxel over each slice and the mean and standard deviation of the myelin water fractions were calculated for each region. Each 32 echo data set was filtered with: No filter (none), Anisotropic diffusion filter 1 iteration (anisol), Anisotropic diffusion filter 3 iteration (aniso3), anisotropic diffusion filter 5 iterations (aniso5), Susan filter (susan), median filter(3x3) (medf ilt), Median filter applied to MWM from 40 CHAPTER 2. MWF NOISE REDUCTION 2.5. SPATIAL FILTERS ON IN VIVO DATA unfiltered data (3x3) (postmed), Wavelet filter (20) (wavelet20), Wavelet filter (25) (wavelet25), and Wavelet filter (30) (wavelet30). 2.5.2 Results The mean and standard deviation of the myelin water fraction were calculated over each region and compared to the mean and standard deviation from the unfiltered data. The mean change in myelin water fraction over all regions, filters and volunteers was -0.3121 with a standard deviation of 0.6793. Therefore, there was very little bias introduced in the mean myelin water fraction as a result of filtering the multi-echo data. Filter fWM genu postWM postIC spl en none 9.6 (4.8) 15.7 (4.1) 12.1 (4.1.) 12.1 (5.7) 15.0 (5.2) anisol 8.9 (3.6) 15.5 (3,6) 11.7. (3.3) 12.2 (6.2) 14.9 (4.8) aniso3 8.1 (2.2) 15.3 (3.5) 12.0 (2.6) 12.9 (6.8) 14.7 (4.0) aniso5 7.8 (2.0) 15.2 (3.0) 11.7 (2.4) 13.0 (6.8) 14.4 (3.4) susan 8.5 (3.9) 15.5 (3.6) 12.3 (3.1) 12.1 (6.0) 15.1 (5.0) medfilt 8.7 (3.0) 15.2 (3.9) 12.2 (3.7) 12.9 (6.6) 16.0 (3.8) postmed 8.8 (2.2) 15.8 (2.9) 12.2 (2.6) 11.5 (4.6) 15.2 (3.0) wavelet20 8.5 (3.4) 15.1 (4.3) 12.0 (3.4) 12.4 (5.0) 14.8 (5.8) wavelet25 8.6 (3.4) 15.3 (4.2) 12.0 (3.2) 12.6 (5.7) 15.0 (5.5) wavelet30 8.8 (3.6) 15.4 (4.2) 12.2 (3.6) 12.3 (5.8) 15.1 (5.1) Table 2.4: Mean myelin water percent for four anatomical regions (standard deviation shown in brackets) from one example volunteer. The regions were: caudate nucleus (caud), frontal white matter (fWM), posterior internal capsules (postIC), and the splenium of the corpus callosum (splen). One measure for assessing the filter is to count the number of times each filter resulted in the myelin water fraction with the smallest standard deviation. This measure was done for all regions and for white matter regions only and the results are shown in Table 2.5. Over all the regions, the myelin water fraction, calculated from the anisotropic diffusion filter (5 iterations), had the smallest standard deviation. 41 CHAPTER 2. MWF NOISE REDUCTION 2.5. SPATIAL FILTERS ON IN VIVO DATA Filter All Regions WM Regions none 9 0 anisol 0 0 aniso3 1 0 aniso5 50 35 susan 0 0 medfilt 4 1 postmed 33 17 wavelet20 5 1 wavelet25 o 0 wavelet30 0 0 Table 2.5: The number of times each filter had the smallest myelin water fraction standard de-viation over all the regions (middle column) and over the white matter (WM) regions (last column). 42 CHAPTER 2. MWF NOISE REDUCTION 2.5. SPATIAL FILTERS ON IN VIVO DATA Figure 2.5: Myelin water maps of the unfiltered (top left), susan (top right), anisol (second row, left), aniso3 (second row, middle), aniso5 (second row, right), median (third row, left), postmedian (third row, middle), wavelet20 (third row, right), wavelet25 (fourth row, left), wavelet30 (fourth row, right). 43 CHAPTER 2. MWF NOISE REDUCTION 2.6. SPATIAL FILTERING VS AVERAGING 2.5.3 Conclusions The anisotropic diffusion filter applied five times is the best filter to use to minimize the standard deviation in the myelin water fraction. The resulting myelin water map is smoother and less blurred then the original myelin water map that is median filtered. The data that is input into the NNLS routine should be as good as possible and therefore the data should be filtered with the anisotropic filter rather than applying the median filter afterward. 2.6 Spatial Filtering vs Averaging Myelin water signal, in the brain, was calculated from an optimized multi-echo spin-echo MRI scan by fitting the decay curves using a non-negative least squares algorithm. The anisotropic diffusion filter is an adaptive smoothing filter that reduces the variability in homogeneous regions without blurring edges. This filter was applied to the multi-echo data to decrease the local variability and therefore increase the local sigrial-to-noise ratio. Forty regions, from five volunteers, were analyzed by three methods of computing the myelin water fraction. A consistent decrease in the myelin water fraction variability with no bias in the mean was found for the 40 regions. The myelin water fraction of white matter was more contiguous and had fewer "holes" than maps from unfiltered echoes. 2.6.1 Methods The in vivo data acquisition is described in Section 2.5.1. 2.6.1.1 Data Filtering The three multi-echo (ME) scans from each volunteer were filtered with 0, 1, 3 and 5 iterations of the anisotropic diffusion filter with parameter K — 3a where a was the standard deviation in the noise on the TE = 10 ms image (measured in the air outside the head) divided by y/2 — -K/2 (for normalization ofthe Rician distributed noise [61] at zero signal). The channel diffusion filter [51] (Section 2.3) was used for all filtering. 44 CHAPTER 2. MWF NOISE REDUCTION 2.6. SPATIAL FILTERING VS AVERAGING 2.6.1.2 Data Analysis Eight regions-of-interest (ROIs) (number pixels average/minimum/maximum was 123/104/168) were drawn in white matter on the 4-average, TE = 80 ms scan of each of the five volunteers: frontal white matter (left and right), internal capsules (left and right), posterior white matter (left and right), genu of the corpus callosum and the splenium of the corpus callosum. The average SNR over all ROIs was 199:1. The T2 distribution was calculated by the regularized NNLS algorithm [40]. The MWF was defined as: where s(T2) was the amplitude of the T2 distribution calculated by NNLS and the maximum T2 was 2000 ms. Two sets of analysis were performed: 1) MWFs calculated per ROI and 2) MWFs calculated voxel-by-voxel (myelin water maps). For the first analysis, three myelin water measures were calculated for each ROI: 1. Average the decay curves in the ROI and apply the NNLS algorithm to the averaged decay curve then calculate the MWF (average-invert). 2. Apply the NNLS algorithm to each decay curve in the ROI, calculate the MWF from each result then average the MWFs (invert-average). 3. Apply the bootstrap algorithm as defined below (bootstrap). The standard MWF measure was defined to be the average-invert, ROI analysis applied to the 4-average, unfiltered scan for each volunteer because simulations showed least bias and lowest standard deviation versus other approaches. , The second analysis was on myelin water maps which were composed of myelin water fractions calculated from T2 distributions calculated voxel-by-voxel by the NNLS algorithm (no averaging of neighboring decay curves). The number of "holes" in an ROI of the myelin water map was defined as the number of voxels in a white matter ROI with a MWF less than 0.03. All data was processed using in-house software and Matlab (The Mathworks, Natick, Mass.) (2.11) 45 CHAPTER 2. MWF NOISE REDUCTION 2.6. SPATIAL FILTERING VS AVERAGING 2.6.1.3 Bootstrap Theory The bootstrap algorithm [62] is a means of estimating confidence intervals and other statistical measures on a set of data. Consider an ROI containing N decay curves. The bootstrap method essentially simulates having many ROIs each containing N decay curves. The standard deviation of the mean MWFs from all the realizations is the standard error of the mean. The bootstrap algorithm is applied as follows: A single bootstrap resampling randomly selects N (same N as defined above) decay curves with replacement from the ROI. Some voxels may be. represented more than once and some may be missing. The mean decay curve is formed from this sample and the MWF found. The bootstrap resampling is repeated 1,000 times and the mean and standard • deviation of the corresponding MWFs found. This mean is very close to the simple mean from . the N decay curves in the ROI. The standard deviation is an estimate of the standard error of the MWF mean. 2.6.2 Results 2.6.2.1 Effect of Filtering Spin-Echo Images To confirm the anisotropic diffusion filter did not introduce a bias in the signal intensity (SI) of the ME data, I compared the mean and variance of the intensities of the anisotropically filtered ME data in the 40 ROIs to the original unfiltered ME data. The maximum SI change in the ME data on the TE = 10 ms image over all ROIs and averages was 0.24%. The mean SI changes, over all ROIs, as a function of the smoothing iterations, were not statistically significant (p > 0.19). The average change in standard deviation of the signal intensities was -1.4% (1 iteration), -4.0% (3 iterations) and -6.0% (5 iterations) over all ROIs and averages relative to the unfiltered spin-echo data. The standard deviation decreased for every ROI but was significant (ratio of variances using the F test) in only 10 ROIs in the 1-average data, 3 ROIs for the 2-average data and 0 ROIs for the 4-average data (all with 5 iterations). 2.6.2.2 Bootstrap Validation of Mean Myelin Water Fraction The bootstrap algorithm was used to estimate the myelin water fraction variability over an ROI. To confirm the mean myelin water fractions were the same as calculated by the bootstrap method a two-tailed t-test was used to compare the MWF from the bootstrap method, /(,, to the "gold 46 CHAPTER 2. MWF NOISE REDUCTION 2.6. SPATIAL FILTERING VS AVERAGING standard" average-invert, fai. The statistic D — h , fai x 100 was calculated for each set of Jai iterations of the filter for each of the averages, and was compared to zero. In every ROI the MWFs calculated from the bootstrap method were found to be statistically the same as those calculated by the average-invert method (minimum p=0.11). For the rest of the work the mean MWF within an ROI was calculated using the bootstrap method unless otherwise stated. 2.6.2.3 Myelin Water Fraction 2.6.2.3.1 Four Average Data The standard MWF measure (MWF calculated from the 4-average, unfiltered ME data) was compared to the MWF calculated from the 4-average, filtered ME data to determine if the MWF measurement had less variability after the application of the filter to the 4-average ME data. The average relative change in mean MWF for the ROIs, calculated with the bootstrap method from the 4-average data, was 1.6% for 1 iteration, 3.6% for 3 iterations and 4.8% for 5 iterations. The mean MWF decreased in 20 of the 40 ROIs for each of 1, 3 and 5 iterations relative to the unfiltered. The number of ROIs with a significant change in the mean MWF as a function of iterations was 6 (1 iteration), 26 (3 iterations) and 28 (5 iterations) based on a t-test of the difference in mean MWF. The mean of the magnitude of the relative decrease in the standard deviation of the MWF over the ROIs was 13% (1 iteration), 32% (3 iterations) and 44% (5 iterations) and the number of ROIs that had a significant decrease in the variance was 21 (1 iteration), 36 (3 iterations) and 40 (5 iterations). 2.6.2.3.2 Two Average Data The 2-average ME data can be collected in half the scan time of the ME data required for the standard (4-averages); therefore, if the MWFs calculated from the 2-average, filtered ME data are as robust as those calculated from the standard, then the 2-average, filtered ME data could be used for future work (which would be a significant savings in scan time). The mean relative difference between the MWF calculated from the 2-average (with and without filtering) multi-echo data and the standard MWF over the 40 ROIs was 26%. Only one 2-average, 5 iteration ROI had a mean MWF statistically the same as the standard. The myelin water map was created by applying NNLS voxel-by-voxel to calculate the myelin 47 CHAPTER 2. MWF NOISE REDUCTION 2.6. SPATIAL FILTERING VS AVERAGING water fraction therefore it was important to compare invert-average to the standard average-invert. The relative difference in mean MWF between the invert-average method (based on the 4-average MWF data) and the standard MWF was 29% (unfiltered), 26% (1 it eration), 20% (3 iterations) and 17% (5 iterations). Only 1/4 ofthe 40 ROIs calculated by the invert-average method were statistically the same as the standard MWF (based on a two tailed t-test). The invert-average myelin water fractions were greater than those calculated using the bootstrap method in 72% of the 40 ROIs. 2.6.2.3.3 Qualitative Myelin Water Maps Figure 2.6 shows the myelin water maps of anisotropically filtered spin-echo MRI scans. The myelin water maps calculated from the 4-average multi-echo data with 5 iterations of the filter were more contiguous in white matter. The average number of holes calculated in each of the ROIs decreased as a function of the number of iterations and as a function of averages (Table 2.6). There was an average of only 2.5 holes in the ROIs in the myelin water maps calculated from the 4-average data filtered with 5 iterations. Similar results were found if the definition of a hole was a voxel with a MWF less than 0.02, 0.04 or 0.05. Averages Unfiltered 1 Iteration 3 Iterations 5 Iterations 1 20 (8) 19 (9) 17 (9) 14 (9) 2 15 (9) 14 (9) 13 (9) 12 (10) 4 8(6) 6(5) 4(3) 3(3) Table 2.6: Number of holes per ROI as a function of iterations of the anisotropic diffusion filter (standard deviations shown in brackets). 48 CHAPTER 2. MWF NOISE REDUCTION 2.6. SPATIAL FILTERING VS AVERAGING Figure 2.6: Myelin water maps calculated from unfiltered (top left), 1 iteration (top right), 3 iterations (bottom left) and 5 iterations (bottom right) of the anisotropic diffusion filter applied to the 4-average, spin-echo data. Note the fewer "holes" in the white matter and more continuity through the maps as the number of iterations is increased. 49 CHAPTER 2. MWF NOISE REDUCTION 2.6. SPATIAL FILTERING VS AVERAGING 2.6.3 Discussion The anisotropic diffusion filter was applied to the original multi-echo data although it could have been applied to the myelin water maps instead. The curve fitting algorithm to compute the myelin water fractions from the decay curves is a non-linear algorithm and, as all curve fitting algorithms, requires a high SNR. Therefore, it was important to apply the curve fitting to data that was the best I could provide which required the application of the anisotropic diffusion filter before the curve fitting. I confirmed no bias in the ME SI was introduced by the anisotropic diffusion filter (a bias could result in artificial changes in the amplitudes or T2 times of the peaks in the calculated T2 distribution) and the standard deviation decreased for each ROI as more iterations of the filter were applied. The number of significant decreases in the standard deviation was more.in the 1-averaged ME data as a result of the lower signal-to-noise ratio. As K was defined to be 3cr, the number of iterations of the filter was the only free parameter. I applied seven and nine iterations of the filter but found both the filtered ME data and myelin water maps were "over filtered" as local areas of smoothly changing intensity became iso-intense. Five iterations of the filter offered a reasonable trade-off between a decrease in local standard deviation without "over filtering". Therefore, the decrease in the standard deviation, without introducing a bias in the Sis should lead to better qualitative MWF results. The MWF calculated from the 4-average, filtered ME data was expected to be more homo geneous than the MWF calculated from the 4-average ME data (standard) because the local variability in the SI of the ME data decreased as more iterations of the filter were applied. As expected, there was only a small change (< 5%) in the mean MWF for all the ROIs. There was no discernible bias in the MWF with filtering because half the ROI's MWF decreased and half increased. The anisotropic filtering of the ME data was shown to improve the MWF ROI measurement compared to the standard. I expected that the mean MWF calculated from the 2-average, multi-echo data with 5 itera tions of the anisotropic diffusion filter would be as accurate as the standard mean MWF calculated from the 4-average, unfiltered spin-echo data. The results did not corroborate this. To understand this, I analyzed MWF histograms in each ROI (for example, see Figure 2.7). The MWF in the 4-average, unfiltered ROI data had a distinct peak (top right plot in Figure 2.7) corresponding 50 CHAPTER 2. MWF NOISE REDUCTION 2.6. SPATIAL FILTERING VS AVERAGING 0.2 1 Average 0.2 2 Averages 0.2 0.4 I .1 0.2 4 Averages 0.4 Figure 2.7: Histograms of myelin water fraction calculated from the 1-average (first column), 2-average (second column) and 4-average (third column) spin-echo data with 0 iterations (first row), 1 iteration (second row), 3 iterations (third row) and 5 iterations (fourth row) of the anisotropic diffusion filter. 51 CHAPTER 2. MWF NOISE REDUCTION 2.6. SPATIAL FILTERING VS AVERAGING to a true MWF that did not vary spatially in the ROI. The MWF histogram of an ROI from the 2-average unfiltered data (top, center) had a broader peak than the histogram of the 4-average unfiltered. The 1-average, unfiltered myelin water fraction histogram had a peak with a mean value higher than the 4-average and, as well, a wider distribution. There are two possible reasons for the broad peak in the histogram, the first due to noise in the original multi-echo data and second due to a variability in the underlying myelin water content in the white matter across the ROI. To verify that the wider distribution was due to noise, I created a white matter model with 10% of the signal at 20 ms and 90% of the signal at 80 ms and then added 1,000 realizations each of SNR =50, 100 and 200 noise. The NNLS algorithm was applied to each model with noise and the mean, standard deviation and skew were calculated for the area under the curve between 0 ms and 50 ms normalized by the total area in the distribution. Based on these simulations, I found similar broad peaks in the 1- and 2-average histograms were due to the noise. In the MRI data the filtering did sharpen the peak in the 1-average ROI but did not shift it lower. I believe that this was a result of using the non-linear curve fitting technique on low SNR data. The invert-average method has been used to create the myelin water maps and there fore is important to validate compared to the standard. Even though there was a large rel ative difference in the mean MWF (16% for 5 iterations for the 4-average data) between the invert-average and bootstrap methods, the relative difference decreased from over 29% for the unfiltered invert-average to 16%. The MWF calculated from the invert-average method tended to be closer to the standard as more iterations of the anisotropic filter were applied to the original ME data. Therefore, the myelin water maps calculated from filtered ME data do tend to be representative of the local MWF. The myelin water fractions were greater in 72% of the 40 ROIs as calculated by the invert-average method compared to the bootstrap method (on the 4-average data). This apparent increase in the myelin water fraction was due to the non-linear NNLS algorithm applied to lower SNR data in the invert-average method. The anisotropic diffusion filter applied to the 4-average, ME data resulted in myelin water maps that had better continuity in all white matter regions. The number of holes within an ROI in the myelin water map decreased in most cases but there were several cases when a hole increased in size. If several voxels clustered together in the original MRI data had a decreased SI then several iterations of the anisotropic filter could decrease the SI of other adjacent voxels and may result in a filtered image that have more local voxels with a lower SI. The NNLS algorithm 52 CHAPTER 2. MWF NOISE REDUCTION 2.7. CONCLUSIONS would then calculate the myelin water fraction as being decreased in the local area. 2.6.4 Conclusions The hypothesis was that the anisotropic diffusion filter applied to spin-echo data would result in more precise quantitative myelin water fractions and better quality myelin water maps. Myelin water fractions calculated from the 1-average and 2-average data, with and without application of the filter, had large biases. There was no bias introduced by the filter but there was a large decrease in the standard deviation for the 4-average images. Therefore, the anisotropic diffusion filter was excellent for 4-average data and reduced the standard deviation by 40% which, effectively, was similar to collecting an 8-average scan. 2.7 Conclusions Reliable estimation of the myelin water fraction is dependent on the underlying SNR of the multi-echo data. The current method to acquire sufficient SNR is to collect multiple averages of the multi-echo data, but this results in a long scan time. A second approach to have higher SNR data is to apply a noise reduction filter to the multi-echo data. The advantage is the filter can be applied after acquiring the data. There are many types of noise reduction filters, but I chose four, each of which designed to retain the edges in the original data. I compared the four filters on simulated and in vivo multi-echo data. For the in vivo data, I acquired 1-average, 2-average and 4-average in vivo, multi-echo data sets on six volunteers. For the simulated and in vivo data, I compared the variability of the MWF in regions. For the simulated data, five iterations of the anisotropic diffusion filter (aniso5) had the greatest reduction in the regional variability. The susan and postmedian filters were similar to the five iteration anisotropic diffusion filter. The application of the filters to multi-echo, in vivo data showed similar results. The anisotropic diffusion filter of five iterations resulted in the smallest variability over different regions of the brain. Again, the susan and postmedian filters had a similar result to aniso5. The MWF estimated from multi-echo data of different number of averages and with the application of no filter, anisol, aniso3, and aniso5 filter were then compared. The expectation 53 CHAPTER 2. MWF NOISE REDUCTION 2.7. CONCLUSIONS was the MWF estimated from filtered, 2-average data would have a similar accuracy and variability to the unfiltered, 4-average data - but this was not the case. The lack of similarity was likely due to a requirement of a lower bound of the SNR to obtain similar results to the 4-average data. The minimum lower bound on the SNR is likely a result of the non-negativity constraint of the NNLS algorithm. The noise reduction filter was good for reducing the variability of the MWF over assumed homogeneous regions. 54 Chapter 3 Myelin Water Fraction from Linear Combination 3.1 Introduction Many tissues in the human body consist of several water compartments (or components) within the same voxel. White matter, for example, has two compartments, a short compartment with T2 ~ 20 ms and medium compartment with T2 ~ 80 ms [18]. The ratio of the areas of the short to medium compartments is approximately 1:6. The short T2 component was shown [18] to be associated with myelin water and is important to quantify for such diseases as multiple sclerosis and adrenoleukodystrophy. The standard method to calculate the short T2 fraction, in vivo, is to collect a decay curve with 32 [18] or more echoes and to fit the curve with a non-negative least squares (NNLS) curve fitting algorithm (Section 1.7.1.2). The 32-echo pulse sequence, defined by MacKay [18], consisted of composite refocusing pulses and alternating descending crushers [31] to reduce the effect of stim ulated echoes. A regularized NNLS algorithm minimized the total signal of the T2 distribution, such that l.02xmin < X2 < ^-•^^Xmin where Xmin was the minimum misfit of the unregularized solution. An alternative method to calculate the short T2 area is to calculate a linear combination of the multi-echo data. The linear combination of the multi-echo data is a filter defined as: where A(T2) is the filter of interest, c is a vector of coefficients. The image corresponding to the (3.1) 55 CHAPTER 3. LINEAR COMBINATION 3.1. INTRODUCTION filter is defined as: N 7L(x) = ^c,/i(x) (3.2) i=i where 7i(x) = 2~^j=i Pj(x) exP (—f /T2;j(x)). The goal, then, is to define coefficients, c, such that A(T2) includes T2 between 10 and 50 ms and excludes all other T2. The method to calculate the coefficients c is the novel work in this chapter. 3.1.1 Linear Combination for T2 Selectivity Two groups investigated the calculation of a good set of coefficients c. Both groups used a linear algorithm to calculate the coefficients. 3.1.1.1 Calculation of Coefficients by Matrix Inversion In 1982 Macovski [63] first described the "projection image" in the context of medical imaging. The projection image was based on a linear combination of images collected from an imaging modality. The technique was used to calculate coefficients to linearly combine images so the resulting projection image enhanced a desired material and eliminated interfering materials. The method was described more fully [64] as a matrix inversion problem I = FL where I was the set of image data at a voxel L (x), F was a matrix of exponential kernels -t;/T2j(x) and L the amplitudes of the T2 distribution. This method is described further in Section 3.6.1. 3.1.1.2 Calculation of Coefficients by Backus-Gilbert Method In 1991 Whittall [65] used the Backus-Gilbert technique [66, 67] to calculate the coefficients. To solve for the coefficients using the Backus-Gilbert method, a misfit equation is defined as a function of a target filter and a linear combination of decaying exponentials. The partial derivative, with respect to the coefficients, of the misfit equation is taken and set to zero. The coefficients, that minimize the misfit, are calculated by solving the derivative of the misfit equation set to zero. The Backus-Gilbert method is described further in Section 3.7.1. Both methods, above, used linear methods to calculate the coefficients. However, use of a non-linear method to calculate the coefficients, enables more constraints to be placed on the shape of the filter, A(T2), and therefore allows a better selection. 56 CHAPTER 3. LINEAR COMBINATION 3.1. INTRODUCTION 3.1.2 Properties of a Linear Combination A linear system must obey two conditions, homogeneity and additivity [68]. • Homogeneity: If the amplitude of the input into the system is scaled by a factor A then the output from the system is scaled by a factor A. /(x) -* s(x) (3.3) A/(x) A5(x) (3.4• Additivity: If two functions are input in to the system then the output response will be a sum of the two individual input responses. /i(x) + /2(x) ^ffl(x) + <?2(x) (3.5) A system that satisfies both conditions can be considered a linear system. The two rules when taken together are referred to as the principle of superposition. A system is shift invariant if the responses to an identical stimulus presented shifted in time and space are the same except for the corresponding shift in time and space. The conditions are important as they define the relationship between the noise in the original images, and the noise in the output image, L. Noise propagates linearly: "1 = E^ (3-6) 1=1 i where Oi is the noise related to image 7, and c, is the coefficient defined above. 3.1.3 Overview I designed a non-linear algorithm to calculate coefficients which define an optimal filter for the short T2 component. The nonlinear algorithm allows a larger variety of constraints on the filter that were not possible using previous methods. As my primary focus for the work was to calcu late the fraction of water with T2 < 50 ms relative to the total water, I calculated two sets of coefficients, one set that selects short T2, cs, and a second set that selects all T2, ca. To verify the linear combination was selecting the appropriate T2 components, a set of nickel/agarose phan toms of known Ti/T2 times were scanned using a multi-echo sequence and images of the short 57 CHAPTER 3. LINEAR COMBINATION 3.2. CALCULATION OF COEFFICIENTS T2 were calculated. To validate the myelin water fraction, I simulated noisy multi-exponential decay curves and compared the myelin water fraction calculated by the new technique to that calculated by the standard NNLS algorithm. Five volunteers were scanned with the multi-echo sequence and the short T2 fraction (in white matter and gray matter regions) calculated by the linear combination method was compared to the standard NNLS method. 3.2 Calculation of Coefficients The goal for each desired filter is to calculate a set of coefficients such that the filter defined by the coefficients is unity in the range of T2 of interest and zero elsewhere. Coefficients for both filters, J4S(T2) and Aa(T2), were calculated in a similar manner, but each with slightly different sets of constraints applied to the desired filters. The algorithm described below was used to calculate the two sets of coefficients: one set, cs, that selects for 10 < T2 < 50 ms and the second set, ca, that selects T2 > 10 ms. The selection filter for the short T2 coefficients was defined as As (T2) = ^2iLi ct exP (~^i/T2) (and similarly for the total T2). The desired filter T!s(T2) is a boxcar: 1 T lower < To < Toupper AS(T2) = { _ 2_ 2 (37) 0 elsewhere and for j4a(T2): , 1 T2 > 10 ms Aa (T2) = { ~ (3.8) 0 T2 < 10 ms though neither definition is able to be constructed from a finite set of exponential decays. Similarly, the desired filter for ^4a(T2) is 1 for T2 > 10 ms (a step function). Therefore, the goal was to create a filter, J45(T2), as broad as possible and near 0 for T2 > 80 ms (to suppress the signal from the large peak at T2 ~ 80 ms). Let a T2 weighted image Ii be defined as: M 7i(x) = kp(x) J^Sj-exp [-<i/T2j-(x)] + a, (x) (3.9) i=i where p (x) is the proton density, k represents RF inhomogeneity and other factors, s, is the frac-58 CHAPTER 3. LINEAR COMBINATION 3.2. CALCULATION OF COEFFICIENTS tional contribution of component T2j and at (x) is the noise. A voxel-by-voxel linear combination of iV T2 weighted images was defined as: TV 4(x) = ^/i(x) (3.10) where Cj are the weights. The noise variance in II is crj = 2~2iLi clai-An effective criterion for calculating the coefficients was to maximize the SNR at a specified T2,0: SNR/L = ^W^i^xP(-^ZlM). (3.n) The inverse of Equation 3.11 was used as the objective function to minimize for the calculation of each set of coefficients. 3.2.1 Coefficients for Short T2 Coefficients, cs, for short T2 selection were calculated by minimizing the inverse of the SNR (Equation 3.11) based on a set of assumptions and constraints. If we assume the noise variance, of in Equation 3.11, is constant for each image, then A;, p(x) and ai can be incorporated into a constant SNRn and Equation 3.11 becomes: SNR< = 8NB.EE•••*«P(-VT.*) (3.12) T2 0 was defined to be 25 ms (as opposed to 20 ms) to force a broader filter in the short T2 range of interest. Four constraints were necessary in the minimization procedure for calculating cs: 1. AS{T2) > 0 for 10 < T2 < 70 ms Force the filter AS(T2) to be positive in the range of interest for the short T2. Previous to adding this constraint, the algorithm tended to find a minimum where AS(T2) fluctuated positive and negative in the range 10 < T2 < 70 ms. 2. AS(T2) < 5 for T2 > 80 ms where 6 = max [A5(T2)] /100 This is the key constraint to suppress signal from T2 outside of the range of interest. Intu-59 CHAPTER 3. LINEAR COMBINATION 3.2. CALCULATION OF COEFFICIENTS itively, the magnitude of the filter response for T2 > 80 ms, relative to the filter response for T2 = 15 ms, must be much less than a factor of six, as the area of the medium T2 component in white matter is approximately six times larger than the area of the short T2 component. The factor of 100 was chosen to provide greater suppression of the medium T2 component. 3. £cs = 0 The filter AS(T2) must go to 0 at long T2. 4. fAs.dT2 = 1 ±0.01 The area under the filter must be near 1 so there is no increase or decrease in ISL from the filter. 0.05 0.04 0.03 '"CM t 0.02 < 0.01 0 -0.01 0 50 100 150 200 T2 (ms) Figure 3.1: Constraints on the short T2 filter. The physical interpretation of the constraints are shown in Figure 3.1. 3.2.2 Coefficients for All T2 Similarly, coefficients to select for T2 > 10 ms, ca, were calculated by minimizing the inverse of the SNR in Equation 3.11 based on a set of assumptions and constraints. 5 As(T2)>0 t 60 CHAPTER 3. LINEAR COMBINATION 3.2. CALCULATION OF COEFFICIENTS The function to maximize was: SNRoEl=i<exp(-tl/T2i0) SNRa = (3.13) where T20 = 25 ms. Two constraints were included in the minimization procedure for calculating ca: 1. Aa(T2) > 0 for T2 > 10 ms Force the filter Aa(T2) to be positive for T2 greater than 10 ms. 2. |Aa(T2) - M\ < S for T2 > 10 ms where M = max L4a(T2)], 6 = M/40 This will constrain the fluctuations in the filter response to be less than 1/40 for T2 > 10 ms. 0.01 0.009 0.008 0.007 ~ 0.006 a< 0.005 0.004 0.003 0.002 0.001 0 Aa(T2)>0 50 100 T2 (ms) 150 200 Figure 3.2: Constraints on the "all" T2 filter. The physical interpretation of the constraints are shown in Figure 3.2. 3.2.3 Description of the Algorithm The minimization algorithm used to calculate the coefficients was f mincon, a standard Levenburg-Marquardt routine in Matlab (Natick, MA). The short T2 fraction was the signal calculated as 61 CHAPTER 3. LINEAR COMBINATION 3.2. CALCULATION OF COEFFICIENTS the linear combination of the short T2 coefficients with the images divided by the signal calculated from linear combination of the total T2 coefficients with the images. The constraints, incorporated into the minimization algorithm, required a T2 filter be created, at each iteration, based on the current set of coefficients. The requirements placed on the T2 filter, for both the short and all filters, required a large number of constraints, one for each point in the filter. The 5 in the above constraints were set to values consistent with what is known about the relative amplitudes of the peaks in the T2 distribution. Small S were tried but the minimization algorithm was not able to find a feasible solution. 3.2.4 Results The coefficients, cs and ca calculated based on the minimization procedure above, are shown iii Table 3.1. The corresponding filters, As (T2) and Aa (T2), are shown in Figure 3.3. t (ms) cs ca t (ms) cs ca 10 3.0820 4.8727 170 0.0576 -0.4878 20 -0.5698 -5.9451 180 -0.1574 -0.4606 30 -1.8381 -1.3997 190 -0.3481 1.1084 40 -1.9002 1.9780 200 -0.5033 1.4641 50 -1.4213 2.5442 210 -0.6275 -0.0698 60 -0.7650 1.2543 220 -0.6940 -0.1675 70 -0.1330 0.6837 230 -0.7184 -0.1988 80 0.3858 -0.4731 240 -0.6962 0.3600 90 0.7583 -1.3316 250 -0.6253 -0.1289 100 0.9756 -1.4437 260 -0.5100 0.4137 110 1.0576 -0.7389 270 -0.3503 0.3065 120 1.0266 -0.9576 280 -0.1441 0.1169 130 0.9113 -0.0542 290 0.1020 -0.0424 140 0.7384 0.1944 300 0.3769 -0.4811 150 0.5188 0.6518 310 0.6918 -0.7215 160 0.2860 0.2573 320 1.0333 -0.1254 Table 3.1: The coefficients for the linear combination method. Coefficients cs were for the short T2 filter and ca were for selecting all T2. 62 CHAPTER 3. LINEAR COMBINATION 3.2. CALCULATION OF COEFFICIENTS 1 0.9 0.8 0.7 CO c 0.6 1 0.4 0.3 0.2 0.1 0 Figure 3.3: T2 filters for short T2 components (AS(T2), solid) and all T2 components (Aa(T2), dashed). The maximum of each curve was normalized to 1 for sake of comparison. 3.2.5 Discussion Objective Function The goal of most imaging is to acquire an image with a signal-to-noise ratio as high as possible, to allow accurate interpretation. This goal was used in the calculation of the coefficients to maximize the signal-to-noise ratio of the myelin water selection. Therefore, the objective function, to maximize, was defined as the signal in the resultant linear combination image (numerator of Equation 3.12) divided by the noise amplification (denominator of Equation 3.12). The goal of the algorithm was to maximize the signal of the myelin water, T2 sa 15 ms, relative to the noise. The variable T2 0 was defined slightly longer than that of myelin water to increase the width of the filter ^4S(T2). Larger values of T20 were tried, but resulted in poor filters AS(T2) or large coefficients and therefore increased noise amplification. Model Independence The power of the linear combination technique, in addition to greatly reduced processing time, is that the solution does not depend on a particular construction technique. This is in contrast T2 (ms) 63 CHAPTER 3. LINEAR COMBINATION 3.3. LINEAR COMBINATION SIMULATIONS to NNLS which must construct a solution based on a model. Therefore, the linear combination method does not assume the number of exponentials in the underlying data. 3.3 Linear Combination Simulations A simulation was created to confirm the mean and standard deviation of the MWF calculated by linear combination were similar to the mean and standard deviation of the MWF calculated by NNLS. 3.3.1 Model A white matter model was created based on the equation: Ii=p[fsexp(-tl/T2s) + fmexp(-tl/T2m)} (3.14) where p = 1000 was the proton density, T2S — 20 ms for the short component, fm = 1 — fs and T2m = 80 ms were the fraction and T2 of the medium component. Seven values for fs were used (/s = 0,0.01,0.02,0.05,0.10,0.20,0.30), to test the linear combination method over a range of typical myelin water fractions. One thousand realizations of Gaussian noise (SNR =100) were added in quadrature to the model, and the estimated short fraction fs was calculated by both the linear combination and the NNLS algorithms. Other parameters included N = 32 echoes, t = 10, 20,... , 320 ms and TR was considered infinite. The estimated short amplitude from both techniques was compared to the true (known) short amplitude by calculating fs — fs. The estimated short fraction fs was compared to the true fraction fs by both a t-test (assumes the data comes from a normally distributed population) and a Wilcoxon signed rank test (no underlying assumption of distribution). A Lilliefors modification of the Kolmogorov-Smirnov test was used to determine whether the resulting fs were normally distributed. 3.3.2 Results and Discussion The mean short T2 fractions calculated by both methods agreed with the true short fraction (Figure 3.4), but the only solution /s statistically equal to the truth was for the NNLS solution 64 CHAPTER 3. LINEAR COMBINATION 3.3. LINEAR COMBINATION SIMULATIONS with the true fs — 0.20. The mean short T2 fraction was within 0.02 of the true fraction. The short T2 fractions calculated by the linear combination method were normally distributed for each of the true fractions, but only when fs > 0.10 were the NNLS results normally distributed. On average the linear combination method calculated approximately 10,000 solutions per second and the regularized NNLS method calculated one solution every two seconds. 0.4 Figure 3.4: Mean short fraction over 1,000 realizations of quadrature noise (standard deviation bars) calculated by the linear combination method (left bar), NNLS (middle bar) and the true short fraction (right bar). 3.3.3 Conclusions The simulations were used to determine the accuracy and precision of the myelin water fraction calculated from the linear combination method to that calculated by the NNLS method. The NNLS method is the current standard method to estimate the myelin water fraction, therefore, if the linear combination method is as accurate, it would be a useful technique for fast myelin water fraction estimation. The mean myelin water fraction from the linear combination had a 65 CHAPTER 3. LINEAR COMBINATION 3.4. MONO-EXPONENTIAL PHANTOMS small positive bias of at most 0.012 compared to the truth. The small positive bias, near zero myelin water fraction, could be attributed to the Rician distributed noise. The linear combination method calculated approximately 10,000 solutions per second, compared to one solution every two seconds for the NNLS method. The linear combination method was shown to be a.very fast and accurate method for estimating the myelin water fraction from a multi-echo decay curve. 3.4 Selectivity of Mono-Exponential Phantoms The selectivity of the short T2 coefficients, in Table 3.1, was analyzed. The short T2 image was calculated as the linear combination of the multi-echo data to show the selectivity of the short T2 phantoms and the suppression of the longer T2 phantoms. 3.4.1 Methods The phantoms (Section 1.8) were scanned using the 32 echo MESE0 pulse sequence (Section 1.6.2). The images were combined linearly using the short T2 filter weights defined in Table 3.1. 3.4.2 Results and Discussion Excellent suppression of the longer T2 phantoms was found as shown in Figure 3.5. The T2 of the nine phantoms surrounding the water bottle were (clock-wise starting from bottom left) T2 = 264,101,25,23,97,363,157,77,23 ms. The short T2 map was analyzed quantitatively by calculating the mean and standard deviation in each of the nine nickel/agarose phantoms, the water bottle and the air (Table 3.2). The mean signal of the short T2 phantoms in the short T2 map was at least fifteen times higher than the mean signal of the longer T2 phantoms. The standard deviation measured from the selected regions was comparable to the standard deviation of the suppressed regions. There was a small increase in the magnitude ofthe signal from phantoms with T2 of approximately 100 ms compared to the phantoms with T2 greater than 200 ms. This was not unexpected as the T2 filter is close to zero near 100 ms, but not zero. Overall the quantitative selectivity of the short T2 phantoms was excellent. 66 CHAPTER 3. LINEAR COMBINATION 3.4. MONO-EXPONENTIAL PHANTOMS Figure 3.5: MRI image (TE =10 ms, left) of nine nickel/agarose phantoms (and the large water bottle) and the short T2 image (right) based on the linear combination of the original 32 spin-echo images. Region Mean Signal T2 (ms) 1~ -0.3 (17.1) 264 2 -12.3 (34.6) 101 3 470.9 (30.8) 25 4 433.5 (19.4) 23 5 -9.9 (15.1) 97 6 -2.7 (15.0) 363 7 -28.9 (22.6) 157 8 -29.1 (20.0) 77 9 547.1 (28.8) 23 Outside 1 0.7 (10.4) Outside 2 -1.4 (10.5) Outside 3 0.8 (10.5) Water 1 6.8 (15.9) 2000 Water 2 4.3 (16.6) 2000 Water 3 8.0 (13.6) 2000 Table 3.2: The mean and standard deviation of the signal in the short T2 image (right image of Figure 3.5). Regions one through nine are the phantoms around the outside and are numbered clock-wise starting in the bottom left hand corner. The three regions labeled "Outside" are regions drawn in the noise area. The three regions labeled "Water" were drawn in the large water phantom. 07 CHAPTER 3. LINEAR COMBINATION 3.5. NNLS VS. T2 FILTER 1.2r T2 (ms) Figure 3.6: The mean signal from each phantom plotted against the T2 of the phantom overlayed on the short T2 filter. 3.4.3 Conclusions The selectivity of the short T2 phantoms, from the linear combination method was excellent as was the suppression of the long T2 phantoms. For mono-exponential phantoms, the filter behaved as well as expected. 3.5 In Vivo Myelin Water Fraction: NNLS vs. T2 Filter The selectivity and suppression of mono-exponential phantoms was excellent. In this section, the linear combination method is applied to a multi-T2 component system. The short T2 fraction calculated by the T2 filter method was compared to those calculated by NNLS (Section 1.7.1.2). Normal volunteers were scanned using the optimized multi-echo pulse sequence (Section 1.6.2) and the myelin water maps were calculated with each method and compared. 68 CHAPTER 3. LINEAR COMBINATION 3.5. NNLS VS. T2 FILTER 3.5.1 Methods Five normal volunteers were scanned with the 32 echo MESE„ pulse sequence (Section 1.6.2). The axial slice was selected to intersect the corpus callosum. Each of the five scans were processed by the linear combination method and NNLS method. Six regions (each left and right) were selected (on the TE = 80 ms image for each scan): cortical gray matter, frontal white matter, genu of the corpus callosum, internal capsules, putamen and posterior white matter. The mean and standard deviation of the calculated short T2 fractions from each technique were compared with a Wilcoxon signed rank test. A Bland and ARman [69, 70] plot, used to compare two measurement techniques, was used to compare the myelin water fraction calculated by the linear combination method to that from NNLS. Given a pair of data, dj calculated by measurement technique 1, and d2 calculated by measurement technique 2, Bland and Altman plotted the difference of the points dj — d2, versus the mean of the points (dj + d2)/2 for all i. 3.5.2 Scan Results Myelin water maps were calculated by each technique and are shown in Figure 3.7. There was very good similarity between the two maps. Note the short T2 phantoms selected on each side of the head. The internal capsules, within the mid portion of the brain, were previously shown to have a higher myelin water fraction than most other regions of the brain. This region, in the map from the linear combination method, shows an isointensity to white matter regions close by. This is likely due to an increased T2 in the short compartment. An increase in the short T2 compartment will lead to a decreased myelin water fraction because of the shape of the filter. The mean and standard deviation of the myelin water fraction, for each region, calculated by the two techniques is shown in Figure 3.8. A Bland-Altman plot was calculated for the white matter and gray matter regions (Figure 3.9). White matter had 87.5% of the differences within 1.96a and a small bias toward small fractions calculated by NNLS. Gray matter had 90% ofthe differences within 1.96CT and a slight bias toward larger fractions calculated by NNLS. For each of the plots there was no significant correlation between difference and mean. The standard deviation was lower for the linear combination method in 85% of the white matter regions (over volunteers and regions). For gray matter, the standard deviation was lower 69 CHAPTER 3. LINEAR COMBINATION 3.5. NNLS VS. T2 FILTER Figure 3.7: Myelin water maps calculated by the linear combination method (left) and from the non-negative least squares solution (right). Data collected from a 32 echo (TE = 10 ms), single slice optimized pulse sequence. Bright circles were the same short T2 nickel/agarose phantoms as shown in Figure 3.5 in 17.5% of the regions for the linear combination method. 70 CHAPTER 3. LINEAR COMBINATION 3.5. NNLS VS. T2 FILTER 0.25 0.2r C g 1 0.15 LL C a> 0.1 0.05 0 I I I 1 LinCom NNLS 2-2.5% Figure 3.8: Myelin water fractions calculated per region from the linear combination method (light bars) and NNLS (dark bars). (Standard deviation shown as whiskers.) Regions were drawn on the left (1) and right (r) sides of: cortical gray matter (cgm), frontal white matter (fwm), genu of the corpus callosum (genu), internal capsules (ic), putamen (pu), and posterior white matter (pwm). 0.06 0.04 0.02 0 £ o 8 -0.02 c 0) CD -0.04 b -0.06 -0.08 -oj ! 1.96'SD • • • ; • • • • Mean > * i» ! ! . -1 96 SO * 0.06 008 0.1 0.12 0.14 0 16 0.18 0i Mean of Pairs 0.06, 0.05 0.04| 0.03 0.02 0.01 0 -0.01 -0.02 -0.0.3 ! > ' 1 i • 1.96 SD • • • - • i • -1.96 SD -0.01 0.01 0.02 0.03 Mean of Pairs 0.04 0.05 Figure 3.9: Bland-Altman plots for white matter regions (left) and gray matter regions (right). Most data lay within one standard deviation and little bias as a function of the size of the myelin water fractions. 71 CHAPTER 3. LINEAR COMBINATION 3.6. COEFFICIENTS BY MATRIX INVERSION 3.5.3 Conclusions Myelin water fractions calculated from linear combinations of T2 decay data is extremely fast and does not rely on assumptions about any underlying models (compared to traditional methods). Myelin water maps were very similar and showed the myelin water fraction even into the white matter folds of the cortex. 3.6 Calculation of Coefficients by Matrix Inversion In 1982, Macovski [63] defined a method to calculate coefficients for T2 selection using a matrix inversion technique. Matrix inversion does not allow constraints on the coefficients nor on the filter. Therefore, the expectation is the filter will be very broad compared to the filters defined in Section 3.2. In this section, I calculated coefficients to select for short T2 (cs, 10 < T2 < 50 ms), medium T2 (cm, 70 < T2 < 120 ms), and long T2 (c(, T2 > 2000 ms). Each set of coefficients was applied to in vivo data and the results were discussed and compared to my technique. 3.6.1 Calculation of Coefficients A matrix inversion problem was defined, I = FL, where I was the vector of image data, F was the matrix of exponential kernels [exp (—£,/T2 j)], and L was the amplitudes which correspond to the T2. The coefficients were the rows of F~l where Fij — [exp (—ij/T2j)] for t = 10, 20,..., 320 ms and T2 = 20, 80,3000 ms. The coefficients calculated are shown in Table 3.3. The coefficients for the short T2, cs, follow a similar trend of positive in the early echoes, then negative, and then positive. The T2 filters, /(T2), shown in Figure 3.3, were constructed by calculating the filter /'(T2) = (F~1)t exp (-i/T2). The short T2 filter, top left filter of Figure 3.10, should be compared with the short T2 filter in Figure 3.3 (solid line). 72 CHAPTER 3. LINEAR COMBINATION 3.6. COEFFICIENTS BY MATRIX INVERSION 1000 T2 (ms; T2 (ms) Figure 3.10: The filters calculated from the matrix inversion technique to calculate coefficients. Short T2 filter (top left), medium T2 filter (top right) and the long T2 filter (bottom left). The short T2 filter calculated by the inversion method {AS'M(T2), dashed line) compared to the short T2 filter (^4S(T2), solid line) calculated by the non-linear method (bottom right, plotted on a log scale to visualize the whole filter). 73 CHAPTER 3. LINEAR COMBINATION 3.6. COEFFICIENTS BY MATRIX INVERSION t (ms) cs c' t (ms) cs c' 10 1.72393 -0.40229 0.04356 170 0.01337 -0.07379 0.05064 20 0.49791 0.08544 -0.01328 180 0.05159 -0.09854 0.05469 30 -0.16482 0.32873 -0.03897 190 0.08533 -0.12034 0.05823 40 -0.49547 0.42999 -0.04681 200 0.11505 -0.13950 0.06134 50 -0.63314 0.45058 -0.04474 210 0.14120 -0.15632 0.06404 60 -0.66120 0.42708 -0.03749 220 0.16418 -0.17107 0.06640 70 -0.62934 0.38113 -0.02781 230 0.18434 -0.18398 0.06845 80 -0.56695 0.32531 -0.01728 240 0.20202 -0.19527 0.07022 90 -0.49114 0.26685 -0.00681 250 0.21749 -0.20513 0.07176 100 -0.41172 0.20972 0.00314 260 0.23102 -0.21372 0.07308 110 -0.33408 0.15599 0.01233 270 0.24283 -0.22120 0.07421 120 -0.26105 0.10661 0.02068 280 0.25313 -0.22768 0.07518 130 -0.19391 0.06188 0.02819 290 0.26209 -0.23330 0.07601 140 -0.13308 0.02176 0.03487 300 0.26987 -0.23815 0.07670 150 -0.07849 -0.01400 0.04080 310 0.27660 -0.24233 0.07728 160 -0.02982 -0.04573 0.04603 320 0.28242 -0.24591 0.07776 Table 3.3: The coefficients for the short, medium and long T2 selection based on the matrix inversion method of calculation. 3.6.2 Results and Discussion The short T2 filter (Figure 3.10, top left) is a very broad peak centered on T2 = 22 ms and crosses 0 at T2 = 80 ms. The primary problem with the short T2 filter is that it extends below 0. The filter which extends below 0 will decrease II (the myelin water contribution) if the voxel has contributions from T2 > 80 ms. If the system has only mono-exponential contributions for all voxels, then this is not a problem, but if multiple T2 compartments are possible per voxel, then this could lead to a mis-interpretation of water compartments 10 < T2 < 80 ms, as the myelin water will appear lower even though the short T2 was not lower. The short T2 filter, calculated from this method, is compared to the short T2 filter, AS(T2) calculated in Section 3.2.4 and is shown in Figure 3.10, bottom right. The filters are plotted on a log-scale to visualize the whole range of the filter. The selectivity of the short T2 filter is very effective if objects in the image are all mono-exponential. In this case, objects with a short T2 near 20 ms will be selected and objects with a T2 > 80 ms will be negative. The medium T2 filter (Figure 3.10, top right), and long T2 filter (Figure 3.10, bottom), both have a very broad selection. 74 CHAPTER 3. LINEAR COMBINATION 3.6. COEFFICIENTS BY MATRIX INVERSION T2 Selected Images The coefficients calculated in Section 3.6.1, were applied to a multi-echo MRI dataset. The short T2 (top left), medium T2 (top right), and long T2 filtered images (bottom left) are shown in Figure 3.11. The image at the bottom right of Figure 3.11 is the myelin water fraction calculated as the ratio of the short T2 image to the sum of the short T2, medium T2, and long T2 images. The myelin water map has very sharp edges but has a darker signal in the posterior white matter. Figure 3.11: Filtered images ofthe short T2 (top left), medium T2 (top right), long T2 (bottom left) and the myelin water fraction (ratio of the short T2 to the sum of the three, bottom right). 75 CHAPTER 3. LINEAR COMBINATION 3.7. BACKUS-GILBERT COEFFICIENTS 3.6.3 Conclusions The method of calculating the coefficients by matrix inversion is extremely fast, but lacks the ability to apply constraints on the coefficients or filter. The filter would be appropriate for use to select T2 when the underlying data are mono-exponential, but would not be appropriate when the underlying data are multi-exponential as is the case for in vivo brain. 3.7 Calculation of Coefficients by Backus-Gilbert The Backus-Gilbert method is a linear algorithm to calculate the coefficients, c, which mini mize a misfit. The misfit function incorporates a desired filter and the filter defined by c. The Backus-Gilbert method is described further and three sets of coefficients are calculated and the corresponding filters are compared to the filter AS(T2) as defined in Section 3.2.4. 3.7.1 Calculation of Coefficients The standard misfit function, to minimize, includes two terms: the first term is an error term between the constructed filter ("averaging function") and a target filter and the second term is a noise amplification term. The two terms are linearly, combined with a trade-off parameter, 6, which can be adjusted to arbitrarily weight one term more than another. The coefficients, c, were calculated to.minimize: /•oo TV <t> = cos(0) [A(T2,T2)0)-n(T2,T2,mm,T2!max)]2dT2 + Sm(0) Vc2cr2 . (3.15) where 6 is the trade-off parameter, A(T2, T2>0) = ^ exp(-TE/T2]0), and n(T2, T2)min, T2,max) is the boxcar function defined as: TTfT T \—) 1 ^ ^2<min - T2 - T2,max (o TC\ J-2,mm> i-2,max) ~ \ • (3.10J I 0 otherwise The coefficients which minimize Equation 3.15 are found by calculating the partial derivative of (j> with respect to each coefficient Cj, and setting the derivative to zero: _0 = O i = l...N (3.17) 76 CHAPTER 3. LINEAR COMBINATION 3.7. BACKUS-GILBERT COEFFICIENTS The linear set of N equations, defined by Equation 3.17, is then solved for c. 3.7.2 Results and Discussion Three sets of coefficients were calculated using the Backus-Gilbert algorithm and are show in Table 3.4. The first set of coefficients, cbgl, had a higher trade-off parameter, the second set, cbg2, had a lower trade-off parameter and third set, cba3, was designed to have a similar noise amplification factor as cs. The short T2 filters, corresponding to each set of coefficients, is shown in Figure 3.13, along with the filter created using the algorithm introduced in this paper. All four filters had a similar selectivity in the T2 of interest. The differences in the filter were around T2 = 500 ms. The filters, corresponding to the coefficients cb91 and cba3, deviated from 0 much more than the filters defined by the coefficients calculated in this paper: This deviation would introduce signal in II from water compartments with 100 < T2 < 1000 ms. The signal from T2 > 100 ms is undesirable as the filter was designed for short T2. The filter, corresponding to coefficients, cb92, had a similar response to the filter defined by the coefficients calculated in his paper, but the noise amplification was 113.4 for, cbg2 and was 29.8 for t (ms) cbgi cbg2 cbg-i t (ms) cbgl cbg2 cbgi 10 3.3504 6.8123 4.0186 3.0820 170 0.1805 -0.1123 0.1191 0.0576 20 -1.6519 -7.5974 -2.3374 -0.5698 180 0.1162 -0.2007 0.0372 -0.1574 30 -1.6415 -2.2797 -2.0012 -1.8381 190 0.0558 -0.2456 -0.0350 -0.3481 40 -1.1697 -0.1489 -1.2865 -1.9002 200 0.0006 -0.3335 -0.0984 -0.5033 50 -0.7032 0.5016 -0.6791 -1.4213 210 -0.0470 -0.3556 -0.1489 -0:6275 60 -0.3204 0.6398 -0.2204 -0.7650 220 -0.0873 -0.3552 -0.1872 -0.6940 70 -0.0287 0.6525 0.1119 -0.1330 230 -0.1196 -0.3735 -0^2113 -0.7184 80 0.1796 0.6430 0.3350 0.3858 240 -0.1420 -0.3452 -0.2212 -0.6962 90 0.3183 0.5675 0.4668 0.7583 250 -0.1561 -0.2870 -0.2187 -0.6253 100 0.3985 0.5403 0.5312 0.9756 260 -0.1611 -0.2300 -0.2021 -0.5100 110 0.4340 0.4623 0.5414 1.0576 270 -0.1575 -0.1456 -0.1741 -0.3503 120 0.4337 0.3956 0.5123 1.0266 280 -0.1456 -0.0473 -0.1314 -0.1441 130 0.4088 0.2962 0.4564 0.9113 290 -0.1259 0.0793 -0.0800 0.1020 140 0.3642 0.2009 0.3812 0.7384 300 -0.0980 0.2389 -0.0171 0.3769 150 0.3078 0.1028 0.2950 0.5188 310 -0.0634 0.3758 0.0556 0.6918 160 0.2455 -0.0174 0.2066 0.2860 320 -0.0228 0.5630 0.1385 1.0333 Table 3.4: Coefficients calculated by the Backus-Gilbert method with with a different set of trade off parameters and, therefore, properties. Coefficients calculated by the method defined in this thesis, cs, are shown for reference. 77 CHAPTER 3. LINEAR COMBINATION 3.7. BACKUS-GILBERT COEFFICIENTS 1 1 I 1 i 1 cbg1 cbg2 - - cb93 ' .... cs • , , , ,lt«* \ 1 \ >^ y -50 100 150 200 TE (ms) 250 300 350 Figure 3.12: Plot of the coefficients calculated from the Backus-Gilbert technique. Coefficients, cs, calculated in Section 3.2 are shown for reference. 78 CHAPTER 3. LINEAR COMBINATION 3.7. BACKUS-GILBERT COEFFICIENTS Figure 3.13: The filters calculated from the Backus-Gilbert technique. Filter, AS(T2), calculated based on the coefficients cs, is shown for reference. 79 CHAPTER 3. LINEAR COMBINATION 3.8. OPTIMAL ECHO TIMES 3.7.3 Conclusions Several sets of coefficients were calculated using the Backus-Gilbert method. The coefficients, and corresponding filters, were compared to those calculated by the method in Section 3.2. The filter, AS(T2), resulted in a smaller contribution for T2 > 100 ms and a small noise amplification. Better filters can be created when non-linear constraints can be used. 3.8 Optimal Echo Times Coefficients defined in the previous sections were specific to 32 echoes spaced 10 ms apart as most data, to date, were collected using an optimized, multi-echo sequence with those parameters. Both sets of parameters, number of echoes and the echo spacing, can be varied and a different set of parameters may result in a better filter for the short T2. Two sets of experiments were done to determine the quality of the short T2 as a function of the number of echoes and the echo spacing. First, coefficients were calculated for equally spaced echoes, and second, coefficients were calculated where the echo spacing was allowed to vary. A set of measures was defined to facilitate the comparison. 3.8.1 Filter Measures A set of measures, calculated from the T2 filters, were defined to allow comparison as the number of echoes and echo spacing were allowed to vary. The goal of the measures was to quantify the width of the filter for 10 < T2 < 50 ms and the T2 at which the filter reaches 1/100 of the filter maximum. The filter is defined as A(T2) = Eili °i exP(_*i/T2), where N is the number of echoes. The goal of the filter is to be as wide as possible for 10 < T2 < 50 ms and near zero for T2 > 80 ms. Figure 3.14 shows a short T2 filter. Therefore, the measure of the goodness of a filter is based on the full width at half the maximum (FWHM), which should be as wide as possible and the T2 that the curve comes down to 1/100 of the maximum height of the curve (full width at 1/100 max [FW100M]). The FW100M is important as the medium T2 peak around 80 ms is at least 6x larger than the peak around 20 ms. The FW100M measure is less important as it is determined, in part, by the constraint in the algorithm. 80 CHAPTER 3. LINEAR COMBINATION 3.8. OPTIMAL ECHO TIMES 0.045 r -0.005" 1 1 1 1 1 : ' 0 20 40 60 80 100 120 T2 (ms) Figure 3.14: The measures of the effectiveness of a short T2 filter. The full width at half the maximum should be as large as possible and the T2 where the filter reaches 1/100 of the maximum should be less than T2 = 80 ms and greater than T2 = 50 ms. 81 CHAPTER 3. LINEAR COMBINATION 3.8. OPTIMAL ECHO TIMES 3.8.2 Equally Spaced Echoes The number of echoes is expected to influence the quality of the filter. Therefore, the calculation of the coefficients was re-run for N = 1 to N = 32 equally spaced echoes. The short T2 filter is analyzed for each N and compared based on measures defined below. 3.8.2.1 Methods For the case of equally spaced echoes, the echo spacing, Ai, was allowed to vary along the coef ficients. The objective function used was the same as the one defined previously. The inverse of the objective function was defined as: with constraints -2 < Cj < 2 for i = 1... N, 10 < At < 1000 ms, and t = At, 2 At,NAt. The 3.8.2.2 Results and Discussion The echo spacing was allowed to vary to determine the best set of equally spaced echo times, the constraints on the echo spacing were implemented to be within a realistic range. The lower bound was set to be the minimum allowed by the MESE0 sequence. The maximum echo spacing was defined but not relevant as the minimization procedure used the minimum echo spacing. The coefficients were calculated given the number of echoes from N = 1 to N = 32 and are shown in Figure 3.15. Coefficients for one and two echoes did not produce an appropriate short T2 filter. Coefficients of three or more echoes were very similar. The first coefficient was positive, the following echoes were negative up to the coefficient that corresponded to TE « 80 ms, and the following coefficients oscillated positive and negative. The top plot of Figure 3.16 shows the width of the FWHM of AS(T2). The width is relatively constant from 4 echoes to 19 echoes. When more than 19 echoes are linearly combined, the width of the filter increases linearly with the number of echoes. Therefore, if having fewer echoes is more important, then one could use 5 echoes with little degradation of the filter. max c,A4 2~2j=i Cjexp (-tj/2b) (3.18) constraints on Cj were defined to constrain the noise amplification ^ 82 CHAPTER 3. LINEAR COMBINATION 3.8. OPTIMAL ECHO TIMES c CD O H— M— CD O O 400 40 0 100 Number of Echoes TE (ms) Figure 3.15: The coefficients as a function of the number of echoes and TE. The line segment is the coefficients for two echoes, and appears out of place only because of the perspective of the axes. 83 CHAPTER 3. LINEAR COMBINATION 3.8. OPTIMAL ECHO TIMES 80 o cn •I 79 o o 78 it 77 76 75_l 10 15 20 25 30 35 Number of Echoes Figure 3.16: The full width at half max (FWHM) of the filter AS{T2) as a function of the number of echoes (top) and the full width at 1/100 max (bottom). The wider the FWHM the better. 84 CHAPTER 3. LINEAR COMBINATION 3.8. OPTIMAL ECHO TIMES One interesting thing is that for evenly spaced echoes it is very important to have data included from near 80 ms. Therefore, for fewer echoes, for example N = 4, the echo spacing must be increased to have an echo near 80 ms. This means the first echo must be near 20 ms. There is a tradeoff between having an echo at shorter TE, near 10 ms, and having an echo near 80 ms. 3.8.2.3 Conclusions The objective function was maximized when the echo spacing, At, was the minimum allowed At = 10 ms. The measures defined for the short T2 filter had little change for 5 < N < 19 and got progressively better for N > 19 echoes. Therefore, if a minimum number of echoes are needed for equidistant sampling, then N = 5 should be sufficient, but if more echoes are possible, the filter would be better with more. 3.8.3 Non-Equidistant Multi-echo data can also be collected with non-equidistant echo spacing. A study, similar to the previous one of equidistant samples, was done with non-equidistant samples. 3.8.3.1 Methods The objective function was the same as that used in the other sections. The inverse ofthe objective function was defined as: with constraints — 2 < Ci < 2 for all i. There was also a constraint on the minimum spacing between consecutive echoes (10 ms) because of hardware constraints. The constraints on Ci were 3.8.3.2 Results The echo times were calculated along with the coefficients. The minimum spacing between con secutive echoes could likely be decreased slightly, but was a reasonable number. A constraint on the maximum echo time, t^, was defined to be 2000 ms as the maximum echo time affects the TR and Ti weighting. max t,c Zi=i^exp(-^/25) (3.19) defined to constrain the noise amplification 85 CHAPTER 3. LINEAR COMBINATION 3.8. OPTIMAL ECHO TIMES The non-equidistant coefficients are shown in Figure 3.17. The coefficients and echo times are similar to the coefficients and echo times calculated for equidistant echo sampling. The FWHM was approximately 2 ms wider, and the FW100M was approximately 2 ms less, when the echo sampling was allowed to vary compared to echoes equally spaced. 400 40 0 100 Number of Echoes TE (ms) Figure 3.17: The coefficients as a function of the number of echoes and TE. 86 CHAPTER 3. LINEAR COMBINATION 3.8. OPTIMAL ECHO TIMES Number of Echoes 821 1 1 r Number of Echoes Figure 3.18: The full width at half max (FWHM) of the filter AS{T2) as a function of the number of echoes (top) and the full width at 100th max (bottom). 87 CHAPTER 3. LINEAR COMBINATION 3.8. OPTIMAL ECHO TIMES The echo times calculated for four equally spaced echoes required a longer echo spacing, and therefore a later first echo, to have an echo near 80 ms. For unequally spaced echos, the time of the first echo was 10 ms, as the other three echoes could be spread out in T2 space. 3.8.3.3 Conclusions The short T2 filter was slightly better when the echo spacing was allowed vary compared to a fixed echo spacing. The echo times tended to be similar to the equidistant echo times. Therefore, to create a better short T2 filter, shorter echo times would be required, which would be difficult on an MRI machine, but would be possible on an NMR spectrometer. 3.8.4 Discussion The selectivity of the short T2 filter, A5(T2), for 10 < T2 < 50 ms is a function of the echo and corresponding coefficients. Figure 3.19 shows the T2 information in spin-echo signal acquired at echo times t = 10, 80, 100, 200, and 800 ms. As expected the later the echo time, t, the less information is available for T2 < t. Therefore, a sharp rise on the left side of AS(T2) requires signal collected at a short echo time. The goal of the coefficients was to select appropriate weightings of the signal, at different echo times, to create a selection in the short T2 range of interest. The signal from each echo contributes different information depending on the echo time of the acquisition. The 32 coefficients, cs, in Table 3.1 vary from positive to negative and back to positive. The first echo, which has a positive coefficient, is the primary echo that provides the sharp increase in the curve (upper solid line). The next six echoes have a negative coefficient and provide the change in slope from positive to negative (dot-dashed curve). The following ten echoes have a positive coefficient and bring the filter positive. Then the next eleven echoes have a negative coefficient. The effect of increasing number of echoes contributed to AS(T2) is shown in Figure 3.20. As more signal from more echoes were added to the selection curve, the selection was better at later echo times. 88 CHAPTER 3. LINEAR COMBINATION 3.8. OPTIMAL ECHO TIMES •? 0.6 • * y • /" / y s * / TE=10 TE=80 - - TE=120 TE=200 TE=800 / / / / / / / 50 100 150 200 250 300 350 T2 (ms) 400 Figure 3.19: The filter, as a function of T2, based on echo time TE = 10 ms (upper solid line), TE = 80 ms (dash-dotted line), TE = 120 ms (dashed line), TE = 200 ms (dotted line), and TE = 800 ms (lower solid line). Each filter is zero at TE = 0 ms and goes to one as T2 —> oo. 0.15 -0.05 -0.1 I 1 After 1 echo / "N\ / ^^^^.c-. ~. ~~ ~ ~~~ After 32 .echoes After 28 echoes \; '"V N. ' V s. I I 11 . . =l=l 0 50 100 150 T2 (ms) Figure 3.20: The short T2 filter as a function of a consecutive subset of echoes based on the thirty-two coefficients. Each consecutive set of echoes constrains later parts of the filter. 89 CHAPTER 3. LINEAR COMBINATION 3.8. OPTIMAL ECHO TIMES 3.8.5 Four Echo Non-equidistant Sampling Fewer echo times require less power deposition and could be interleaved so multiple slices could be acquired in such a way the MT effect could be minimal. A set of four coefficients was calculated and are shown in Table 3.5. The echo times, along with the coefficients, were calculated as part TE (ms) cs ca 10 3.1872 2.5642 30 -4.7952 -2.5133 120 2.3515 1.8698 270 -0.7539 -1.1654 Table 3.5: Coefficients for short and all T2 filter based on four unevenly spaced echoes. of the minimization procedure and are approximately logarithmically spaced. Each of the echo times are near an important region ofthe short T2 filter. The first echo time, TE = 10 ms, defines the rapid increase in the portion of the curve for 0 < T2 < 18 ms. The second and third echoes, at TE = 30 ms and TE = 120 ms, bracket a portion of the T2 filter where the filter changes from decreasing to flat. The fourth echo, at TE = 270 ms forces the tail of the filter to be near zero at long T2. The short T2 filter (Figure 3.21) calculated from the four echoes followed the T2 filter cal culated from 32 echoes. The ascending portion of the filter (0 ;$ T2 ;$ 15 ms) was exactly the same in each as that portion of the curve is determined by the TE = 10 ms contribution only. The descending portion of the filter (15 ~ T2 ;$ 100 ms) was higher for the filter based on the four echoes as there were fewer degrees of freedom to be modified to force the filter lower. The short T2 filter, from 4 echoes, was not as close to zero as the short T2 filter from 32 echoes for T2 > 200 ms. This was due to the lack of the number of data points at long echo times, which contribute to making the filter as flat as possible. Myelin water maps were calculated from the four unevenly spaced echoes and is shown in Figure 3.21 along with a myelin water map from 32 echoes. There was little difference in the myelin water maps. The original multi-echo data were filtered with five iterations ofthe anisotropic diffusion filter. The short T2 filter had a very similar response with as few as four unevenly spaced echoes. Therefore, the primary information required to create such a selective filter was required to have 90 CHAPTER 3. LINEAR COMBINATION 3.8. OPTIMAL ECHO TIMES Figure 3.21: Short T2 filter calculated from four echoes (red dashed line) and from 32 echoes (blue solid line). The height of the curve was normalized to one to allow a visual comparison of two curves. The filter is plotted on a log axis to visualize over a large range of T2. several echoes at small TE, one echo near the T2 where the filter dropped to zero and one echo at long TE. 91 CHAPTER 3. LINEAR COMBINATION 3.9. CONCLUSIONS Figure 3.22: Myelin water map calculated by the linear combination method from four echoes of the 32 echo sequence (left) and the linear combination method from all 32 echoes (right). The two bright phantoms on the right side of each head are the nickel/agarose phantoms with a T2 ~ 20 ms. 3.9 Conclusions Fast accurate estimation of myelin water is important for assessing diseases such as multiple sclerosis. The current method to quantify the myelin water fraction in vivo is to acquire multi-echo images and compute the distribution of T2 voxel-by-voxel. The myelin water fraction is calculated as the ratio of the signal of water with 10 < T2 < 50 ms to the total signal in the distribution. The calculation of the T2 distribution for each voxel in an image is a slow process. Another method to estimate the signal of water within a range of T2 is to linearly combine the multi-echo images with optimized coefficients. I designed a non-linear algorithm to calculate two sets of optimized coefficients, one set to select water with a 10 < T2 < 50 ms (short T2) and a second set to select water with T2 > 10 ms (all T2). The coefficients to select for short T2 were applied to multi-echo images of a set of phantoms. The resulting short T2 image had very good selection of the phantoms with T2 near 20 ms and good suppression of phantoms with T2 > 70 ms. The NNLS algorithm and linear combination method were applied to a set simulated, noisy bi-exponential decay curves, of white matter model. The mean myelin water fraction estimated from each of the algorithms was within 10% of the true myelin water fraction and the noise standard deviations were similar. Maps of 92 CHAPTER 3. LINEAR COMBINATION 3.9. CONCLUSIONS the in vivo myelin water fraction, from five volunteer MRI scans, were calculated by the two algorithms and had very similar detail. The linear combination of multi-echo data was an efficient algorithm to calculate the myelin water fraction based on the short T2. It was 20,000 times faster than the regularized NNLS algorithm and the myelin water fraction was as accurate and precise. The strength of the. linear combination method lies in the speed and that it does not depend on a construction method. Initial work showed fewer echoes may be required to estimate the myelin water fraction with the same accuracy and precision as the current technique. 93 Chapter 4 Non-180 Refocusing Pulses 4.1 Introduction 4.1.1 Motivation Accurate T2 estimation traditionally depends on accurate refocusing pulses. Many studies [31, 35, 71, 72] showed large deviations in the measured T2 when the refocusing pulse flip angle deviates from 180°. For example, Majumdar [71] measured biases in R2 (R2 = 1/T2) of greater than 10% for 20% deviations from an ideal pulse and biases of over 100% for deviations of greater than 35%. The biases result from a change in shape of the decay curve due to stimulated echoes superimposed on the primary echo and as the refocusing pulse flip angle a decreases from 180°. There are many methods to reduce the stimulated echo "artifacts" during acquisition of T2 decay curves. Levitt [30] proposed using composite refocusing pulses (QO^-lSO^-QOx) to correctly rephase spins in the presence of moderate Bi inhomogeneity. Zur and Stokar [73] analyzed artifacts due to inhomogeneous Bi and developed a phase cycling scheme to dephase spurious echoes. Several other groups [31, 35, 74] proposed crusher patterns to dephase spurious echoes. Each method reduces artifacts well but increases power deposition or increases the minimum echo spacing. Majumdar [71] modeled the magnetization at the ith echo, measured from a train of N re focusing pulses, to be attenuated by the ith power of a function / (a) for a pulse sequence that crushes stimulated echoes. The function / (a) was derived analytically for the first 3 echoes of pulse sequences consisting of simple RF pulses and of composite RF pulses. In 2000, Sled and Pike [75] described a method to correct inaccurate mono-exponential T2 measurements due to imperfect refocusing pulses. A multi-echo pulse sequence was used to calculate the T2 and a separate acquisition was required to calculate the Bi map. The calculated T2 was divided by f (a) to obtain.a "corrected" mono-exponential T2- Another approach by Lepage [76] corrected 94 CHAPTER 4. NON-180 REFOCUSING PULSES 4.1. INTRODUCTION the experimentally determined T2 of a gel phantom based on a Bj map calculated from a dif ferent homogeneous gel phantom. The techniques by Sled and Lepage each required a separate acquisition to calculate the Bi map. 4.1.2 Previous Work The evolution of primary, stimulated and indirect echoes acquired from a train of refocusing pulses of flip angle a is usually calculated by the application of rotation matrices to magnetization vectors [37, 71, 77, 78, 79]. Woessner [37] derived a set of equations, based on the Bloch equations, that define magnetization immediately after a hard refocusing pulse given the magnetization immediately before: F{t + tw) = F(t) cos 2 (a/2)+F*(t) sin 2 (a/2)-iMz(t) sin (a) (4.1) Mz(t + tw) = Mz(t) cos (a) - ^i[F(t) - F*(t)} sin (a) where F(t) is the dephasing transverse magnetization, F*(t) is the rephasing transverse magne tization, Mz(t) is the longitudinal magnetization, tw is the pulse duration, and a is the flip angle of the pulse. These equations describe three echo pathways: a primary echo due to the phase reversal of the magnetization; a stimulated echo due to longitudinal magnetization rotated back into the transverse plane; and the "virtual" stimulated echo due to magnetization unaffected by the first pulse, inverted by the second and later rephased to form an echo. Hennig [78] extended Woessner's work to track groups of spins with the same phase coher ence through a multi-pulse experiment. Magnetization after a non-1800 refocusing pulse is dis tributed into the Z axis or remains in the transverse plane either defocusing or refocusing. To account for this, Hennig introduced magnetization sub-states defined by a net magnetization, phase accumulation, and direction of spins (refocusing or defocusing). Constant time between refocusing pulses greatly reduced the number of possible sub-states and allowed the phase accu mulation to be grouped in an integral number of phase intervals over which phase is accumulated. Therefore, a multi-pulse experiment corresponds to an application of a rotation matrix, based on Equation 4.1, applied to each set of magnetization sub-states followed by an inter-pulse re laxation of exp(-2r/T2) for F(t) and F*(t), exp(-2r/Ti) for Mz(t), and a generating term MQ (1 — exp(—2r/Ti)) (where 2r is the echo spacing). For accurate T2 quantification, Hennig 95 CHAPTER 4. NON-180 REFOCUSING PULSES 4.1. INTRODUCTION suggested all magnetization from stimulated and indirect echoes must be removed by spoilers. On the contrary, accurate and precise T2 quantification is possible from decay curves collected with stimulated and indirect echoes superimposed on the primary echoes. 1 Mx, My F* RF Pulses Figure 4.1: Phase coherence diagram based on a train of five refocusing pulses. The phase of the magnetization would follow the lines between F\ and F* around Mx,My if the refocusing pulses were 180°. 96 CHAPTER 4. NON-180 REFOCUSING PULSES 4.1. INTRODUCTION 4.1.2.1 Mathematical Description of Woessner's Paper Woessner [37] defines the Bloch equations in terms of U (in phase component), V (out-of-phase component) and Mz (longitudinal component) as: jU(t) = -(Uo-u)V-U(t)/T2 jV{t) = {LjQ-u3)U-V(t)lT2-uiMz(t) ±M.(t) = _^LZp)+WlVW (4.2) where u>o is the nuclear precessional frequency, OJ is the frequency of the rotating coordinate system, and UJ\ = 7B1 is the frequency of the excitation pulse. The variable F(t) = U(t) + iV(t) is the rotating component of the magnetization. After a rotation of 6 = u)\tw (where tw is the width of the RF pulse): U{t + tw) = U(t) V(t + tw) = V(t) cos (0) - Mz sin (0) Mz(t + tw) = V(<)sin(0) + M2cos(6») (4.3) and F{t + tw) = U (t + tw) + iV (t + tw) = U (i) + iV (t) cos (9) - iMz sin (9) (4.4) 97 CHAPTER 4! NON-180 REFOCUSING PULSES 4.1. INTRODUCTION Equation 4.4 is then written as a function of F(t), F*(t), M(t), and M*(t): F(t + tw) = U(t) +iV(t) cos {9) - iMz (t)sin(0) (4.5) = U(t) [cos2 (9/2) + sin2 (9/2)] + iV(t) [cos2 (9/2) -sin2 (9/2)] - iMz (t) sin (9) (4.6) = U(t) cos2 (9/2) +U(t) sm2 (9/2)+ iV(t) cos2 (9/2)-iV(t) sin2 (9/2) - iMz (t) sin (9) (4.7) = U(t) cos2 (0/2) + iV(t) cos2 (0/2) + r/(i) sin2 (9/2) -iV(t) sin2 (9/2)-iMz(t) sin (9) • (4.8) = [U(t) + iV(t)} cos2 (9/2) + [C/(t) - iV(t)] sin2 (0/2) - iM2 (t) sin (0) (4.9) = F(i)cos2(0/2) + F*(t)sin2(0/2)-iMz(i)sin(0) (4.10) by using the cosine double angle formula: cos (0) = cos2 (0/2) - sin2 (0/2). (4.11) and the identity: cos2 (0/2) + sin2 (0/2) = 1. (4.12) Then, F(t + tw) and Mz(t + tw), written as functions of F(t) and Mz(t), will be: F(t + tw) = F(i)cos2(0/2) + F*(i)sin2(0/2)-iM2(i)sin(0) (4.13) F*(t + tw) = F(t)sin2(9/2) + F* (t) cos2 (912)+iM*z(t) sin (9) (4.14) Mz(t + tw) = Mz(t)cos(9)-^[F(t)-F*(t)}sin(9) (4.15) M*(t + tw) = M*(t) cos (9) - \[F(t) - F*(t)]sin(9) (4.16) 4.1.3 Implementation of Algorithm The algorithm toxcalculate the signal of each echo from a train of RF pulses is shown in Algo rithm 4.1. The first section creates the rotation matrix for the magnetization sub-states. The sub-state vector, B, is preallocated and set to —i, dephasing magnetization. Then, the algo-98 CHAPTER 4. NON-180 REFOCUSING PULSES 4.1. INTRODUCTION rithm loops over each echo, applies the rotation matrix to each magnetization sub-state, saves the echo signal intensity, applies the relaxation and then applies the phase evolution. The Ti (Algorithm 4.2) and T2 (Algorithm 4.3) relaxation are applied to the longitudinal and transverse magnetization, respectively. The phase evolution (Algorithm 4.4) updates the magnetization sub-state vector B. Each defocusing sub-state is defocused by one time unit, each refocusing sub-state is defocused one time unit, and the magnetization sub-state that forms an echo becomes a defo cusing sub-state. The phase evolution bookkeeping was simplified greatly due to the assumption the refocusing pulses were evenly spaced. The primary signal attenuation is due to the dephasing of the spins and is accounted for by the T2 decay, but there is also a weak signal dependence on the Ti relaxation time. For the primary work in this chapter, involving the brain, the Ti « 650 ms [80]. The Ti relaxation is applied to magnetization stored in the longitudinal axis by the equation exp(—r/Ti). The important aspect, for this work, is exp(—r/Ti) « exp (—5/650) = 0.9923. Therefore, the decay curve will have a very small modulation because r/Ti is not zero. The decay curve dependence on Ti could possibly be exploited as a means of calculating the Ti along with T2, a, and p. The signal intensity, y, was calculated (Algorithm 4.1) given T2, Ti, a, r, and the number of refocusing pulses. The algorithm calculated the 32-echo decay curve in approximately 25 ms on a P4 1.5 GHz computer with 1 Gigabyte of RAM. 99 CHAPTER 4. NON-180 REFOCUSING PULSES 4.1. INTRODUCTION Algorithm 4.1 Signal calculation from a train of refocusing pulses function [decay] = stimulate (alpha, num_echoes, Tl, T2, tau) % alpha = rotation vector in x'y'z frame of reference. % num.echoes = number of echoes. % Tl = Tl of tissue (ms) % T2 = T2 of tissue (ms) % tau = time between 90-180 (ms) jay = sqrt(-l); % Preallocate the signal intensity vector decay = zerosG, num_echoes) ; % Preallocate the array of num^echoes substates [F F* Z Z*] and initialize B = zeros(num_echoes*4, 1); B(l) = 0 + jay*l; % Precalculate rotation matrix Tone for one substate block Tone = [[ (cos(alpha/2))"2, (sin(alpha/2))"2, -jay*sin(alpha), 0]; ... [ (sin(alpha/2))~2, (cos(alpha/2))"2, 0, jay*sin(alpha)]; ... [ (-jay)*l/2*sin(alpha), (jay)*l/2*sin(alpha), cos(alpha), 0] ; ... [ (-jay)*l/2*sin(alpha), (jay)*l/2*sin(alpha), 0, cos(alpha)]]; % Calculate the sparse rotation matrix of "num-echoes" blocks Tp = kron(sparse(diag(ones(l.nunuechoes))), create.Tp(alpha)); % Apply initial T2 decay and Tl regrowth (between 90-180 pulses). B = t2relax(B, tau, T2); B = tlrelax(B, tau, Tl); for echo = 1:nunuechoes % Apply rotation matrix. % Each block of the rotation matrix acts on [F(k), F*(k), Z(k), Z*(k)] B = (Tp * B); % Save echo intensity, this is Fl* decay(echo) = abs(B(2)) * exp(-tau/T2); % Apply transverse decay for time 2*tau. B = t2relax(B, 2*tau, T2); % Apply longitudinal regrowth for time 2*tau. B = tlrelax(B, 2*tau, Tl); % Let the spins dephase for time 2*tau B = evolution(B); end 100 CHAPTER 4. NON-180 REFOCUSING PULSES 4.1. INTRODUCTION Algorithm 4.2 Ti relaxation function [B] = tlrelax(B, tau, Tl) % The Z and Z* substates will be modified indicies = [ 3:4:length(B) 4:4:length(B) ]; % Apply the Tl relaxation to the Z and Z* substates B(indicies) = B(indicies) * exp(-tau/Tl); Algorithm 4.3 T2 relaxation function [B] = t2relax(B, tau, T2) % The F and F* substates will be modified indicies = [ 1:4:length(B) 2:4:length(B) ]; % Apply the T2 relaxation to the Z and Z* substates B(indicies) = B(indicies) * exp(-tau/T2); Algorithm 4.4 Phase evolution of magnetization function [Bnew] = evolution(B) % This function is a book-keeping method to track the phases of the spins % as they are flipped and as they progress through the inter-pulse time. Bnew = B; % Number of elements in B L = length(B); % All F(n) go to F(n+1). Bnew([ 5:4:L ]) = B([ l:4:L-4 ]); % F*(l) goes to F(l). Bnew(l) = B(2); % All F*(n) go to F*(n-1). Bnew([ 2:4:L-4 ]) = B([ 6:4:L ]); 101 CHAPTER 4. NON-180 REFOCUSING PULSES 4.2. MONO-EXPONENTIAL FIT 4.1.4 Overview of Work I created a method to fit decay curves acquired from a multi-echo pulse sequence with non-180° refocusing pulses. This method was implemented in two algorithms: first, a non-linear curve fitting algorithm to estimate mono-exponential decay curve parameters (Section 4.2) and second, a non-negative least squares algorithm to estimate multi-exponential decay curve parameters (Section 4.3). A set of simulations were used to determine the accuracy and precision of the parameter estimation using each fitting algorithm. Scans of phantoms, from the MESEC sequence with a train of a pulses, were used to verify the parameter estimation. A set of in vivo scans were acquired from volunteers and myelin water maps calculated from the MESEC sequence were compared to myelin water maps calculated from the MESE0. Several extensions to the work will be shown, including: slab selected, multi-slice acquisition (Section 4.4), the effect of the non-90° excitation pulse (Section 4.5), and trains of non-constant refocusing pulses (Section 4.6). 4.2 Mono-Exponential T2 Decay Curve Fit An algorithm to fit mono-exponential decay curves is described. Noise simulations were used to assess the accuracy and precision of the fit algorithm. Finally, a set of nickel/agarose phantoms, of known T2 values, were scanned with different sets of refocusing pulses. The estimated T2 was compared to the T2 estimated from the MESE0 sequence and the refocusing pulse were compared to the prescribed. 4.2.1 Mono-Exponential Fit Algorithm Each decay curve was fit to a function of proton density p, spin-spin relaxation time T2, baseline offset d, and refocusing pulse flip angle a. A non- linear least squares algorithm (lsqcurvef it in Matlab [The Mathworks, Natick, MA]) minimized the misfit x2 = YliLiiVi ~ Vi)2la1 where y and a were the mean and standard deviation of the measured data and y was the decay curve estimated by the algorithm. Initial conditions for the curve fitting were: p = 1200, T2 = 100 ms, a — 160°, d = 5. The reconstructed decay curve was calculated at each step of the least squares algorithm by calculating the primary, stimulated, and indirect echoes as defined in Hennig [78]. TR was considered infinite for all fitting. 102 CHAPTER 4. NON-180 REFOCUSING PULSES 4.2. MONO-EXPONENTIAL FIT 4.2.2 Monoexponential Fit — Noise Simulation A set of noise simulations was used to assess the decay curve parameter accuracy and precision from mono-exponential decay curves based on non-1800 refocusing pulses. The mean and standard deviation were calculated for each parameter over the noise realizations. Methods A noiseless decay curve was created based on the model: y = pexp(-t/T2) (4.17) where p = 1000 and T2 = 20,80,140 ms. The decay curve was created based on flip angles of a = 90°, 100°,..., 180° and echo times t = 10, 20,... , 320 ms. One thousand realizations of quadrature noise (SNR =200) were added. TR was set to infinity. The T2 decay curve parameters used in the simulation are similar to values for brain tissue. The mean and standard deviation of the estimated parameters were calculated and compared to the true values. Results and Discussion Each parameter estimated from the decay curve, for T2 = 20,80,140 ms was within 5% of the true parameter. Table 4.1 shows the results for T2 = 80 ms. The parameters estimated from the decay curve were within 2% of the true parameter for T2 = 80,140 ms. The increased accuracy, as a function of T2, was expected as more echoes were able to contribute information to the solution as T2 increased. The variability in the estimated parameters was higher for the lower flip angles (for example, variability in T2 in Table 4.1). To determine if the increased variability is related to the decrease in the flip angle, I re-ran the simulations and scaled the 1,000 realizations of noise by, sin2 (a/2), on the relative attenuation due to the flip angle. For example, I multiplied the random noise by sin2 (a/2) = 1.0 for a = 180° and by sin2 (a/2) = 0.5 for a = 90°. In this way the noise decay curves estimated with a = 90° had the same SNR as the curves simulated with a = 180°. Then any discrepancies in the standard deviations of the parameters, were not be due to SNR issues. 103 CHAPTER 4. NON-180 REFOCUSING PULSES 4.2. MONO-EXPONENTIAL FIT "true (°) A(°) P T2 ( ms) x1 090 89.836 (0.964) 1001.220 (12.455) 79 296 (1.266) 28 413 (8 050) 100 99.848 (1.015) 1000.919 (11.116) 79 407 (1.139) 28 391 (8 017) 110 109.849 (1.117) 1000.774 (10.058) 79 492 (1.035) 28 285 (8 025) 120 119.851 (1.148) 1000.622 (8.817) 79 557 (0.939) 28 333 (8 070) 130 129.839 (1.269) 1000.551 (8.040) 79 594 (0.867) 28 270 (8 035) 140 139.848 (1.450) 1000.370 (7.292) 79 630 (0.799) 28 185 (8 058) 150 149.841 (1.669) 1000.269 (6.490) 79 650 (0.735) 28 137 (7 992) 160 159.851 (1.984) 1000.144 (5.717) 79 669 (0.686) 28 141 (8 Oil) 170 169.992 (2.726) 1000.027 (4.823) 79 681 (0.638) 28 177 (8 049) 180 177.474 (2.925) 1000.391 (4.179) 79 636 (0.605) 28 642 (8 030) Table 4.1: The estimated p and T2 (standard deviation in brackets) given the true p = 1000 and true T2 = 80 ms. "true (°) a{°) P T2 ( ms) X1 090 89.900 (0.464) 1000.469 (5.860) 79 532 (0.658) 7.132 (1.753) 100 99.902 (0.587) 1000.287 (6.211) 79 588 (0.709) 9.794 (2.384) 110 109.912 (0.670) 1000.094 (6.026) 79 638 (0.711) 12.822 (3.137) 120 119.932 (0.887) 999.904 (6.408) 79 678 (0.755) 15.928 (3.916) 130 129.982 (1.026) 999.612 (6.367) 79 707 (0.756) 19.117 (4.719) 140 139.988 (1.245) 999.597 (6.037) 79 704 (0.728) 22.015 (5.450) 150 149.964 (1.592) 999.713 (5.996) 79 689 (0.719) 24.446 (6.032) . 160 160.159 (1.969) 999.381 (5.327) 79 711 (0.696) 26.423 (6.618) 170 170.412 (2.770) 999.490 (4.645) 79 693 (0.672) 27.776 (6.919) 180 177.861 (2.579) 1000.032 (4.348) 79 634 (0.661) 28.757 (6.986) Table 4.2: The estimated p and T2 from a simulation with SNR scaled as a factor of sin2 (a/2) to determine whether the increased variability in the estimated parameters was due to the decrease in a. 104 CHAPTER 4. NON-180 REFOCUSING PULSES 4.2. MONO-EXPONENTIAL FIT The estimated parameters, derived from the variable noise decay curves, are shown in Ta ble 4.2. The standard deviation of p and T2 remain relatively constant as a decreased to 90°. The constant variability was true for T2 = 140 ms as well. When T2 = 20 ms, the variability was still a function of a, though decreased by approximately half when the SNR was scaled by the effect of the refocusing pulse. Conclusions Mono-exponential T2 parameters were accurately and precisely estimated from decay curves sim ulated from non-1800 refocusing pulses. The accuracy and precision were found to depend on a, and the short T2 was affected more by the decreased flip angle than the longer T2 components. 4.2.3 Phantom Verification A set of mono-exponential phantoms were used to confirm the accuracy of T2 from a set of scans. A set of six nickel/agarose phantoms of known T2 were scanned (Section 1.8) with a 16 echo pulse sequence with prescribed refocusing pulse flip angles. Methods The data was acquired with the single slice, 16 echo MESEC pulse sequence with parameters: echo spacing = 11.1 ms, TR = 3000 ms and slice selective refocusing pulses. Regions were drawn on the TE = 11.1 ms image for each of the phantoms. The signal was averaged, within the region, for each echo to create one decay curve for each phantom. The decay curves were fit using the mono-exponential algorithm (Section 4.2.1) in three ways: flip angle set to the prescribed flip angle ap, flip angle set to aiso = 180°, and flip angle estimated to minimize misfit (abest)-Results and Discussion The measured data and curve fits are shown in Figures 4.2, 4.3, 4.4, and 4.5, for one phantom (T2 = 90.8 ms). Curves fit with a as a parameter in the fitting algorithm were better representative of the subtle changes in the curve as seen in the decreased residuals (small panels of Figures 4.2, 4.3, 4.4, and 4.5). In comparison, the curves fit with a fixed to 180° had large residuals that 105 CHAPTER 4. NON-180 REFOCUSING PULSES 4.2. MONO-EXPONENTIAL FIT alternated positive and negative for successive echoes. The curve fitting algorithm worked well over a range of refocusing pulse flip angles and T2 (not shown). The estimated parameters of the mono-exponential fits to phantoms are shown in Table 4.3. The estimated flip angle a was always lower than the prescribed flip angle and was significantly lower when ap = 180°. The T2 was typically within ±10% of the T2 estimated from the 48 echo MESE0 data. OLp "best P180 PP Pbest T2,180 t2,p t2,best Xl80 4 Xbest 90 85 498.2 792.7 842.7 136.4 88.5 , 83.4 .141.1 4.0 1.9 110 103 596.8 777.2 816.5 123.7 96.4 87.8 81.8 3.3 0.6 130 123 687.3 783.1 813.0 110.4 97.4 92.4 49.9 3.4 0.9 150 137 750.3 786.2 826.2 101.7 97.2 92.5 23.7 6.5 0.8 170 151 784.3 788.5 819.2 98.3 97.7 94.0 7.5 5.7 1.3 180 156 800.2 800.2 819.9 97.0 97.0 92.0 5.4 5.4 1.0 Table 4.3: Phantom decay curve parameters estimated from mono-exponential fits to decay curves collected using a train of 16 ap refocusing pulses. Parameters with a subscript of p were estimated with a fixed to ap. Parameters with a subscript of best were estimated allowing a to vary. Results are from one phantom (T2 = 90.8 ms). Flip angles are in degrees and T2 in milliseconds. The measured refocusing pulse flip angles were always lower than expected. This is likely due to the shape of the refocusing pulse slice profile. A good slice selection profile requires a high power refocusing pulse, but the power of the refocusing pulse is limited as the pulse sequences used to acquire multi-echo data requires a train of RF pulses. Therefore, the slice profile will have a transition from no refocus to a = 180°. This transition zone will.affect the integrated flip angle. The mono-exponential curve fit was accurate at reproducing T2 and a. To assess the effect of initial conditions on the final solution, a set of six other initial conditions were used one with each initial parameter above and then below the original initial condition. For every initial condition, the T2, p and a from each phantom decay curve were estimated to be the same as the result from the original initial condition. Therefore, the curve fits were reliable and not dependent on initial conditions. 106 CHAPTER 4. NON-180 REFOCUSING PULSES 4.2. MONO-EXPONENTIAL FIT ro TE (ms) Figure 4.2:, Decay curves calculated as the mean signal intensity over a region of one of the phan toms (error bar is one standard deviation). The large panel shows the measured data (circles), the mono-exponential least squares fit based on the prescribed refocusing pulse flip angle ap — 180° (solid line), the curve where the flip angle was estimated as a parameter of the fitting method (dashed line) and the fit in which the flip angle was assumed to be 180° (dash-dotted line). The small panel shows the residuals of the corresponding fit. T2 = 90.8 ms for the phantom. 107 CHAPTER 4. NON-180 REFOCUSING PULSES 4.2. MONO-EXPONENTIAL FIT 700 TE (ms) Figure 4.3: Decay curves calculated as the mean signal intensity over a region of one of the phan toms (error bar is one standard deviation). The large panel shows the measured data (circles), the mono-exponential least squares fit based on the prescribed refocusing pulse flip angle av = 150° (solid line), the curve where the flip angle was estimated as a parameter of the fitting method (dashed line) and the fit in which the flip angle was assumed to be 180° (dash-dotted line). The small panel shows the residuals of the corresponding fit. T2 = 90.8 ms for the phantom. 108 CHAPTER 4. NON-180 REFOCUSING PULSES 4.2. MONO-EXPONENTIAL FIT Figure 4.4: Decay curves calculated as the mean signal intensity over a region of one of the phan toms (error bar is one standard deviation). The large panel shows the measured data (circles), the mono-exponential least squares fit based on the prescribed refocusing pulse flip angle ap = 110° (solid line), the curve where the flip angle was estimated as a parameter of the fitting method (dashed line) and the fit in which the flip angle was assumed to be 180° (dash-dotted line). The small panel shows the residuals of the corresponding fit. T2 = 90.8 ms for the phantom. 109 CHAPTER 4. NON-180 REFOCUSING PULSES 4.2. MONO-EXPONENTIAL FIT n 3 TJ W O tr 200 -200 1 1 1 0 50 100 150 200 Figure 4.5: Decay curves calculated as the mean signal intensity over a region of one ofthe phan toms (error bar is one standard deviation). The large panel shows the measured data (circles), the mono-exponential least squares fit based on the prescribed refocusing pulse flip angle ap — 90° (solid line), the curve where the flip angle was estimated as a parameter of the fitting method (dashed line) and the fit in which the flip angle was assumed to be 180° (dash-dotted line). The small panel shows the residuals of the corresponding fit. T2 = 90.8 ms for the phantom. 110 CHAPTER 4. NON-180 REFOCUSING PULSES 4.3. MULTI-EXPONENTIAL FIT Conclusions The decay curves from the mono-exponential phantoms were reliably fit using a mono-exponential curve fit algorithm designed to calculate the refocusing pulse flip angle during the fit. The decay curve fits that estimated the refocusing pulse flip angle were better representative of the acquired decay curve as seen by the decrease in the x2 misfit. 4.3 Multi-Exponential T2 Decay Curve Fit Similar to the previous section, an algorithm is defined for fitting multi-exponential decay curves collected from a pulse sequence with a train of a° refocusing pulses. The algorithm will be validated based on a set of simulated and in vivo multi-exponential decay curves. 4.3.1 Multi-Exponential Fit Algorithm 4.3.1.1 T2 Decay Curve Fit Using NNLS Multi-exponential fits of the in vivo data were done with the non-negative least squares (NNLS) algorithm [18, 39, 40]. NNLS solves ^4s = y where Aij = exp (-TEi/T2j), y is the measured decay curve and Sj are the unknown amplitudes that correspond to T2j. The solution minimizes the x2 misfit. In the standard NNLS algorithm, each column Aj, of the matrix A, is a mono-exponential decay curve with T2 = T2j. To account for a, I created a matrix Aa where each column was a decay curve, calculated using Hennig's method (Section 4.1.2), with T2 = T2 j and flip angle a. Therefore, the objective was to find the minimum x2 ia) misfit of ^4Qs = y over the range of a. To do this, I created a new algorithm, NNLSa, to find the minimum x2 (a) by solving Aas = y for a between 40° and 180°. The multi-exponential fit algorithm I introduced, NNLS", minimized x2 as a function of a. Figure 4.6 shows three example curves of x2 as a function of a for ap = 90° (solid line), ap = 150° (dashed line) and ap = 180° (dash-dotted line). Each of the curves has a distinct minimum at the true a. Therefore, the goal is to find a fast, robust method to find the refocusing pulse flip angle that minimizes the x2 misfit. The x2 misfit increases quickly as the a° refocusing pulse flip angle deviates from the true 111 CHAPTER 4. NON-180 REFOCUSING PULSES 4.3. MULTI-EXPONENTIAL FIT 200 180 160 140 £ 120 ro co 100 i O 80 60 40 20 !—' / 1 1 r 1 1 1 / True a = 90 - - • Truea=150 ' -- Truea=180 1 :/:*': 1 / 1 ; / 1 i 1 : / i 1 1 : : / i . i I : : / I \ \ / ' > \: ;/ \ • I : / i : \ \ 7 *: V • / • / / / ' I t \ / \ / \ \ \ ' \ * / ^ / * y 1 1 ' 1 ' 1 1 1 60 80 100 120 140 Flip Angle 160 180 Figure 4.6: The x2 misfit plotted as a function of refocusing pulse flip angle a for a true flip angle of a = 90° (solid line), a = 150° (dashed line), and a = 180° (dash-dotted line). The goal of the outer loop of the NNLS" algorithm is to find the minimum x2 over a. 112 CHAPTER 4. NON-180 REFOCUSING PULSES 4.3. MULTI-EXPONENTIAL FIT refocusing pulse flip angle (as seen in Figure 4.6). The goal is to create a fast and reliable algorithm to determine the refocusing pulse flip angle a with the minimum x2 • Two types of algorithms were assessed for speed and robustness. The first type of algorithm samples several points along the x2/a curve then fits the points and determines the minimum x2 from the curve fit. Two curve fit algorithms were tried, a quadratic fit and a spline fit. The second type of algorithm is the standard minimization technique where the minimum x2 is found by minimizing the function X2 {a) over a. The second technique was inherently slower, but more robust. 4.3.1.2 Analysis ofthe Aa Matrix Introduction The decay curves collected with a < 180° have more structure than decay curves with a = 180°. Therefore, it may be the accuracy of quantitative T2 measurements is higher when a < 180° than when a = 180°. Decay curves near 150° tend to have more structure and therefore the columns of A may be more linearly independent. Based on that theory, I defined an Aa matrix which was calculated for a — 1°, 2°,..., 270°. Two measures were calculated for each matrix 1) the rank and 2) the inverse of the condition number. The measures were calculated in Matlab (The Mathworks, Natick Mass) by the built-in ; functions rank and svd, respectively. The rank provides an estimate of the number of linearly , independent rows or columns. The condition number is defined as the ratio of the largest singular value to the smallest, but for this work I used the ratio of the smallest singular value to the largest. Two sets of Aa matrices were calculated. The first set was based on 32 echoes, TE = 10 ms, Ti = 800 ms. The second set was based on 200 echoes, TE = 5 ms, Ti = 800 ms. The Ti was chosen to be 800 ms as that is approximately the Ti of white and gray mater. Results 32 Echoes The rank (Figure 4.7) was found to be a maximum (rank of 31) between 136° < a < 158° and was near the lowest (rank of 18) at a = 180°. Similar results were found for the inverse of the condition number (Figure 4.7). The largest 113 CHAPTER 4. NON-180 REFOCUSING PULSES 4.3. MULTI-EXPONENTIAL FIT Figure 4.7: The rank (left) and inverse of the condition number (right) of• Aa based on a 32 echo decay with TE = 10 ms and TR = 800 ms. ratio (smallest to largest singular values) was 105.2686 x 10~17 at a = 158° and the lowest was 0.4996 x 10~17 at a = 180°. 200 Echoes Similar results were found with 200 echoes. The rank (Figure 4.8) was found to be a maximum (rank of 50) at a = 156° and was the lowest (rank of 24) at a = 180°. Similar results were found for the inverse of the condition number (Figure 4.8). The largest ratio (smallest to largest singular values) was 0.8023 x 10-17 at a = 150° and the lowest was 0.0021 x 10~17 at a = 180°. The small differences in the inverse of the condition number on opposite sides of a = 180° was found to be due to very small differences in the Aa matrices. For example, the matrix with the largest inverse condition number was A208, but the difference between The mono-exponential decay curves in A208 and A152 was less than IO-13. Conclusions The matrix A180 has the lowest rank and the highest condition number and therefore is about the worst matrix with which to calculate T2 distribution. By measuring the decay curve with a refocusing pulse flip angle of a = 159° a higher rank (almost double) and a lower condition 114 CHAPTER 4. NON-180 REFOCUSING PULSES 4.3. MULTI-EXPONENTIAL FIT 300 Flip Angle (degrees) Flip Angle (degrees) Figure 4.8: The rank (left) and inverse of the condition number (right) of Aa based on a 200 echo decay with TE = 5 ms and TR = 800 ms. number (two orders of magnitude) should allow for better T2 quantification. Another benefit of a multi-echo pulse sequence using 159° non-composite pulses is that the SAR would be lower than that of a sequence using 180° pulses by a factor of (180/159)2 = 1.28. Even though there is a large relative change in the inverse of the condition number from a = 159° to a = 180°, the magnitude of the inverse of the condition number for both flip angles is very large. Therefore, the theoretical gain in the matrix based on a = 159° would likely not be realized given the SNR of the data from typical in vivo MRI data. 4.3.2 Bi-Exponential Parameters From Incorrect a To determine the miscalculation of the bi-exponential parameters that result when an incorrect flip angle is assumed. A bi-exponential model was used with a short T2 compartment of p s = 200 and T2S = 20 ms and p m = 800 and T2m = 80 ms. A true decay curve was created using the bi-exponential model parameters and refocusing pulse flip angles of atrUe = 90°, 100°,..., 180°. For each true decay curve one thousand realizations of SNR — 200 quadrature noise (Gaussian on two channels) were added. The non-negative least squares algorithm was applied to each decay curve, assuming the refocusing pulse was a = 180°, and the decay curve parameters were calculated. 115 CHAPTER 4. NON-180 REFOCUSING PULSES 4.3. MULTI-EXPONENTIAL FIT The myelin water fraction, proton densities and T2s had less than 10% error for a > 170° (Table 4.4). All parameters were incorrectly determined when a < 170°. Therefore, accurate and precise quantification of T2 parameters requires the flip angle be calculated as well. Q'true fs fm T m 2 x1 90 0.000 (0 000) 0.000 (0.000) 1.000 (0 000) 107.707 (0.956) 1123.617 100 0.000 (0 000) 0.000 (0.000) 1.000 (0 000) 98.013 (1.204) 955.397 110 0.000 (0 000) 0.000 (0.000) 1.000 (0 002) 89.886 (1.212) 767.491 120 0.000 (0 000) 0.000 (0.000) 0.984 (0 029) 82.727 (1.783) 579.248 130 0.003 (0 029) 0.716 (5.805) 0.919 (0 083) 75.387 (4.323) 405.741 140 0.150 (0 147) 25.070 (22.278) 0.838 (0 142) 82.402 (9.003) 256.958 . 150 0.255 (0 076) 34.956 (7.518) 0.744 (0 074) 86.947 (4.419) 142.229 160 0.219 (0 048) 26.921 (5.214) 0.779 (0 047) 82.931 (3.244) 66.923 170 0.209 (0 048) 21.753 (5.071) 0.790 (0 047) 80.270 (3.217) 32.087 180 0.204 (0 044) 19.804 (4.620) 0.794 (0 044) 79.259 (2.962) 26.948 Table 4.4: Measured parameters of a bi-exponential T2 decay curve when the nominal flip angle is assumed to be a = 180°, but isn't. The true parameters are: fs = 0.20, dns = 200, iy = 20, fm = 0.80, dnm = 800, T2m = 80. 4.3.3 NNLSQ Multi-Exponential Simulations A set of noise simulations were used to assess the accuracy and precision of the parameters of the T2 distribution and refocusing pulse flip angle. A bi-exponential T2 distribution as in Section 4.3.2 was used. The true parameters were fs = 0.20, T2S = 20 ms, fm = 0.80, and T2m = 80 ms. The mean and standard deviation of the estimated parameters are shown in Table 4.5. As in the mono-exponential results, the accuracy of the parameters was within 10% for a > 100°. The standard deviation increased as the refocusing pulse decreased. As in the mono-exponential case, the simulations were repeated with the noise scaled by 1/sin2 (a/2). When the noise was scaled, the standard deviations of the parameters were similar as a function of a. Noise simulations assessed the variability of the decay curve parameters for a bi-exponential distribution simulated with refocusing pulses from 90° to 180°. Table 4.5 shows the means and standard deviations of the parameters estimated by the NNLSQ algorithm. All estimated means were within 10% of truth for a > 100°. In the multi-exponential noise simulations, the mean a was 2.5° lower than the true a and the variability was higher compared to the variability of the other true a. This was likely due 116 CHAPTER 4. NON-180 REFOCUSING PULSES 4.3. MULTI-EXPONENTIAL FIT to noise in the decay curve causing noise in the curve near 180° and the minimization of the outer loop settling at slightly less than a — 180°. The NNLS° algorithm was robust in the presence of noise and was able to find the best refocusing pulse flip angle. ^true a fs T2), fm T2,m 90 89.87 (1.44) 0 20 (0 13) 15.90 (10.63) 0 79 (0 17) 73.33 (15.66) 100 99.81 (1.47) 0 20 (0 12) 16.65 (9.67) 0 79 (0 15) 74.50 (13.89) 110 109.83 (1.57) 0 21 (0 11) 17.75 (9.17) 0 78 (0 12) 76.12 (12.15) 120 119.81 (1,63) 0 20 (0 10) 18.06 (8.20) 0 79 (0 10) 76.55 (9.98) 130 129.80 (1.68) 0 21 (0 09) 18.32 (7.70) 0 78 (0 10) 76.77 (9.78) 140 139.72 (1.72) 0 20 (0 08) 18.33 (6.43) 0 79 (0 08) 77.06 (7.98) 150 149.88 (1.96) 0 20 (0 07) 18.61 (6.11) 0 79 (0 07) 77.12 (7.44) 160 159.80 (2T0) 0 20 (0 07) 18.29 (6.29) 0 79 (0 07) 77.16 (7.03) 170 170.01 (2.98) 0 20 (0 07) 17.81 (6.87) 0 79 (0 07) 77.04 (7.21) 180 176.76 (2.97) 0 21 (0 07) 17.97 (6.62) 0 79 (0 07) 77.61 (6.78) Table 4.5: Noise simulation where the decay curve parameters and refocusing pulse flip angle were estimated from 1,000 realizations of noise for a = 90°, 100°,..., 180°. The true bi-exponential distribution was /s = 0.2, T2]S = 20 ms, fm = 0.8 and T2,m = 80 ms. Flip angles are in degrees and T2 in milliseconds. 4.3.4 Phantoms The multi-exponential algorithm, NNLS", to calculate the multi-component T2 was validated based on a set of MRI scans of the nickel/agarose phantoms. The phantoms were scanned with the MESEC sequence, from which, the T2 were compared to T2 calculated from the MESEG sequence. The refocusing pulse flip angles were calculated from the MESEC sequence, as well, and compared to flip angles calculated from a technique by Stollberger [81]. The Ti relaxation times were estimated from a pulse-saturation experiment for reference. 4.3.4.1 Methods Data Acquisition The T2 measurements were done with a 32 echo MESE0 pulse sequence with TE = 10, 20,... , 320 ms. Other parameters were TR = 3000 ms, 2 averages. The T2 measurements were done with a 32 echo MESEC pulse sequence with and echo spacing of 12.1 ms. Other parameters were TR = 3000 ms, 2 averages. The Bi measurements were calculated from a pair of spin echo images [81], the 117 CHAPTER 4. NON-180 REFOCUSING PULSES 4.3. MULTI-EXPONENTIAL FIT first image was a standard spin-echo image, the second had an excitation pulse of 180°. Other parameters included TE = 11.1 ms, TR = 3000 ms, and no signal averaging. Data Analysis For each multi-echo data acquisition the decay curve was calculated as the mean signal intensity over the phantom and was fit using NNLS. The Bi, from the pair of spin-echo images, was calculated as [81]: where Ia was the standard spin-echo image with a prescribed excitation pulse of a = 90°, fya was the spin-echo image with a prescribed excitation pulse of 180°, and 0:5 was the spatially varying flip angle calculated by the Stollberger method. The excitation pulses will vary spatially and therefore, the prescribed a = 90° will not be 90° throughout the slice. SI(TRi) = p[l -exp(-TR/Ti)]. 4.3.4.2 Results T2 Calculation The T2 relaxation times, shown in Table 4.6, were calculated from the MESE0 sequence (for reference) and from the MESEC sequence. The T2 relaxation times calculated from the decay curve collected by the MESEC sequence were within 8% of the relaxation times calculated from the decay curve collected by the optimized MESE0 sequence. The Ti relaxation times, calculated from the pulse-saturation sequence, are shown in Table 4.8, for reference. (4.18) For each phantom, the mean signal intensity over the phantom was fit to the equation 118 CHAPTER 4. NON-180 REFOCUSING PULSES 4.3. MULTI-EXPONENTIAL FIT Phantom MESEQ T2 (ms) MESEC T2 (ms) bottomJeft 24.3 22.9 bottom-middle 86.4 80.7 bottom_right 26.1 24.2 top-left 316.6 304.6 top-middle 99.3 91.7 topjight 110.6 102.8 Table 4.6: T2 time calculated from the MESE0 sequence and the MESEC (with the flip angle calculated). The refocusing pulse flip angle was calculated by the Stollberger method (as) and from the decay curves acquired with the MESEC sequence (a) and are shown in Table 4.7. The refocusing pulses calculated from the decay curve were within 6% of the refocusing pulses calculated by the Stollberger method. Phantom a (°) as (°) bottom-left 163.0 165.7 bottom_middle 159.0 166.3 bottom_right 163.0 167.7 topJeft 156.0 164.8 top-middle 158.0 165.1 top_right 160.0 166.9 Table 4.7: The refocusing pulse flip angle was calculated by the Stollberger method (as) and from the decay curves acquired with the MESEC sequence (a). Phantom Ti (ms) bottom-left 343.4 bottom_middle 352.9 bottom-right 690.7 top-left 730.7 top-middle 712.0 top-right 1843.9 Table 4.8: Estimated Ti from curve fits of five points (TR = 100, 300, 800, 2000, 3000 ms) from a pulse saturation experiment. 119 CHAPTER 4. NON-180 REFOCUSING PULSES 4.3. MULTI-EXPONENTIAL FIT 4.3.4.3 Conclusions The nickel/agarose phantoms were scanned using the MESEC sequence from which the T2 and refocusing pulse flip angle a were estimated. The T2 estimated from each phantom was within 10% of the estimated T2 from the MESE o sequence. The estimated refocusing pulse flip angle was also within 10% of the flip angle estimated by the Stollberger method. 4.3.5 In Vivo Multi-Exponential Fits Two sets of in vivo scans were acquired, of a volunteer, with different setting of the refocusing pulse flip angle. For one scan, the flip angle was set to anom = 180° and for the second scan, the flip angle was set to anom = 150°. From these two scans, the myelin water map was calculated, as was a map of the refocusing pulse flip angle and the mean T2. 4.3.5.1 Methods Two sets of scans were done on a normal volunteer using the 16-echo MESEC with TR — 3 s, TE = 11.096 ms, 256 x 128 freqxphase, 24 cm FOV and 4 averages. The first scan was acquired as per normal with anom = 180°. The second scan was acquired with anom = 150°, by scaling the amplitude of the refocusing pulses by a factor of 150/180. Due to the long scan time of the two sequences, a MESE0 scan was not acquired. The T2 distribution was estimated voxel-by-voxel in each scan as well as the refocusing pulse flip angle using the algorithm described in Section 4.3.1. Maps of the myelin water fraction, refocusing pulse flip angle and mean T2 calculated from each scan were compared visually. For the first analysis, the T2 information was calculated along with the flip angle at every voxel ("Variable Flip Angle") in the image and for the second analysis the flip angle was assumed a = 180° at every voxel ("Fixed Flip Angle"). 4.3.5.2 Results Variable Flip Angle Maps of the myelin water were estimated from the multi-echo data with anom = 180° and with anom = 150° and are shown in Figure 4.9. The map of the myelin water was qualitatively very similar to other myelin water maps acquired with the optimized MESE0. The map from the 120 CHAPTER 4. NON-180 REFOCUSING PULSES 4.3. MULTI-EXPONENTIAL FIT Figure 4.9: The myelin water maps calculated from the pulse sequence with a refocusing pulse flip angle prescribed as 180° (left) and 150° (right). c*nom = 180° (Figure 4.9, left) was brighter in the top part of the image compared to the bottom part of the image. This shading was not seen on the myelin water map from anom = 150° (Figure 4.9, right). The histogram of the myelin water fraction, within the brain, was narrower for anom = 150° than for anom = 180°. As well, there was more myelin water structure visible in the middle portion of the myelin water map acquired with an0m = 150°. The bright phantoms on the right side of each myelin water map are the short T2 nickel/agarose phantoms. Region MWF (anom = 150°) MWF (anom = 180°) caudate (left) 3.6 (1.7) 0.8 (1.5) caudate (right) 2.4 (2.2) 0.1 (0.2) putamen (left) 3.2 (2.7) 1.6 (1.7) putamen (right) 3.9 (2.3) 0.1 (0.0) frontal (left) 16.4 (2.3) 20.9 (2.0) frontal (right) 15.3 (1.9) 17.8 (4.1) genu 18.1 (5.9) 24.6 (4.4) ic (left) 20.8 (1.5) 24.4 (4.5) ic (right) 16.9 (2.5) 22.4 (1.7) post (left) 13.0 (3.6) 13.3 (2.6) post (right) 10.8 (4.3) 8.3 (2.0) Table 4.9: Percent myelin water estimated from scans with anom = 150° and anom = 180°. 121 CHAPTER 4. NON-180 REFOCUSING PULSES 4.3. MULTI-EXPONENTIAL FIT The maps of the refocusing pulse flip angles were calculated voxel-by-voxel from the two scans, one with anom = 180° and the other with anom = 150°, and are the top pair of images in Figure 4.10. The mean refocusing pulse flip angle from MESEC with «nom = 150° was 138.2° and the mean refocusing pulse flip angle with anom = 180° was 156.3°. The bottom pair of images in Figure 4.10 are the refocusing pulse flip angle maps normalized by the mean refocusing pulse flip angle. Figure 4.10: The refocusing pulse flip angle maps calculated from the pulse sequence with a refocusing pulse flip angle prescribed as 180° (left) and 150° (right). The top plot is in degrees and the bottom plot was normalized by the mean flip angle over the whole map. 122 CHAPTER 4. NON-180 REFOCUSING PULSES 4.3. MULTI-EXPONENTIAL FIT Region & (anom = 150°) a (ttnom = 180°) caudate (left) 140.0 (1.9) 157.5 (1.6) caudate (right) 142.3 (0.9) 158.6 (1.1) putamen (left) 139.9 (1.3) 157.2 (1.0) putamen (right) 140.8 (1.1) 156.6 (2.4) frontal (left) 139.8 (1.2) 159.3 (1.5) frontal (right) 141.5 (1.1) 161.4 (1.5) genu 142.2 (4.4) 163.6 (3.9) ic (left) 143.4 (1.4) 161.1 (1.3) ic (right) 143.4 (2.1) 158.1 (2.7) post (left) 139.1 (2.0) 159.1 (2.0) post (right) 138.4 (1.8) 158.7 (1-7) Table 4.10: Refocusing pulse flip angle calculated over the regions from the scans with anom = 150° and anom = 180°. 120 100 Figure 4.11: The arithmetic mean T2 maps calculated from the pulse sequence with a refocusing pulse flip angle prescribed as 180° (left) and 150° (right). 123 CHAPTER 4. NON-180 REFOCUSING PULSES 4.3. MULTI-EXPONENTIAL FIT The geometric mean T2 of the medium T2 component estimated from MESEC with ttn0m = 180° and aMm = 150° are the left and right images, respectively, in Figure 4.11. The mean T2 of the medium component was calculated from several regions within the brain and is shown in Table 4.11. The T2m from the two scans were within 4%. Region T2m (ms) T2m (ms) caudate (left) 81.1 (2 1) 79.0 (1.6) caudate (right) 80.5 (1 5) 79.4 (1.1) putamen (left) 71.6 (1 6) 70.3 (0.6) putamen (right) 71.6 (1 7) 70.5 (1.4) frontal (left) 76.9 (1 3) 77.5 (2.4) frontal (right) 76.8 (1 1) 75.4 (1.7) genu 74.6 (5 4) .74.1 (2.0) ic (left) 77.3 (0 8) 79.3 (5.2) ic (right) 76.9 (1 0) 78.4 (4.2) post (left) 78.4 (1 9) 81.3 (4.1) post (right) 76.9 (1 2) 77.7 (1.2) Table 4.11: Geometric mean T2 of the medium component from select regions throughout the image. The x2) calculated during the estimation of the T2 distribution, is shown in Figure 4.12. The regions of high \2•, and therefore, poor curve fit, are primarily in the cerebrospinal fluid (CSF). There are two possible reasons. First, the CSF is flowing and therefore may cause artifacts on the decay curves that can not be accounted for by echo formation. Second, the T2 of CSF is very long, but for this data only 16 echoes were collected with an echo spacing of 11.1 ms, therefore, the decay curve may not be sufficiently sampled to accurately determine T2 distribution of CSF. 124 CHAPTER 4. NON-180 REFOCUSING PULSES 4.3. MULTI-EXPONENTIAL FIT Figure 4.12: The x2 maps calculated from the pulse sequence with a refocusing pulse flip angle prescribed as 180° (left) and 150° (right). Fixed Flip Angle For comparison, the T2 distribution was calculated voxel-by-voxel from the two scans assuming anom = 180°. Based on the results in Section 4.3.2, I expect the myelin water fraction to decrease in regions where a < 180°. Figure 4.13 shows the myelin water fraction calculated from the two scans assuming an0m = 180°. As expected, the regions of the image with a < 180°, for example, the lower part of the image, show a greatly reduced myelin water fraction. 125 CHAPTER 4. NON-180 REFOCUSING PULSES 4.4. MULTI-SLICE, MULTI-ECHO Figure 4.13: The myelin water maps calculated from the pulse sequence with a refocusing pulse flip angle prescribed as 180° (left) and 150° (right). The NNLS algorithm was applied assuming, incorrectly, the refocusing pulse flip angle was 180°. 4.3.5.3 Conclusions The myelin water map estimated from both scans was similar to myelin water maps estimated from an optimized MESE0 sequence. Quantitatively the estimated percent myelin water, refocusing pulse flip angle, and geometric mean T2 from the two scans corresponded very well. The voxels in the brain with CSF tended to be poorly fit (high x2 misfit) which could be due to flow effects or an under-sampled T2 decay curve. 4.4 Multi-Slice, Multi-Echo 4.4.1 Introduction The current multi-echo acquisition method, MESE0, is able to collect only a single slice within the given scan time, due to the hard, composite pulses. The MESEC pulse sequence can be used to acquire multiple slices within one slab with the same SNR as the MESE„ sequence. This is important when viewing pathology, for example multiple sclerosis lesions, as they extend across several slices. Acquiring more slices will enable better interpretation of the MS lesions. The ability to collect multiple slices in a clinically feasible scan time will be an important step in the 126 CHAPTER 4. NON-180 REFOCUSING PULSES 4.4. MULTI-SLICE, MULTI-ECHO use of the myelin water imaging. I scanned a volunteer using the 3D MESEC sequence. Six slices were acquired, within the brain, as well as a slice acquired with the MESE0 for comparison. Myelin water maps were calculated for each slice of the MESEC sequence and the single slice of the MESE0. The myelin water maps were compared as well as the myelin water fraction from select regions. 4.4.2 Methods Two sets of data were collected: first, a standard 32 echo MESE0 with 4 averages and TE = 10 ms; and second, a 6 slice slab, 42 echo 3D-MESEC sequence with 1 average and TE = 10.1 ms. Other parameters for both sequences were: FOV = 24 cm, TR = 3000 ms, and 5 mm thick slices. Phantoms were placed on either side of the head (top and bottom phantoms on right side of image were the short T2 phantoms). Both sets of multi-echo data were filtered with the anisotropic diffusion filter algorithm (3 iterations) and the myelin water maps were calculated, per voxel, with the small norm NNLS model described in Section 1.7.1.2. A regional T2 analysis compared decay curve parameters estimated by a mono-exponential and NNLS from the MESE0 dataset and the corresponding slice of the 6 slice slab MESEC data set. Finally, the myelin water map was calculated such that the refocusing pulse flip angle was calculated per voxel. The resulting myelin water map was compared to that calculated assuming a = 180° throughout the image. 4.4.3 Results and Discussion 4.4.3.1 Images Figure 4.14 shows the TE = 10 ms image (top) and the TE = 10.1 ms images (second and third rows) of the brain. The outer slices of the six slice slab (first image of row 2 and last image of row 3) have aliasing artifact due to the slab selective excitation pulse. The aliasing artifact occurs in the slice phase encode direction as the head extends beyond the FOV and spins are not mapped to the correct location. The images are normally not reconstructed during a clinical scan, but were in this work for completeness. The myelin water maps estimated from the MESE0 sequence and each slice of the MESEC 127 CHAPTER 4. NON-180 REFOCUSING PULSES 4.4. MULTI-SLICE, MULTI-ECHO sequence are shown in Figure 4.15. As can be seen in Figure 4.16 there is a small increase in the myelin water fraction in the top right portion of the image. 128 CHAPTER 4. NON-180 REFOCUSING PULSES 4.4. MULTI-SLICE, MULTI-ECHO Figure 4.14: The TE = 10 ms image from the four average, 32 echo MESE0 sequence (top row). The six slices of the six slice slab MESEC acquisition (bottom two rows). The outside two slices of the slab selective sequence (second row, left image and third row, right image) have an aliasing artifact which results in a ghost-like overlay of another slice. 129 CHAPTER 4. NON-180 REFOCUSING PULSES 4.4. MULTI-SLICE, MULTI-ECHO Figure 4.15: Myelin water maps calculated using the small norm NNLS solution based on data from the 6 slice, 42 echo pulse sequence MESEC. The two outside slices (top left and bottom right) were not expected to be reasonable as they had significant (as expected) aliasing artifacts. The images were thresholded based on other parameters calculated (global variance, x2 and global density) which is why the two phantoms on the right side were thresholded. 130 CHAPTER 4. NON-180 REFOCUSING PULSES 4.4. MULTI-SLICE, MULTI-ECHO CHAPTER 4. NON-180 REFOCUSING PULSES 4.4. MULTI-SLICE, MULTI-ECHO 4.4.3.2 Regional Analysis Regions were selected in white matter and gray matter areas of the earliest echo image. The signal intensity of each region was averaged, for each echo, to calculate a T2 decay curve, the NNLSQ algorithm was applied to the decay curve. The short T2 and medium T2 parameters were calculated for various regions of interest (left and right): frontal white matter (frontaLwm), genu of the corpus callosum (genu_wm), posterior white matter (post_wm), mid-brain white matter (mid_wm), peri-ventricular gray matter (perivent_gm), and posterior gray matter (post-gm). The myelin water fraction calculated from data acquired with the MESEC sequence, was always higher than the myelin water fraction calculated from the data acquired with the MESEG sequence (Table 4.12). The decay curves were further assessed by comparing curves acquired from the two pulse sequences. One example T2 decay curve is shown in Figure 4.17. In general, the T2 decay curve acquired with the slab selected, multi-slice MESEC sequence decayed faster than the T2 decay curve acquired with the MESE0 sequence. The MESEC sequence may have subtle differences in the acquisition method to accommodate the acquisition of multiple slices within the slab. 132 CHAPTER 4. NON-180 REFOCUSING PULSES 4.4. MULTI-SLICE, MULTI-ECHO NNLS Structure Sequence Short p Short T2 Med p Med T2 frontal_wm_l MESE0 (180) 47.258 ( 4.6%) 15.000 931.387 (91.5%) 79.337 frontaLwmJ MESEC (180) 78.021 ( 8.0%) 36.109 863.296 (88.1%) 80.566 frontal_wm_l MESEC (174) 83.473 ( 8.5%) 34.430 873.974 (88.8%) 81.122 frontal_wm_r MESE0 (180) 36.290 ( 3.5%) 15.000 977.200 (93.7%) 79.755 frontaLwm_r MESEC (180) 138.347 (13.5%) 24.656 885.655 (86.5%) 83.284 frontal_wm_r MESEC (180) 138.347 (13.5%) 24.656 885.655 (86.5%) 83.284 genu_wmJ MESE0 (180) 94.612 ( 9.2%) 22.819 908.342 (88.0%) 77.340 genu-wmJ MESEC (180) 150.983 (15.2%) 37.115 837.326 (84.4%) 78.994 genu_wmJ MESEC (180) 150.983 (15.2%) 37.115 837.326 (84.4%) 78.994 genu_wm_r MESE„ (180) 99.003 ( 9.3%) 15.127 947.429 (89.4%) - 78.345 genu_wm_r MESEC (180) 179.500 (17.8%) 35.845 798.541 (79.3%) 81.258 genu_wm_r MESEC (173) 179.886 (17.8%) 35.065 799.187 (78.9%) 80.791 post_w'mJ MESE0 (180) 88.007 ( 8.5%) 17.971 942.529 (91.5%)- 92.041 post_wm_l MESEC (180) 244.997 (24.4%) 35.921 686.287 (68.4%). 96.742 post_wm_l MESEC (180) 244.997 (24.4%) 35.921 686.287 (68.4%) 96.742 post_wm_r MESE0 (180) 70.068 ( 6.7%) 15.291 808.996 (77.5%) 87.281 , post_wm_r MESEC (180) 224.328 (22.3%) 40.148 589.941 (58.5%) 86.406 post_wm_r MESEC (180) 224.328 (22.3%) 40.148 589.941 (58.5%) 86.406 mid_wm_l MESE0 (180) 177.964 (17.2%) 29.219 607.018 (58.8%) 82.458 mid_wmJ MESEC (180) 215.597 (22.4%) 39.133 500.218 (52.0%) 79.342 mid_wmJ MESEC (180) 215.597 (22.4%) 39.133 500.218 (52.0%) 79.342 mid_wm_r MESE0 (180) 168.491 (16.7%) 29.417 606.233 (60.0%) 85.955 mid_wm_r MESEC (180) 228.307 (23.1%) 30.727 567.696 (57.3%) 86.775 mid_wm_r MESEC (178) 229.005 (23.1%) 30.709 567.653 (57.3%) 86.831 perivent_gm_l MESE0 (180) 0.000 ( 0.0%) — 1161.671 (100.0%) 76.568 perivent_gm_l MESEC (180) 0.000 ( 0.0%) — 1196.713 (100.0%) 68.768 perivent-gmJ MESEC (172) 0.000 ( 0.0%) — 1200.621 (100.0%) 68.533 perivent_gm_r MESE0 (180) 17.381 ( 1.5%) 15.000 1173.212 (98.5%) 81.559 perivent_gm_r MESEC (180) 0.000 ( 0.0%) — 1163.000 (100.0%) 73.643 perivent_gm_r MESEC (180) 0.000 ( 0.0%) — 1163.000 (100.0%) 73.643 post_gm_l MESE0 (180) 0.000 ( 0.0%) — 1018.847 (100.0%) 86.499 post_gm_l MESEC (180) 197.091 (19.1%) 42.066 745.337 (72.2%) 79.251 post_gm_l MESEC (159) 139.822 (13.1%) 22.560 926.926 (86.9%) 79.282 post_gm_r MESE0 (180) 0.000 ( 0.0%) — 873.843 (83.6%) 84.274 post_gm_r MESEC (180) 199.367 (19.3%) 42.229 607.816 (59.0%) 79.542 post_gm_r MESEC (172) 209.809 (20.3%) 41.925 599.110 (58.0%) 80.114 Table 4.12: Decay curve parameters estimated for white matter regions throughout the brain. The small norm NNLS solution was calculated from the MESE0 and MESEC data. For each triple of lines, the parameters were estimated based on the assumption that the refocusing pulse flip angle was a = 180°, the parameters on the bottom line were estimated based a best estimate of the flip angle (shown in brackets). The NNLS short and medium parameters were estimated from bins 10 < T2 < 50 ms and 50 < T2 < 120 ms, respectively. 133 CHAPTER 4. NON-180 REFOCUSING PULSES 4.4. MULTI-SLICE, MULTI-ECHO 100 200 300 TE (ms) 400 500 co D •g CO CD CC 10 0 -10 ( 1000 800 £ 600 CD CO |> 400 co 1 1 / \ : i i i i i ! 1 1 V ,4t/«>^^.w-.,v..J,-1 1 1 50 100 150 200 250 300 350 400 450 200 m u a 1 1000 800 O Measured _ . MESEc Fitted (a=180) .... MESE Fitted (a=176.0) 600 400 20 40 8»»p«88tetaaa 60 "0 50 100 150 200 250 300 350 400 450 TE (ms) Figure 4.17: Region (genu_wm_l) shown on the top figure, T2 decay curves acquired by the MESE0 and MESEC (middle figure). NNLS fits to the MESEC data, one assuming the re focusing pulse flip angle was 180° and the other where the optimal flip angle was calculated (bottom figure). 134 CHAPTER 4. NON-180 REFOCUSING PULSES 4.4. MULTI-SLICE, MULTI-ECHO 4.4.3.3 Flip Angle Estimation The slice from the slab MESEC that corresponded to the slice from the single-slice MESE0 was further analyzed. The myelin water map was calculated per voxel in two ways: 1) assuming that the refocusing pulse flip angle was a = 180° everywhere (mwm180) and 2) estimating the refocusing pulse flip angle per voxel (mwme). Figure 4.18: Maps of the refocusing pulse flip angle calculated from the 6 slice, 42 echo pulse se quence MESEC. The two outside slices (top left and bottom right) were not expected to be reasonable as they had significant (as expected) aliasing artifacts. 4.4.4 Conclusions Moving from single-slice to multiple slices (slab selection) per acquisition will greatly improve the utility of the myelin water mapping. To acquire multiple slices per acquisition slice-selective refocusing pulses are necessary. The flip angle of slice-selective refocusing pulses is inaccurate and must be accounted for during T2 estimation from a multi-echo decay data. A set of six slices, each with 42 echoes, was acquired. Myelin water maps and flip angle maps were calculated from each of the six slices. Myelin water maps of the inner four slices were similar to myelin water 135 CHAPTER 4. NON-180 REFOCUSING PULSES 4.5. NON 90° EXCITATION PULSE maps calculated from MESE0. The T2 estimated from the decay curves of regions within the six slices had a faster decay than expected, but is believed to be a problem with the pulse sequence and not the theory. 4.5 Non 90° Excitation Pulse The standard MRI pulse sequence contains two types of RF pulses, an excitation pulse, used to rotate the magnetization away from the longitudinal axis, and a set of refocusing pulses, used to refocus the magnetization in the transverse plane. Both sets of pulses are affected by imperfections in the hardware. The signal intensity along a decay curve is a function of the excitation pulse (ae), the refocusing pulse (ar) and the tissue parameters, p, T2, Ti. The goal of this work is to determine the effect of an "imperfect" excitation pulse on quantification of the tissue parameters and flip angles. 4.5.1 Decay Curve Calculation Given the tissue parameters T2, Ti, p and the pulse sequence parameters TE, ae and ar, the decay curves were calculated by using a modified version of Hennig's work [78]. The magnetization vector M = [F F* Z Z*]T where F and F* are transverse magnetization components, Z and Z* are longitudinal components. The components without a superscript * are dephasing magnetization and components with a superscript * are rephasing magnetization. The initial magnetization is defined to be M = [0 0 1 0]T and was then rotated by an angle a based on the rotation matrix: R(cx) = cos2 (a/2) sin2 (a/2) sin (a) 0 sin2 (a/2) cos2 (a/2) 0 sin (a) -1/2 sin (a) 1/2 sin (a) cos (a) 0 1/2 sin (a) -1/2 sin (a) 0 cos (a) (4.19) 4.5.2 Simulation Methods Three sets of simulations were done to determine the effect of knowing or not knowing the exci tation pulse flip angle ae. The first simulation (Section 4.5.3) had an excitation pulse that varied and the decay curve parameters were calculated assuming the excitation pulse was ae = 90° and 136 CHAPTER 4. NON-180 REFOCUSING PULSES 4.5. NON 90° EXCITATION PULSE the refocusing pulse was calculated. The second simulation (Section 4.5.4) set ae = ar/2 but assumed the excitation pulse ae = 90°. The third simulation (Section 4.5.5) had the excitation and refocusing pulses varied (but still ae = ar/2) and calculate both ae and aT. 4.5.3 Simulation 1: Inaccurate Excitation and Accurate Refocusing, Assume both Accurate This simulation will assess the effect of having an inaccurate excitation pulse but a perfect refo cusing pulse ar = 180°. Methods A set of true decay curves were calculated with the excitation pulse flip angle set' to ae — 90°,85°,...,50°. Other parameters were: ar = 180, 32 echoes, Ti = 1000 ms, T2 = 80 ms, TR = oo, echo spacing = 10 ms and p = 1000. For each setting of the excitation pulse, 500 realizations of quadrature noise (Gaussian on two channels) with SNR = 200 were added to the true decay curve. For each decay curve with noise, ar was calculated, as well as, p, T2 and x2 (the misfit). Results and Discussion This simulation assessed the parameter estimation when only the excitation pulse was varied and the curve fitting was done assuming only the refocusing pulse was variable. The primary result is that all parameters were estimated as accurately when ae ^ 90° as when ae = 90°, except for p. The average p decreased as a function of sin(ae). The x2 misfit remained constant at approximately 27, therefore the decrease in p corresponds very well to a decrease in ae. This would imply that if the only ae or p changed, they could not be distinguished. The acquisition of data with an imperfect excitation pulse and perfect refocusing pulse would be unlikely. This simulation was done for completeness and for comparison with Simulation 2 and Simulation 3. ' 137 CHAPTER 4. NON-180 REFOCUSING PULSES 4.5. NON 90° EXCITATION PULSE True ae (°) ar (°) P T2 ( ms) X2 90 176 (4.3) 1011.2 (17.2) 79.2 (2.3) 27.9 (7.5) 85 176 (4.5) 1008.9 (18.8) 79.2 (2.0) 27.6 (7.5) 80 176 (4.5) 996.4 (16.8) 79.1 (2.0) 27.9 (7.5) 75 176 (4.5) 978.3 (18.2) 79.0 (1.7) 27.4 (7.0) 70 176 (4.5) 951.4 (16.7) 79.0 (1.9) 27.7 (7.7) 65 176 (5.0) 917.4 (17.5) 79.1 (2.6) 27.8 (7.7) 60 176 (4.8) 877.6 (18.3) 79.0 (2.7) 26.8 (6.9) 55 176 (4.6) 830.4 (17.7) 79.0 (2.6) 27.3 (7.2) 50 175 (5.2) 777.8 (17.4) 78.9 (3.5) 26.8 (7.0) Table 4.13: The excitation pulse was varied from 90° to 50° and the parameters were calculated assuming only the refocusing pulse was incorrect. The primary effect of an imperfect excitation pulse is an overall decrease in the decay curve. 4.5.4 Simulation 2: Inaccurate Excitation and Inaccurate Refocusing, Assume Accurate Excitation This simulation will assess the effect of having ae = aT/2 but assuming ae = 90° for the curve fitting. This simulation closer to reality when acquiring multi-exponential data, as spatial nonuni-formities in the transmit power will affect both the excitation and refocusing pulses. Methods The parameters for the simulation are: 32 echoes, Tj = 1000 ms, T2 = 80 ms, TR = oo, echo spacing =10 ms and p — 1000. The excitation pulse was varied: ae = 90°, 85°,..., 50° (and ar = 2ae in every case). For each setting of ae-ar, 500 realizations of quadrature noise (Gaussian on two channels) with SNR = 200 were added to the true decay curve. For each decay curve with noise, ar was calculated as well as p, T2 and x2 (the misfit). The excitation pulse flip angle was assumed to be ae = 90°. Results and Discussion The estimated parameters are shown in Table 4.14. The refocusing pulse flip angle aT was estimated reasonably accurately, but all other parameters had problems. As in the first simulation, p decreased almost as a function of sin (ae). As well, T2 decreased as ae decreased. The x2 misfit increased as ae decreased, therefore the curve fits got progressively worse as the excitation and 138 CHAPTER 4. NON-180 REFOCUSING PULSES 4.5. NON 90° EXCITATION PULSE refocusing pulses decreased. True ae/ar (°) ar (°) P T2 (ms) x2 90/180 176 (4.3) 1011.2 (17.2) 79.2 (2.3) 27.9 (7.5) 85/170 173 (6.0) 1004.4 (17.3) 79.6 (2.4) 27.9 (7.6) 80/160 167 (9.6) 983.5 (16.5) 80.1 (2.4) 38.2 (9.2) 75/150 151 (6.3) 974.5 (20.7) 79.4 (2.0) 57.9 (11.7) 70/140 138 (2.3) 955.0 (15.2) 78.9 (2.7) 77.0 (16.0) 65/130 127 (1.9) 925.7 (9.2) 78.6 (2.0) 92.1 (17.4) 60/120 117 (1.7) 902.2 (8.5) 77.2 (1.3) 103.1 (18.7) 55/110 107 (1.4) 871.6 (9.0) 76.3 (1.7) 113.6 (18.7) 50/100 97 (1.3) 852.8 (10.2) 74.0 (0.9) 123.1 (20.4) Table 4.14: The excitation and refocusing pulses were varied with ae = 90°, 85°,..., 50° and ar = 2ae. The curve fitting calculated ar, T2, p and x2 assuming ae = 90°. The true T2 = 80 ms and p = 1000. 4.5.5 Simulation 3: Inaccurate Excitation and Inaccurate Refocusing The realistic scenario is the excitation pulse and refocusing pulse are both inaccurate and need to be calculated during T2 estimation. This simulation will assess the effect of having ae = ar/2 and calculating both ae and ar along with the other parameters. Methods The parameters for the simulation were: 32 echoes, Ti = 1000 ms, T2 = 80 ms, TR = oo, echo spacing = 10 ms and p = 1000. The excitation pulse was varied: ae = 90°, 85°,..., 50° (and ar = 2ae in every case). For each setting of ae, 500 realizations of quadrature noise (Gaussian on two channels) with SNR = 200 were added to the. true decay curve. For each decay curve with noise, the ar, ae, p, T2 and %2 (the misfit) were calculated. The decay curve was fit assuming the ratio ae/ar = 1/2. Results and Discussion All parameters (shown in Table 4.15) were estimated accurately, including ae and ar. The stan dard deviation of the estimated T2 and p increased by a factor of two as ae = 90° to ae = 50°. The standard deviation of the estimated parameters doubled as the excitation/refocusing pulse flip angle halved. 139 CHAPTER 4. NON-180 REFOCUSING PULSES 4.5. NON 90° EXCITATION PULSE True ae/ar (°) a S(°) ar (°) P T2 (ms) x1 90/180 88 (2.4) 176 (4.8) 1013.1 (17.8) 79.2 (1.8) 28.0 (7.6) 85/170 86 (2.3) 172 (4.5) 1007.2 (18.0) 79.3 (2.1) 27.7 (7.6) 80/160 80 (1.7) 161 (3.3) 1004.9 (17.4) 79.5 (1.8) 27.8 (7.5) 75/150 76 (0.8) 152 (1.6) 998.1 (14.7) 79.7 (1.7) 28.1 (7.6) 70/140 70 (1.0) 140 (2.0) 1013.9 (24.8) 78.8 (1.7) 27.9 (7.4) 65/130 65 (0.6) 130 (1.1) 1008.3 (18.2) 79.3 (2.7) 27.6 (7.3) 60/120 60 (0.7) 120 (1.5) 1011.5 (23.5) 79.1 (2.4) 27.5 (7.6) 55/110 55 (0.5) 110 (1.0) 1005.4 (20.2) 79.5 (3.0) 28.0 (7.5) 50/100 50 (0.6) 100 (1.3) 1011.3 (28.3) 79.3 (3.6) 26.9 (6.8) Table 4.15: The excitation and refocusing pulses were varied with ae = 90°, 85°,..., 50° and ar = 2ae. The curve fitting calculated ae, ar, T2, p and x2 assuming ae/ar = 1/2. 4.5.6 Discussion The accuracy of the excitation and refocusing pulse each affect the estimation of the decay curve parameters. The accuracy of the refocusing pulse flip angle is more important and more difficult to attain than the accuracy of the excitation pulse flip angle. The primary reason the refocusing pulse flip angle is more important is there is an accumulation of errors in the magnetization along a train of refocusing pulses. It is exactly this accumulation of errors which I have exploited to estimate the refocusing pulse flip angle. The excitation pulse is applied once at the beginning of the train and therefore there is not an accumulation effect. The estimated T2, when the excitation pulse flip angle was assumed to be ae = 90°, was as accurate and precise (Table 4.14) as the T2 estimated when the excitation pulse flip angle was estimated along with the refocusing pulse flip angle (Table 4.15). A second reason the excitation pulse is less important than the refocusing pulse is related to the linearity of the Bloch Equations. The magnetization profile deviates from the Fourier transform of the pulse profile when the flip angle exceeds 30° [82]. The linearity of the Bloch Equations is reasonable up to 90°, but larger errors in the shape of the slice profile occur as a increases to 180°. Therefore, excitation pulses near 90° produce a slice profile closer to the desired than refocusing pulses near 180°. 140 CHAPTER 4. NON-180 REFOCUSING PULSES 4.6. NON-CONSTANT TRAIN 4.5.7 Conclusions The excitation pulse flip angle, as well as the refocusing pulse flip angle, is affected by the in-homogeneity and imperfections in the hardware. The inhomogeneities and imperfections lead to spatially varying flip angles. The simulations in this section were used to assess the accuracy and precision of the decay curve parameter estimation in the presence of inaccurate pulses. As was seen, the accuracy and precision of the decay curve parameter estimations, if the flip angles were fit, were similar to the accuracy and precision from a sequence with ae = 90° and ar = 180°. 4.6 Non-Constant Train of Flip Angles The work previously showed that T2 decay curve parameters are accurately and precisely deter mined when the train of refocusing pulses is less than 180°. As well, as the refocusing pulse flip angle decreased, the signal to noise in the data decreased. This decrease in SNR will lead to increased variability of the measured decay curve parameters. The short T2 is of primary interest and is determined from the signal in the decay curve with TE < 50 ms. Therefore, the refocusing pulses with TE < 50 ms should be near 180° to retain the signal to noise, but the refocusing pulses where TE > 50 ms can be decreased. One method to do this is to decrease the amplitude of latter pulses and have the earlier pulses remain at 180°. This would provide a reasonable trade-off between decreased power deposition and high SNR in the early echoes. The decrease in amplitude was accomplished by defining a scaling vector, sa, in which each element can vary between 0 and 1. Each element in the scaling vector sQ is then multiplied by 180° to create the vector of refocusing pulse flip angles. There are an infinite number of combinations to scale the flip angle of the refocusing pulses within a pulse train. For the simulations and scanner results of a 16 echo pulse train, I analyzed four possibilities: 16 • sLi}='i.o,. r., i.o 16 • si2) = 0.5,. , 0.5 13 • s^3) = 1.0,1.0,1.0,0.5,0.5* .. ,0.5 141 CHAPTER 4. NON-180 REFOCUSING PULSES 4.6. NON-CONSTANT TRAIN • s{a] = 1.0000, 0.9048, 0.8187, 0.7408, 0.6703, 0.6065, 0.5488, 0.4966, 0.4493, 0.4066, 0.3679, 0.3329, 0.3012, 0.2725, 0.2466, 0.2231. The first set of scalings, Sa \ is a standard pulse train of 180° refocusing pulses. The second set (2) of scalings, sa , is a pulse train of 90° refocusing pulses. Both of these were used for reference when testing and s$. The third set of scalings (Figure 4.19), is an amalgamation of the first two, it retains 180° refocusing pulses for the first part of the train to retain the higher SNR and the last part of the train is 90° to decrease the power deposition. The fourth set of scalings (Figure 4.20), SQ\ was defined as Sa\i) = exp(—i/10) for i = 0,1,..., 15. The denominator within the exponential was chosen so the decay in the scalings was not too fast. This last one resulted in a smooth decrease in the amplitude of the refocusing pulses throughout the train. ,(4) 1 1 1 1 1 :A""." A—A .A—-jx—/L (1 -4 U J m m L 4 1 1 „ i fl r I fl fl i M 1 M Ml II fl fl • ... u 1 RF G G DACQ 10 20 30 Time (ms) 40 50 60 70 Figure 4.19: A multi-echo, spin echo pulse sequence in which the first three refocusing pulses have a flip angle of 180° and the remaining thirteen have a flip angle of 90° (only six of the 16 shown). A set of simulations were used to assess the accuracy and precision of decay curve parameter estimation from decay curves collected with non-constant refocusing pulse trains. A pulse sequence was modified to allow arbitrary scaling of the refocusing pulses through the train and was used to assess the effect on the phantoms (Section 1.8). 142 CHAPTER 4. NON-180 REFOCUSING PULSES 4.6. NON-CONSTANT TRAIN RF G • \~ 1 ! 1 1 1 ;A /A—„..A.. . A—-A M fl L -4 L -4 f--4 i--4 l_ -4 I \ I „ 1 , I „ l i r I •fl fl II M I M II Ml u ../ v • • v -DACQ 10 20 30 Time (ms) 40 50 60 70 Figure 4.20: A multi-echo, spin echo pulse sequence in which the refocusing pulse amplitudes are decaying exponentially. 4.6.1 Simulations Determine the accuracy and precision of calculating the proton density, T2 and flip angle from a decay curve that was "collected" with non-constant refocusing pulses through the train. 4.6.1.1 Methods I created a set of simulations with p = 1000, T2 = 80 ms, ae = 90°. The four sets of refocus ing pulse trains, defined above, were used. For each decay curve with noise, the decay curve parameters were calculated and compared to truth. For the decay curves with non-constant refocusing pulses, the relative scale between each pulse was assumed to be known. For example, the relative scale between pulses in the third refocusing pulse train was assumed to be 1.0,1.0,1.0,0.5,... , 0.5. 4.6.1.2 Results The parameters a, p, and T2 were estimated from the simulated data for each of the scalings and are shown in Table 4.16. There was negligible difference between refocusing pulse train s|,3' and 143 CHAPTER 4. NON-180 REFOCUSING PULSES 4.6. NON-CONSTANT TRAIN SQ ^ compared to the constant refocusing pulse train of ar = 180°. As expected, p and T2 had a (2) higher standard deviation in simulation sa . Simulation ar p T2 (ms) x2 176 (4.5) 1010.6 (19.6) 78.7 (2.6) 13.0 (5.4) sk2) 90 (1.4) 1005.2 (27.3) 78.4 (5.5) 12.8 (4.9) s(a] 176 (4.5) 1010.6 (19.8) 78.8 (2.8) 12.6 (4.8) sL4) 176 (4.4) 1010.0 (20.1) 78.8 (2.9) 12.8 (5.1) Table 4.16: Decay curve parameters calculated for each of the four refocusing pulse trains. For simulation s^ and s^, the refocusing pulse flip angle, ar in the second column is to be multiplied by the assumed relative scale to obtain the train of refocusing pulse flip angles. The true T2 was 80 ms and p = 1000. 4.6.2 Scanning The four sets of scalings were implemented in a single slice, 16 echo MESEC pulse sequence. The nickel/agarose phantoms were scanned with the pulse sequence and the estimated T2 were compared to the T2 estimated from the MESE„ sequence. 4.6.2.1 Methods Scans of the nickel/agarose phantoms were acquired with each of the set of scalings of the refocus ing pulse flip angles. Other scan parameters were: FOV = 16 cm, TE = 12.12 ms, TR = 3000 ms, and 16 echoes. The scalings were defined in a file on the scanner computer. At the beginning of the scanO routine, the amplitude of the ith RF pulse (ia_rf 2) was scaled by s\ defined in the file. 4.6.2.2 Results and Discussion The decay curves for each phantom were fit using NNLS and the results for two phantoms are shown in Figures 4.21 and 4.22. The fitted decay curve was very close to the measured decay curve for both s^3) and There was a small discrepancy between the measured and fitted decay curves, in particular with s^, at the fourth echo, which is the first echo after the large change in refocusing pulse flip angle within the train. 144 CHAPTER 4. NON-180 REFOCUSING PULSES 4.6. NON-CONSTANT TRAIN a: 180 16 a: 90 16 800 ^600 CO c CD | 400 co c a) ^ 200 0 600 "0 50 100 150 TE (ms) a: 180, 162.9, 40 200 50 100 150 TE (ms) a: 180, 180, 180, 90 200 13 800 ^600 CO c 5 400 CO c 200 0 800 ^600 c/> c 5 400 cvj c g> ^ 200 "0 50 100 150 TE (ms) 200 0 \\ c 0 50 100 150 TE (ms) 200 Figure 4.21: Decay curve fit of the bottom, middle phantom, for refocusing pulse trains of 180°y (top, left), 90°N (top, right), exponential decay of refocusing pulse flip angles (bottom, left), and 180°, 180°, 180°, 90°,..., 90° (bottom, right). Measured data shown with whiskers and fitted data shown as circles. 145 CHAPTER 4. NON-180 REFOCUSING PULSES 4.6. NON-CONSTANT TRAIN a: 180 16 a: 90 16 600 500 -t—t cn c 400 CD I * _c 300 "co o)200 CO 100 0 600 500 >^ "w c 400 CD | • 300 "CO D)200 cn 100 0 "0 50 100 150 TE (ms) a: 180, 162.9, 40 200 350 300 | 250 I 200 _c To 150 c 100 50 0 0 50 100 150 TE (ms) a: 180, 180, 180, 90 •(t-C^-JHll 600 500 -4—* 'co c 400 0 c 300 IB S)200 CO 100 "0 50 100 150 TE (ms) 200 0 200 13 50 100 150 TE (ms) 200 Figure 4.22: Decay curve fit of the bottom, right phantom, for refocusing pulse trains of 180°y (top, left), 90^ (top, right), exponential decay of refocusing pulse flip angles (bottom, left), and 180°, 180°, 180°, 90°,..., 90° (bottom, right). Measured data shown with whiskers and fitted data shown as circles. 146 CHAPTER 4. NON-180 REFOCUSING PULSES 4.6. NON-CONSTANT TRAIN The T2 of each phantom was calculated from each pulse sequence and is shown in Table 4.17. The exponentially decaying refocusing pulse train typically underestimated the T2 by approxi mately 10%. Phantom Sequence T2 (ms) x2 bottomJeft 24.6 162.0 2.0 bottomJeft 20.9 81.0 1.0 bottomJeft s(3) 24.6 154.0 7.0 bottomJeft S(4) 21.2 142.0 2.8 bottom_middle s(D 83.2 157.0 1.0 bottom_middle S(2) 70.0 85.0 0.6 bottom_middle S(3) 82.8 156.0 38.9 bottom_middle S(4) 75.0 145.0 7.7 bottom_right S(D 23.0 166.0 1.5 bottom_right 17.3 80.0 1.1 bottom_right S(3) 22.5 156.0 3.4 bottom_right S(4) 20.3 148.0 3.6 topJeft SW 299.2 156.0 0.4 topJeft S(2) 109.9 83.0 1.4 topJeft S(3) 374.9 155.0 58.1 topJeft S(4) 280.4 146.0 16.0 Table 4.17: The T2 and refocusing pulse flip angle (for the first echo) calculated from each train of refocusing pulses. The T2, from the MESE0 sequence were: 25 ms (bottom, left), 85 ms (bottom, middle), 24 ms (bottom, right), and 280 ms (top, left). One thing that was a little odd was the curves for s^4) did not follow as well as I would have expected. I noticed from previous work that the calculated refocusing pulse flip angle given ao = 180° is approximately acaic = 157° and given ao — 90° is approximately ftca]c — 83°. The ratio 157°/180° = 0.8722 does not coincide with 83°/90° = 0.9222. This would imply that even though I scale the amplitude of the RF pulses by 1.0 and 0.5 when I have two different scalings in the same pulse train, I don't know if I can expect to calculate the same relative scalings. Another way to put it, I don't believe the amplitude scaling of the RF pulses results is linearly correlated with the actual applied RF pulse flip angle. 147 CHAPTER 4. NON-180 REFOCUSING PULSES- 4.6. NON-CONSTANT TRAIN S90 x2 T2 (ms) 0.50 3.434 22.47 0.53 2.551 22.39 0.55 1.851 22.38 0.57 1.387 22.48 0.60 0.987 22.50 0.62 0.735 22.54 0.65 0.603 22.74 0.68 0.596 22.81 0.70 0.647 23.00 0.72 0.780 23.21 0.75 0.969 23.53 Table 4.18: The resulting x2 misfit given the scaling of the RF pulses was 1.0,1.0,1.0,590 from s^4). The small misfit was when the scaling of the last 13 echoes was set to 0.68 rather than 0.5, as expected. Row with the minimum x2 is in bold. Figure 4.23: The measured data (blue circles), fitted data assuming the scaling of the last 13 echoes was 0.5 (red dashed line) and fitted data assuming scaling was 0.68 (black dot-dashed line). Residuals of fit to measured data shown in top panel. 148 CHAPTER 4. NON-180 REFOCUSING PULSES 4.7. CONCLUSIONS 4.6.3 Conclusion The amplitude of the RF pulses in the train can be scaled to have any amplitude and the T2 tends to be reasonably accurate. There does seem to be a non-linear correlation between the input scale parameter and the calculated scaling which may cause some problems. Possibly, this could be included as another step in the fitting process or could be calibrated. 4.7 Conclusions Methods for T2 correction due to Bi inhomogeneities have relied on the acquisition of a Bi map or based on uniform phantoms. Other methods required two acquisitions to calculate the corrected mono-exponential T2, one acquisition to calculate the multi-echo decay curve and the second acquisition to estimate the Bi field. The method I propose requires only one acquisition, the multi-echo decay curve is sufficient to calculate T2 and a. Therefore, the T2 and Bi are inherently registered and the Bi accounts for coil loading. As well, the method can be incorporated into any curve fitting technique and, as shown above, is able to accurately quantify multi-exponential decay curves. High field systems (> 3T) have inhomogeneous Bi fields, due to the dielectric effect [83], and require high power deposition for 180° refocusing pulses [84]. Gareau et al. [85] compared T2 distributions of guinea pig brain at 1.5T and 4T and found both the short and medium T2 peaks had a faster decay at the higher field strength. The decrease in T2 at higher field would require shorter echo times. I have shown that decay curve parameters can be calculated in the presence of refocusing pulses with a < 180°. Another equally important observation is that quantitative T2 measurements can be estimated from a train of a refocusing pulses having a <C 180°. There has been a restriction on the length of the train of 180° refocusing pulses due to the power deposition in vivo. The methodology described enables quantitative T2 measurements to be done on high field systems and will account for the echo artifacts. 149 Chapter 5 Conclusions 5.1 Filtering Image processing, noise reduction filters are a simple and fast method to decrease the variability in local areas of an image. The goal of the work was to decrease the variability in local areas of the myelin water map by filtering the multi-echo data acquired by an MRI machine. Variability in the myelin water fraction over a region was a logical and simple measure to analyze the filters. Simulated multi-echo data was an important step in assessing the filters as there were no artifacts due to imaging, for example flow, head motion, and Fourier transform reconstruction. The in vivo myelin water maps looked much better if the multi-echo data was processed using a noise reduction filter. For myelin water maps calculated from the non-linear curve fitting algorithms, it is important to apply the noise reduction filter to the original multi-echo data. The image processing filters, presented in this work, are the more popular noise reduction filters used in MRI today. Image processing is an evolving field and there will likely be other filters which could prove to be better. The primary goals for a noise reduction filter, applied to myelin water quantification, should include: ability to smooth over homogeneous areas and retain the boundaries between the tissues, produce little or no bias in the mean myelin water fraction over a region, and reduce the variability of the myelin water fraction within a local region. 5.2 Linear Combination The linear combination method to calculate the myelin water fraction was shown to be robust and very fast, in addition to its extreme ease of implementation after the coefficients are determined. The short T2 filter was to select for T2 near 20 ms only, and suppressed all T2 above 70 ms. The myelin water maps were comparable to myelin water maps calculated by the NNLS algorithm. The linear combination method to analyze T2 from a multi-echo MRI dataset is restricted primarily to 150 CHAPTER 5. CONCLUSIONS 5.2. LINEAR COMBINATION the short T2 range. The restriction is due to the broad selection if longer T2 are of interest. The broad selection would negate the use of linear combination for T2 longer than 50 ms. Therefore, the curve fitting algorithms which provide information about the T2 distribution are still required if the T2 of interest is greater than 50 ms. The coefficients calculated using the non-linear algorithm, applied to the calculation of the myelin water fraction, were shown to provide robust myelin water fraction measures. Though, other sets of coefficients could provide more robust myelin water fraction measures. Other types of algorithms could be used to create another set of coefficients. One method is to fit the filter to a pre-defined filter and minimize the misfit between the two. Another method is to change the constraints, used in the algorithm presented, to create a more appropriate filter. One must remember that the underlying data is a set of exponentials, and exponentials are very non-orthogonal as a function of T2. Therefore, it may not be possible to create a better filter than what is presented. The short T2 filtered associated with the coefficients calculated in this chapter is slightly asymmetric in that the filter for 20 < T2 < 70 ms is less sharp than the filter for 10 < T2 < 20 ms. This provides a difficulty that a myelin water compartment which has an increase in the T2 is indistinguishable from a myelin water compartment that has only a decrease in the amplitude of the T2 peak. But, this has provided some interesting insight into one region of the brain. The internal capsules (at the level of the genu and splenium of the corpus callosum) typically appear brighter on a myelin water map, which implies a higher myelin water fraction. Myelin water maps, from the linear combination method, appear to have a decreased myelin water fraction in the same region. This would seem to imply there is a shift in the myelin water T2 toward higher T2 and is corroborated by results based on the NNLS algorithm. The shift toward higher T2 would be accounted for by a greater distance between the lipid bilayers compared to the distance between the lipid bilayers in other white matter. The linear combination for myelin water quantification provides a fast and reliable method to assess myelin water. The technique applied to 32 echoes proved to be robust relative to the NNLS technique, though, using fewer echoes may prove to be an important step. Shorter trains of RF pulses are more desirable than longer trains as there is less power deposition, but also, the "dead time" may be used in creative ways. A four-echo, interleaved in time and space may allow myelin water quantification from a volume of the brain. 151 CHAPTER 5. CONCLUSIONS 5.3. Bi/T2 5.3 B1/T2 The T2 distribution is calculated from a multi-echo decay curve. The decay curve is a combination ofthe primary echo as well as stimulated and indirect echoes. The stimulated and indirect echoes are a function of the refocusing pulse flip angle (a). Therefore, a modified T2 curve fitting algorithm should be able to estimate the T2 distribution, as well as the refocusing pulse flip angle, from a multi-echo decay curve. I presented an algorithm to estimate T2 and a from a mono-exponential decay curve and a modified NNLS algorithm, NNLSa, to estimate the T2 distribution and a from a multi-exponential decay curve. Noise simulations showed that the accuracy and precision were similar to the stan dard method to estimate T2. A decrease in the precision of the estimation was shown to be a function of the decreased a. Several extensions to the method were discussed: slab selected multi-slice, train of non-constant refocusing pulses and the effect of a non-90° excitation pulse. The method enables the use of slice-selective refocusing pulses, as the accuracy of the refocusing pulse is less impor tant as it is calculated during the curve fit. Therefore, multiple slices, within a slab, can be excited and acquired. This is an important extension to the current single-slice method. Multiple sclerosis lesions, for example, may span several slices, and therefore, to track lesions over time, it is important to obtain a view of normal and diseased tissue above and below a lesion. A sec ond strength of the method is the curve fit algorithm is able to estimate the T2 from a train of non-1800 refocusing pulses. A pulse sequence, with a train of non-1800 refocusing pulses deposits less power into a patient. Therefore, T2 estimation on a high field system (e.g., 3T) would be possible with a train of refocusing pulses less than 180°, though a trade-off may happen due to the decrease in T2 of the short and medium components. The amplitude of a refocusing pulse is not the only property that can be changed, another option is to change the phase of the pulse. Many groups have worked on phase schemes for more robust T2 decay curves, which have resulted in the Meiboom-Gill [86] extension to the original Carr-Purcell [87] pulses, and MLEV phase cycling [30]. The effectiveness of adjusting the phase of the RF pulse, rather than the amplitude, is that adjusting the phase can be done more precisely. The relationship between the amplitude of an RF pulse and the flip angle is based on a small tip angle assumption, which is weakly valid for a > 30° and less-so for a > 90° (this is likely 152 CHAPTER 5. CONCLUSIONS the cause of some discrepancies seen in Section 4.6, tailored RF pulses, such as those from the Shinnar-Le Roux algorithm, seek to work around these assumptions). Three weaknesses exist in the standard method to estimate the myelin water fraction: 1) low SNR, 2) single slice, and 3) long scan time. The first and third weaknesses are tightly coupled. The primary method to obtain a higher SNR is to do multiple averages. Decreasing the scan time, without a large decrease in the SNR, would increase the clinical feasibility of myelin water imaging. The standard assumption in quantitative T2 is the pulse sequence must have a long TR to allow magnetization relax back to thermal equilibrium along the longitudinal axis. Currently, the time between the final echo and the following excitation is at least one second. More efficient use of the dead time is required. One possibility is to decrease the dead-time between the final echo and the next excitation pulse. If possible, it would be best to find a method to force all magnetization back to the longitudinal axis. A second possibility, would be to apply the phase coherence information in Section 4.1.2 and to track the magnetization not only through a single pulse train, but through successive pulse trains. 5.4 Preferred T2 Quantification Protocol The quantification of myelin water is important for continued basic science research. The current protocol, using MRI, to quantify the myelin water fraction is to use an optimized 32-echo spin-echo pulse sequence to collect the multi-echo data, then quantify the short T2 using the NNLS algorithm. The multi-echo sequence and NNLS algorithm were shown, previously, to be a good combination to quantify the myelin water fraction. There are three motivations to change the current protocol: first, higher SNR data is required for more robust quantification of the myelin water fraction; second, averaging of single slice data is an inefficient use of time; third, more slices are necessary to assess Multiple Sclerosis lesions and their change over time. Therefore, given the motivations to change the current protocol, the following protocol is suggested. The preferred method to calculate the myelin water fraction is to collect a multi-slice, multi-echo MRI dataset with a long TR. The dataset would be processed with the 3D anisotropic diffusion filter which averages the signal based on edge information derived over the channels (channel filtering). The non-negative least squares algorithm, modified to calculate the Bi at each voxel, would be applied voxel-by-voxel to the noise-reduced, 3D dataset. 153 Bibliography [1] F. Bloch. Nuclear induction. Phys. Rev., 70(7,8):460-474, 1946. [2] N. Bloembergen, E. M. Purcell, and R. V. Pound. Relaxation effects in NMR absorption. Phys. Rev., 73:679-712, 1948. [3] P. Morell and W.T. Norton. Myelin. Scientific American, 242(5):88-117, 1980. [4] R.B. Dietrich and C. Hoffmann. Myelination and dysmyelination. In D.D. Stark and Jr. W.G. Bradley, editors, Magnetic Resonance Imaging, chapter 23, pages 1057-1080. Mosby Year Book, Inc., St. Louis, Missouri, second edition, 1992. [5] E.R. Kandel, J.H. Schwartz, and Jessel T.M., editors. The Principles of Neural Science. Prentice Hall, 3rd edition, 1985. [6] A.J. Barkovich. Concepts of myelin and myelination in neuroradiology. Am J Neuroradiol, 21:1099-1109, 2000. [7] P. Morell, editor. Myelin. 2nd edition, 1985. [8] T.J. Copeland. http://www.albany.net/~tjc/glossl-m.html#M. [9] G. R. Cherryman and F. W. Smith. Nuclear magnetic resonance in adrenoleukodystrophy: report of a case. Clin Radiol, 36(5):539-40, Sep 1985. [10] J. C. Masdeu, J. Moreira, S. Trasi, P. Visintainer, R. Cavaliere, and M. Grundman. The open ring, a new imaging sign in demyelinating disease. J Neuroimaging, 6(2): 104-7, Apr 1996. [11] M. H. Teicher, C. M. Anderson, A. Polcari, C. A. Glod, L. C. Maas, and P. F. Renshaw. Functional deficits in basal ganglia of children with attention-deficit/hyperactivity disorder 154 Bibliography shown with functional magnetic resonance imaging relaxometry. Nat Med, 6(4):470-3, April 2000. [12] V. Vasilescu, V. Simplaceanu E. Katona and, and D. Demco. Water compartments in the myelinated nerve. III. Pulsed NMR results. Experientia, 34:1443-1444, 1978. [13] RM Kroeker and Henkelman RM. Analysis of biological NMR relaxation data with continuous distributions of relaxation times. J. Magn. Reson., 69:218-235, 1986. [14] R. S. Menon and P. S. Allen. Application of continuous relaxation time distributions to the fitting of data from model systems and excised tissue. Magn Reson Med, 20:214-227, 1991. [15] A. L. MacKay, K. P. Whittall, K. S. Cover, D. K. B. Li, and D. W. Paty. In vivo visualization of myelin water in brain from mr spin-spin relaxation measurement. In Proceedings of the RSNA, page 142. RSNA, 1991. [16] A. L. MacKay, K. P. Whittall, K. S. Cover, D. K. B. Li, and D. W. Paty. In vivo T2 relaxation measurements of brain may provide myelin concentration. In Proceedings of the SMRM, page 917. SMRM, 1991. [17] W. A. Stewart, A. L. MacKay, K. P. Whittall, G. R. W. Moore, and D. W. Paty. Spin-spin relaxation in experimental allergic encephalomyelitis, analysis of CPMG data uing a non-linear least squares method and linear inverse theory. Magn. Reson. Med., 29:767-775, 1993. [18] A. MacKay, K. Whittall, J. Adler, D. Li, D. Paty, and D. Graeb. In vivo visualization of myelin water in brain by magnetic resonance. Magn Reson Med, 31(6):673-7, Jun 1994. [19] K. P. Whittall, A. L. MacKay, D. Graeb, R. Nugent, D. K. B. Li, and D. W. Paty. In vivo measurement of T2 distributions and water contents in normal human brain. Magn. Reson. Med., 37:34-43, 1997. [20] C Beaulieu, FR Fenrich, and PS Allen. Multicomponent water proton transverse relaxation and T2-discriminated water diffusion in myelinated and nonmyelinated nerve. Magn Reson Imaging, 16(10):1201-1210, 1998. 155 Bibliography [21] IM Vavasour, KP Whittall, AL MacKay, DK Li, G Vorobeychik, and DW Paty. A comparison between magnetization transfer ratios and myelin water percentages in normals and multiple sclerosis patients. Magn Reson Med, 40(5):763-768, 1998. [22] GJ Stanisz, A Kecojevic, MJ Bronskill, and RM Henkelman. Characterizing white matter with magnetization transfer and T2. Magn Reson Med, 42(6): 1128-1136, 1999. [23] P:J. Gareau, B.K. Rutt, S.J. Karlik, and J.R. Mitchell. Magnetization transfer and multi-component T2 relaxation measurements with histopathologic correlation in an experimental model of MS. J. Magn. Reson. Imag., ll(6):585-595, 2000. [24] G.J. Stanisz and R.M. Henkelman. Diffusional anisotropy of T2 components in bovine optic nerve. Magn. Reson. Med., 40:405-410, 1998. [25] R. S. Menon and P. S. Allen. Application of continuous relaxation time distributions to the fitting of data from model systems and excised tissue. Magn. Reson. Med., 20:214-227, 1991. [26] D. L. Phillips. A technique for the numerical solution of certain integral equations of the first kind. Journal of the ACM, 9:84-97, 1962. [27] P. C. Hansen. Numerical tools for analysis and solution of Fredholm integral equations of the first kind. Inverse Problems, 8:849-872, 1992. [28] C. L. Lawson and R. J. Hanson. Solving least square problems. Prentice Hall, Englewood Cliffs NJ, 1974. [29] DB Twieg. The k-trajectory formulation of the NMR imaging process with applications in analysis and synthesis of imaging methods. Med. Phys., 10:610-621, 1983. [30] M.H. Levitt and R. Freeman. Compensation for pulse imperfections in NMR spin-echo experiments. J. Magn. Reson., 43:65-80, 1981. [31] CS. Poon and R.M. Henkelman. Practical T2 quantitation for clinical applications. J. Magn. Reson. Imag., 2:541-553, 1992. [32] R. M. Henkelman. Measurement of signal intensities in the presence of noise in MRI images. Med. Phys., 12(2):232-233, 1985. 156 Bibliography [33] E. W. Weisstein. Rice distribution. Eric Weisstein's World of Mathematics, 2000. http://mathworld.wolfram.com/RiceDistribution.html. [34] James A. Roberts. http://people.eecs.ku.edu/~jroberts/private/Appendix_A.pdf. [35] A.P. Crawley and R.M. Henkelman. Errors in T2 estimation using multislice multiple-echo imaging. Magn. Reson. Med., 4:34-47, 1987. [36] I. M. Vavasour, K. P. Whittall, D. K. Li, and A. L. MacKay. Different magnetization transfer effects exhibited by the short and long T2 components in human brain. Magn Reson Med, 44(6):860-6, December 2000. [37] D.E. Woessner. Effects of diffusion in nuclear magnetic resonance spin-echo experiments. J. Chem. Phys., 34(6):2057-61, 1961. [38] J.H. Jensen, R. Chandra, and H. Yu. Quantitative model for the interecho time dependence ofthe CPMG relaxation rate in iron-rich gray matter. Magn. Reson. Med., 46:159-165, 2001. [39] CL. Lawson and R.J. Hanson. Solving least square problems. Prentice Hall, Englewood Cliffs NJ, 1974. [40] K.P. Whittall and A.L. MacKay. Quantitative interpretation of NMR relaxation data. J. Magn. Reson., 84:134-152, 1989. [41] K. P. Whittall, M. J. Bronskill, and R. M. Henkelman. Investigation of analysis techniques for complicated NMR relaxation data. J. Magn. Reson., 95:221-234, 1991. [42] R. M. Kroeker and R. M. Henkelman. Analysis of biological NMR relaxation data with continuous distributions of relaxation times. J. Magn. Reson., 69:218-235, 1986. [43] K. P. Whittall and A. L. MacKay. Quantitative analysis of NMR relaxation data. J. Magn. Reson., 84:134-152, 1989. [44] S. J. Graham, P. L. Sanchev, and M. J. Bronskill. Criteria for analysis of multicomponent tissue T2 relaxation data. Magn. Reson. Med., 35(3):370-378, 1996. 157 Bibliography [45] K. P. Whittall, A. L. MacKay, and D. K. B. Li. Are mono-exponential fits to a few echoes suf ficient to determine T2 relaxation for in vivo human brain? Magn. Reson. Med., 41 (6): 1255-1257, 1999. [46] K. A. Kraft, P. P. Fatouros, G. D. Clarke, and P. R. S. Kishore. An MRI phantom material for quantitative relaxometry. Magn. Reson. Med., 5:555-562, 1987. [47] J. O. Christoffersson, L. E. Olsson, and S. Sjoberg. Nickel-doped agarose gel phantoms in mr imaging. Acta Rad., 32:426-431, 1991. [48] M. D. Mitchell, H. L. Kundel, L. Axel, and P. M. Joseph. Agarose as a tissue equivalent phantom material for NMR imaging. Magn. Reson. Imag., 4:263-266, 1986. [49] A.L. MacKay, K.P. Whittall, J. Adler, D.K.B. Li, D.W. Paty, and D. Graeb. In vivo visu alization of myelin water in brain by magnetic resonance. Magn. Reson. Med., 31:673-677, 1994. [50] P. Perona and J. Malik. Scale space and edge detection using anisotropic diusion. Proc. IEEE Computer Society Workshop on Computer Vision, 1987. [51] G. Gerig, O. Kubler, R. Kikinis, and F. Jolesz. Nonlinear anisotropic filtering of MRI data. IEEE Transactions on Medical Imaging, 11(2):221-231, 1992. [52] J. R. Mitchell, S. J. Karlik, D. H. Lee, M. Eliasziw, G. P. Rice, and A. Fenster. Quantification of multiple sclerosis lesion volumes in 1.5 and 0.5 T anisotropically filtered and unfiltered MR exams. Med Phys, 23(l):115-26, 1996. [53] B. Mackiewich. Intracranial boundary detection and radio frequency correction in magnetic resonance images. Master's thesis, Simon Fraser University, 1995. [54] M. M. Orkisz, C. Bresson, I. E. Magnin, O. Champin, and P. C. Douek. Improved ves sel visualization in MR angiography by nonlinear anisotropic filtering. Magn Reson Med, 37(6):914-9, 1997. [55] R. Zingale and A. Zingale. Detection of MRI brain contour using isotropic and anisotropic diffusion filter. A comparative study. J Neurosurg Sci, 42(2):lll-4, 1998. 158 Bibliography [56] S.M. Smith and J.M. Brady. SUSAN - a new approach to low level image processing. Int. J. Comput. Vis., 23:34ff, 1997. [57] L. Woog. Adapative waveform algorithms for denoising. PhD thesis, Yale University, 1996. [58] R.R. Coifman, M.V. Wickerhauser, and L. Woog. Adaptive design tookit software. Fast Mathematical Algorithms and Software Corporation, 1997. [59] J.C. Wood and K.M. Johnson. Wavelet packet denoising of magnetic resonance images: Importance of rician noise at low SNR. Magn Reson Med, 41:631-635, 1999. [60] Ingrid Daubechies. Ten lectures on wavelets. Society for Industrial and Applied Mathematics, 1992. [61] H. Gudbjartsson and S. Patz. The Rician distribution of noisy MRI data. Magn. Reson. Med., 34:910-914, 1995. [62] B. Efron and R. Tibshirani. Bootstrap methods for standard errors, confidence intervals, and other measures of statistical accuracy. Statistical Science, l(l):54-77, 1986. [63] A. Macovski. Selective projection imaging: Applications to radiography and NMR. IEEE Trans, on Med. Imag., l(l):42-47, July 1982. [64] T. Brosnan, G. Wright, D. Nishimura, Q. Cao, A. Macovski, and F.G. Sommer. Noise reduction in magnetic resonance imaging. Mag Reson Med, 8:394-409, 1988. [65] K.P. Whittall, M.J. Bronskill, and R.M. Henkelman. Investigation of analysis techniques for complicated NMR relaxation data. J Magn Reson, 95:221-234, 1991. [66] G.E. Backus and J.F. Gilbert. The resolving power of growth earth data. Geophys. J. Roy. Astron. Soc, 16:169-205, 1968. [67] G.E. Backus and J.F. Gilbert. Uniqueness in the inversion of inaccurate gross earth data. Phil. Trans. Roy. Soc. London Ser. A, 266(123-192), 1970. [68] David Heeger. Linear Systems, http://white.stanford.edu/~heeger/linear-systems/linear-systems.html. 159 Bibliography [69] J.M. Bland and D.G. Altman. Statistical method for assessing agreement between two meth ods of clinical measurement. The Lancet, i:307—310, 1986. [70] J.M. Bland and D.G. Altman. Measuring agreement in method comparison studies. Statistical Methods in Medical Research, 8:135-160, 1999. [71] S. Majumdar, S.C. Orphanoudakis, A. Gmitro, M. O'Donnell, and J.C. Gore. Errors in the measurements of T2 using multiple-echo MRI techniques I. Effects of radiofrequency pulse imperfections. Magn. Reson. Med., 3:397-417, 1986. [72] G.P. Liney, A.J. Knowles, D.J. Manton, and A. Horsman. T2 mapping with the fast spin-echo and conventional spin-echo sequences in the presence of BI inhomogeneities. In Proceedings of the ISMRM, page 2071, 1997. [73] Zur Y. and Stokar S. A phase-cycling technique for cancellling spurious echoes in NMR imaging. J. Magn. Reson., 71:212-228, 1987. [74] J.H. Duijun, J.H.N Creyghton, and J. Smidt. Suppression of artifacts due to imperfect TT pulses in multiple echo Fourier imaging. In Proceedings of the Third Annual Society for Magnetic Resonance, number 197-198, 1984. [75] J.G. Sled and G.B. Pike. Correction for BI and BO variations in quantitative T2 measure ments using MRI. Magn. Reson. Med., 43:589-593, 2000. [76] M. Lepage, P.S. Tofts, S.A.J. Back, P.M. Jayasekera, and C. Baldock. Simple methods for the correction of T2 maps of phantoms. Magn. Reson. Med., 46:1123-9, 2001. [77] R. Kaiser, E. Bartholdi, and R.R. Ernst. Diffusion and field-gradient effects in NMR fourier spectroscopy. J. Chem. Phys., 60(8):2966-2979, 15 April 1974. [78] J. Hennig. Multiecho imaging sequences with low refocusing flip angles. J. Magn. Reson., 78:397-407, 1988. [79] J Simbrunner. Generalization of the partition method for calculating echo magnitudes. J. Magn. Reson. A, 109:117-120, July 1994. 160 Bibliography [80] C. Seong, D. Jones, W.E. Reddick, R.J. Ogg, and R.G. Steen. Establishing norms for age-related changes in proton Tl of human brain tissue in vivo. Magn Reson Imag, 15(10):1133-1143, 1997. [81] R Stollberger and P Wach. Imaging of the active Bl field in vivo. Magn Reson Med, 35(2):246-51, 1996. [82] E. T. Lebsack and S. M. Wright. Iterative RF pulse refinement for magnetic resonance imaging. IEEE Trans Med Biol Eng, 49(l):41-48, 2002. [83] M. Alecci, C. M. Collins, M. B. Smith, and P. Jezzard. Radio frequency magnetic field mapping of a 3 Tesla birdcage coil: Experimental and theoretical dependence on sample properties. Magn Reson Med, 46(2):379-385, 2001. [84] D.L Hoult. Signal and Power in High Field Engineering. In ISMRM Workshop on High Field Engineering, pages 2-3, 14 October 1999. [85] P.J. Gareau, B.K. Rutt, CV. Bowen, S.J. Karlik, and J.R. Mitchell. In vivo measure ments of multi-component T2 relaxation behaviour in guinea pig brain. Magn. Reson. Imag., 17(9):1319-1325, 1999. [86] S. Meiboom and D. Gill. Modified spin-echo method for measuring nuclear relaxation times. Phys. Rev., 29:688-691, 1958. [87] H.Y. Carr and E.M. Purcell. Effects of diffusion on free precession in nuclear magnetic resonance experiments. Phys. Rev., 94:630-639, 1954. 161 Appendix A Related Publications • Chapter 2: Myelin Water Fraction Noise Reduction - Paper: Jones, C.K., Whittall, K.P., Mackay, A.L., "Robust Myelin Water Quantifi cation: Averaging vs Spatial Filtering", Magn Reson Med (in press), 2003. - Presentation: Jones, C.K., Whittall, K.P., MacKay, A.L., "Validation of Anisotropic Filtering for Myelin Water Maps". Proceedings of the International Society of Magnetic Resonance in Medicine. Glasgow, Scotland, April 21-27, (2001). Page 820. • Chapter 3: Myelin Water Fraction from Linear Combination - Paper: Jones, C.K., Xiang, Q.S., Whittall, K.P., Mackay, A.L., "Short T2 Selection From Linear Combination of Images", Magn Reson Med (accepted), 2003. - Presentation: Jones, C.K., Xiang, Q.S., Whittall, K.P., Mackay, A.L., "Multi-Echo Linear Combination for Myelin Water Imaging". Proceedings of the International Society of Magnetic Resonance in Medicine. Honolulu, HI, May 18-24, (2002). Page 2289 • Chapter 4: Non-180 Refocusing Pulses - Paper: Jones, C.K., Xiang, Q.S., Whittall, K.P., Mackay, A.L., "Calculating T2 and B\ from Decay Curves Collected with non-1800 Refocusing Pulses", Magn Reson Med (submitted), 2003. - Presentation: Jones, C.K., Xiang, Q.S., Whittall, K.P., Mackay, A.L., "Calculating T2 and B\ from Decay Curves Collected with non-1800 Refocusing Pulses", Proceedings of the International Society of Magnetic Resonance in Medicine. Toronto, ON, July 10-16, (2003). 162
- Library Home /
- Search Collections /
- Open Collections /
- Browse Collections /
- UBC Theses and Dissertations /
- T2 decay curve acquisition and analysis in MRI noise...
Open Collections
UBC Theses and Dissertations
Featured Collection
UBC Theses and Dissertations
T2 decay curve acquisition and analysis in MRI noise considerations, short T2, and B1 field encoding Jones, Craig K. 2003-11-13
pdf
Page Metadata
Item Metadata
Title | T2 decay curve acquisition and analysis in MRI noise considerations, short T2, and B1 field encoding |
Creator |
Jones, Craig K. |
Date Issued | 2003 |
Description | The myelin sheath is a lipid bilayer wrapped around axons in the human brain. The myelin allows faster conduction and uses less energy than along non-myelinated axons. Diseases, such as multiple sclerosis, break down the myelin sheath resulting in cognitive and physical disability. Water trapped between the myelin bilayers has a T2 ≈ 15 ms compared to intra/extra-cellular water with a T2 ≈ 80 ms and cerebrospinal fluid which has a T2 ≈ 2000 ms. A multi-echo MRI pulse sequence is used to acquire a T2 decay curve which is fit using a non-negative least-squares (NNLS) curve fitting algorithm to calculate the signal as a function of the T2 (T2 distribution). The myelin water fraction (MWF) is calculated as the signal with a T2 < 50 ms divided by the total signal in the T2 distribution. Quantification of the MWF in vivo, is slow, prone to errors due to artifacts in the T2 decay curve, and requires many data acquisition averages to attain the necessary high signal-to-noise ratio. Each of these areas were addressed in this thesis. First, I compared the accuracy and precision of MWF estimates from a set of simulated and in vivo multi-echo data with and without noise reduction filters applied. The MWF estimated from anisotropically filtered multi-echo data had the smallest variability, over homogeneous regions, compared to MWF estimates from other filtered data. Second, I created a new technique to optimize coefficients that when linearly combined with multi-echo data, result in fast, accurate estimates of the MWF. Simulations and in vivo estimates of the MWF were as accurate as those from the NNLS algorithm and 20,000 times faster. Finally, the standard multi-echo sequence was optimized to remove artifacts due to inaccurate refocusing pulses. I created a novel T2 curve fitting algorithm, based on NNLS, to accurately estimate the T 2 distribution and refocusing pulse flip angle from a T2 decay curve. Simulations and data acquired on phantoms showed the new technique was as accurate as the NNLS algorithm in quantifying the T2 distribution parameters. In vivo myelin water maps were as good as those estimated from the optimized multi-echo pulse sequence. |
Extent | 8349221 bytes |
Genre |
Thesis/Dissertation |
Type |
Text |
File Format | application/pdf |
Language | eng |
Date Available | 2009-11-13 |
Provider | Vancouver : University of British Columbia Library |
Rights | For non-commercial purposes only, such as research, private study and education. Additional conditions apply, see Terms of Use https://open.library.ubc.ca/terms_of_use. |
DOI | 10.14288/1.0085721 |
URI | http://hdl.handle.net/2429/14941 |
Degree |
Doctor of Philosophy - PhD |
Program |
Physics |
Affiliation |
Science, Faculty of Physics and Astronomy, Department of |
Degree Grantor | University of British Columbia |
Graduation Date | 2003-05 |
Campus |
UBCV |
Scholarly Level | Graduate |
Aggregated Source Repository | DSpace |
Download
- Media
- 831-ubc_2003-859312.pdf [ 7.96MB ]
- Metadata
- JSON: 831-1.0085721.json
- JSON-LD: 831-1.0085721-ld.json
- RDF/XML (Pretty): 831-1.0085721-rdf.xml
- RDF/JSON: 831-1.0085721-rdf.json
- Turtle: 831-1.0085721-turtle.txt
- N-Triples: 831-1.0085721-rdf-ntriples.txt
- Original Record: 831-1.0085721-source.json
- Full Text
- 831-1.0085721-fulltext.txt
- Citation
- 831-1.0085721.ris
Full Text
Cite
Citation Scheme:
Usage Statistics
Share
Embed
Customize your widget with the following options, then copy and paste the code below into the HTML
of your page to embed this item in your website.
<div id="ubcOpenCollectionsWidgetDisplay">
<script id="ubcOpenCollectionsWidget"
src="{[{embed.src}]}"
data-item="{[{embed.item}]}"
data-collection="{[{embed.collection}]}"
data-metadata="{[{embed.showMetadata}]}"
data-width="{[{embed.width}]}"
async >
</script>
</div>
Our image viewer uses the IIIF 2.0 standard.
To load this item in other compatible viewers, use this url:
http://iiif.library.ubc.ca/presentation/dsp.831.1-0085721/manifest