UBC Theses and Dissertations

UBC Theses Logo

UBC Theses and Dissertations

T2 decay curve acquisition and analysis in MRI noise considerations, short T2, and B1 field encoding 2003

You don't seem to have a PDF reader installed, try download the pdf

Item Metadata

Download

Media
ubc_2003-859312.pdf
ubc_2003-859312.pdf [ 7.96MB ]
ubc_2003-859312.pdf
Metadata
JSON: 1.0085721.json
JSON-LD: 1.0085721+ld.json
RDF/XML (Pretty): 1.0085721.xml
RDF/JSON: 1.0085721+rdf.json
Turtle: 1.0085721+rdf-turtle.txt
N-Triples: 1.0085721+rdf-ntriples.txt
Citation
1.0085721.ris

Full Text

T2 Decay Curve Acquisi t ion and Analysis in M R I Noise Considerations, Short T 2 , and B x F ie ld Encoding by Craig K . Jones M . S c , The University of Western Ontario, 1997 B . S c , Simon Fraser University, 1992 A THESIS S U B M I T T E D IN P A R T I A L F U L F I L M E N T OF T H E R E Q U I R E M E N T S F O R T H E D E G R E E O F D O C T O R OF P H I L O S O P H Y in The Faculty of Graduate Studies (Department of Physics and Astronomy) We accept this thesis as conforming to the required standard T H E U N I V E R S I T Y OF BRITISH C O L U M B I A 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 M R I 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 M W F 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 M W F estimates from a set of simulated and in vivo multi-echo data with and without noise reduction filters applied. The M W F estimated from anisotropically filtered multi-echo data had the smallest variability, over homogeneous regions, compared to M W F 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 M W F . Simulations and in vivo estimates of the M W F were as accurate as those from the N N L S 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 N N L S , 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 N N L S 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. i i Table of Contents Abstract i i Table of Contents i i i List of Tables vii List of Figures ix List of Algorithms x i Acknowledgements x i i Glossary xiv 1 Introduction 1 1.1 Magnetic Resonance Imaging 1 1.2 T 2 and T i 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 T 2 Decay Curve 7 ' 1.5.1 Mono-Exponential T 2 Relaxation 7 .1.5.2 Multi-Exponential T 2 Relaxation . . 8 1.6 T 2 Decay Curve Acquisition 9 1.6.1 K-Space . . . 9 1.6.2 Optimized Multi-Echo Pulse Sequence 10 1.6.3 Standard Multi-Echo Pulse Sequence 11 1.6.4 Noise in Decay Curve 14 1.6.5 Artifacts Influencing the T 2 Decay Curve 16 1.7 T 2 Decay Curve Analysis 18 1.7.1 Decay Curve Inversion (Decay Curve to T 2 Distribution) 18 1.7.2 T 2 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 27 2.2 SNR Effect on Myelin Water Fraction 28 2.3 Description of Spatial Filters 29 2.3.1 Anisotropic Diffusion Filter 30 2.3.2 Susan Filter 31 2.3.3 Median Filter 31 2.3.4 Wavelet Filter 32 2.4 Spatial Filters on Simulated Data . . . 32 2.4.1 Methods 32 2.4.2 Results 36 2.4.3 Conclusions 40 2.5 Spatial Filters Applied to Zn Vivo Data , . 4 0 2.5.1 Data Acquisition and Analysis 40 2.5.2 Results 41 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 Conclusions 53 3 Myelin Water Fraction from Linear Combination 55 3.1 Introduction 55 3.1.1 Linear Combination for T2 Selectivity 56 3.1.2 Properties of a Linear Combination . . 57 3.1.3 Overview 57 3.2 Calculation of Coefficients 58 3.2.1 Coefficients for Short T 2 59 3.2.2 Coefficients for Al l T 2 60 3.2.3 Description of the Algorithm '. . . . 61 3.2.4 Results 62 3.2.5 Discussion 63 3.3 Linear Combination Simulations 64 3.3.1 Model 64 3.3.2 Results and Discussion 64 3.3.3 Conclusions 65 3.4 Selectivity of Mono-Exponential Phantoms 66 3.4.1 Methods 66 3.4.2 Results and Discussion 66 3.4.3 Conclusions 68 3.5 In Vivo Myelin Water Fraction: NNLS vs. T 2 Filter 68 3.5.1 Methods 69 3.5.2 Scan Results 69 3.5.3 Conclusions 72 iv Table of Contents 3.6 Calculation of Coefficients by Matrix Inversion 72 3.6.1 Calculation of Coefficients 72 3.6.2 Results and Discussion 74 3.6.3 Conclusions 76 3.7 Calculation of Coefficients by Backus-Gilbert 76 3.7.1 Calculation of Coefficients : 76 3.7.2 Results and Discussion 77 3.7.3 Conclusions 80 3.8 Optimal Echo Times . 80 3.8.1 Filter Measures 80 3.8.2 Equally Spaced Echoes 82 3.8.3 Non-Equidistant . . . . : 85 3.8.4 Discussion. . . 88 3.8.5 Four Echo Non-equidistant Sampling . 90 3.9 Conclusions • • • • . - 9 2 4 Non-180 Refocus ing Pulses 94 4.1 Introduction . . 94 4.1.1 Motivation 94 4.1.2 Previous Work 95 4.1.3 Implementation of Algorithm 98 4.1.4 Overview of Work 102 4.2 Mono-Exponential T 2 Decay Curve Fi t 102 4.2.1 Mono-Exponential Fit Algorithm 102 4.2.2 Monoexponential Fi t - Noise Simulation 103 4.2.3 Phantom Verification 105 4.3 Multi-Exponential T 2 Decay Curve Fi t : ,111 4.3.1 Multi-Exponential Fit Algorithm :. . . . ' T i l 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 126 4.4.2 Methods 127 4.4.3 Results and Discussion '. . . . 127 4.4.4 Conclusions 135 4.5 Non 90° Excitation Pulse - - 136 4.5.1 Decay Curve Calculation 136 4.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 F l ip Angles 141 4.6.1 Simulations 143 4.6.2 Scanning 144 4.6.3 Conclusion 149 4.7 Conclusions 149 5 Conclusions 150 5.1 Filtering 150 5.2 Linear Combination 150 5.3 B ! / T 2 152 5.4 Preferred T 2 Quantification Protocol 153 Bibliography . . . . . . . . . 154 A Related Publications • 162 List of Tables 1.1 T i and T 2 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 37 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 ; 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 T 2 and all T 2 linear combination method 62 3.2 Mean and standard deviation of the signal in the short T 2 phantom image . . . . . . 67 3.3 Coefficients for the short, medium and long T 2 selection based on matrix inversion. 74 3.4 Coefficients calculated by the Backus-Gilbert method 77 3.5 Coefficients for short and all T 2 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 4.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 T 2 time calculated from the M E S E Q sequence and the M E S E C (with the flip angle calculated) . 1 1 9 4.7 Refocusing pulse flip angle calculated from multi-echo decay curve •'. 119 4.8 Estimated T i from curve fits of five points (TR = 100,300,800,2000,3000 ms) from a pulse saturation experiment . 1 1 9 4.9 Percent myelin water estimated from scans with (xnom = 150° and a n o m = 180°. . 121 4.10 Refocusing pulse flip angle calculated over the regions from the scans with <xnom — 150° and a n o m = 180° . 1 2 3 4.11 Geometric mean T 2 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°, 8 5 ° , . . . ,50° and ar = 2ae 139 4.15 The excitation and refocusing pulses were varied with ae = 90°, 8 5 ° , . . . , 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 T 2 and refocusing pulse flip angle (for the first echo) calculated from each train of refocusing pulses 147 4.18 The resulting x2 misfit given the scaling of the R F pulses was 1.0,1.0,1.0, S90 from s<4) 148 vi i i 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 M R I image and corresponding k-space ' 10 1.5 Optimized single-slice, multi-echo pulse sequence ( M E S E 0 ) 11 1.6 Conventional multiple echo 2D spin-echo pulse sequence (MESE C ) . . . . . . . . . . 12 1.7 Conventional multiple echo 3D spin-echo pulse sequence (MESE C ) . . . . . . . . . . 13 .: 1.8 Rice distribution as a function of the mean . .14 1.9 Decay and T 2 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 M R I data 43 2.6 Myelin water maps from unfiltered and filtered (anisotropic diffusion filter) multi- echo data . 4 9 2.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" T 2 filter 61 3.3 T2 filters for short T 2 components and all T 2 components 63 3.4 Myelin water fraction noise simulation to compare linear combination and N N L S . . 65 3.5 Linear combination short T 2 selectivity of phantoms 67 3.6 Mean signal from phantoms in selected image versus T 2 filter 68 3.7 In vivo myelin water maps calculated by linear combination and N N L S 70 3.8 Myelin water fractions calculated per region from the linear combination method and N N L S 71 3.9 Bland-Altman plots for white matter regions and gray matter regions 71 3.10 T2 filters based on coefficients calculated by matrix inversion 73 3.11 Brain images T 2 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 T E 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 T E 86 3.18 Width of T 2 filter for unequally spaced echo times 87 3.19 The T 2 filter, exp ( - T E / T 2 ) as a function of T E 89 3.20 The short T 2 filter as a function of a consecutive subset of echoes based on the thirty-two coefficients 89 3.21 Short T 2 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 echoes 92 4.1 Phase coherence diagram 96 4.2 Decay curves, collected with ap — 180°, fit to obtain T 2 and a 107 4.3 Decay curves, collected with ap = 150°, fit to obtain T 2 and a 108 4.4 Decay curves, collected with ap = 110°, fit to obtain T 2 and a 109 4.5 Decay curves, collected with ap = 90°, fit to obtain T 2 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 T 2 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 T E = 10 ms image from the four average, 32 echo M E S E G sequence and six slices of the six slice slab M E S E C acquisition 129 4.15, Myelin water maps calculated from a 6-slice slab M E S E C . . . . . . . 130 4.16 Difference between the myelin water map calculated from the M E S E C sequence and the myelin water map calculated from the M E S E C sequence 131 4.17 In vivo decay curve from M E S E C and M E S E 0 from the genu ofthe corpus callosum (WM) 134 4.18 Maps ofthe refocusing pulse flip angle calculated from a 6-slice slab M E S E C . . . . 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 T 2 = 85 ms 145 4.22 Decay curve fit of decay curves collected with non-constant refocusing pulse train, phantom with T 2 = 23 ms 146 4.23 Measured and fitted data from non-constant train of refocusing pulses 148 x List of Algorithms 1.1 Small norm N N L S regularization 23 4.1 Signal calculation from a train of refocusing pulses . 100 4.2. T i relaxation 101 4.3 T 2 relaxation . . . : . . 101 4.4 Phase evolution of magnetization . 101 xi Acknowledgements When I first moved to Vancouver in 1997, I began working for Dr. Don Paty, Dr. David L i , and Andrew Riddehough at the M S / M R I 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. L i , were overwhelmingly supportive. Dr. Paty's comment was, "I was wondering when you would come and talk to me about this." Wi th a lot of work from Andrew, an arrangement was worked out such that I worked in the M S / M R I 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 M S / M R I 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 t c 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 M R I 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. M y 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... xi i i Glossary Abbreviation Meaning B G Backus-Gilbert CSF Cerebrospinal Fluid F O V Field of View (cm) F W H M Full-Width-at-Half-Max G M Gray Matter ' M E " " • - •  Multi-echo M E S E Multi-echo Spin-Echo Pulse Sequence M R I Magnetic Resonance Imaging MS Multiple Sclerosis M T Magnetization Transfer M W F Myelin Water Fraction N N L S Non-Negative Least Squares R F Radiofrequency ROI Region-of-Interest R M S Root-Mean-Square S A R Specific Absorption Rate SI Signal Intensity SNR Signal to Noise Ratio T Echo spacing (ms) T E Echo Time (ms) T i Spin-Lattice Relaxation Time (ms) T R Repetition Time (ms) T 2 Spin-Spin Relaxation Time (ms) xiv Glossary Glossary W M White Matter xv Chapter 1 Introduction In magnetic resonance imaging, information about the local water environment is quantified by T i and T 2 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 T 2 , 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 T 2 to assess tissue and pathology (Section 1.4). The T 2 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 R F signal from the protons. The spatially encoded R F 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 (T 2 ) and the spin-lattice relaxation time (Ti) . 1 C H A P T E R 1. I N T R O D U C T I O N 1.2. T 2 AND Ti RELAXATION 1.2 T 2 and T i Relaxation The T 2 relaxation time calculated from a multi-echo decay curve acquired using the M R I 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 T i 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 M R I voxel, typically contains several distinct water compartments. Each of the water compartments may have a different T 2 , therefore, a modified set of magnetization' equations must be analyzed to calculate the T2 of each compartment. 1.2.1 Spectral Density The T i and T2 relaxation times are related to the Larmor frequency of the proton, wo, and the correlation time, T c , 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/T 2 = k \jM(0) + 5-jU("o) + \jW(2u0) (1.3) where (w) is defined in Equation 1.1 and k = (^f 7 4 / i 2 § J (I + 1) [2]. 2 C H A P T E R 1. I N T R O D U C T I O N 1.2. T 2 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 ) x B ] - ^ - ^ i y + T l * ( L 4 ) where M ( i ) is the magnetization vector; B is the applied magnetic field; T 2 is the spin-spin relaxation time; T i is the spin-lattice relaxation time; 7 = 42.577 M H z / T is the gyromagnetic ratio of a proton; M o 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) = 7 M , 5 0 - ^ (1.5) jMy{t) = l M x B 0 - ^ (1.6) ±MM (,7) The 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) = M 0 + [ M z ( 0 ) - M 0 ] e x p ( - i / T 1 ) . (1.10) where 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) = M 0 e ^ 2 e1Wot (1.11) Mz(t) = M 0 ( l - e - ' / T l ) . (1.12) C H A P T E R 1. I N T R O D U C T I O N 1.3. CENTRAL NERVOUS SYSTEM where Mxy(0) = 1 and M 2 (0 ) = 0. The frequency of precession OJQ is typically demodulated out of the transverse magnetization to obtain: Mxy(t) = M 0 e - * / T 2 (1.13) Mz(t) = M0 ( l - e - * / T i ) . (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 R F coil is digitized and filtered to obtain the signal intensity. The primary interest, for this thesis, is T 2 relaxation, therefore, the acquisition and analysis of T 2 decay curves wil l 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. M R I is well suited for this and therefore, will continue to be an important factor in assessment, research and treatment of disease. 4 C H A P T E R 1. I N T R O D U C T I O N 1.3. C E N T R A L 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, T 2 , which is a function of the local water environment. The T 2 decay curve acquisition, from a stan- dard M R I pulse sequence, results in signal which decays as an exponential function of T 2 . Simple systems decay as a function of a single T 2 (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 M R I (ss 1 mm in-plane resolution) and therefore, the signal decays as a sum of single T 2 exponentials (multi-exponential decay). 5 C H A P T E R 1. I N T R O D U C T I O N 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 T 2 decay curve analysis of frog sciatic nerve. They found three distinct T 2 water compartments based on a "peeling-off" mathematical procedure applied to multi-echo data. A slowly relaxing compartment, T 2 ~ 300 — 500 ms, attributed to water in the intra-cellular space. A n intermediate relaxing component, T 2 ~ 80 ms ascribed to axoplasmic water. And a fast relaxing component, T 2 ~ 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 C O N T I N 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 T 2 peak (T 2 ~ 15 ms) was attributed to water trapped between the myelin bilayers. Menon fit T 2 decay curves, using N N L S , acquired from myelinated cat brain nerves and found C H A P T E R 1. I N T R O D U C T I O N 1.5. T 2 DECAY CURVE four distinct components: component near T 2 = 1 ms attributed to phospholipid protons in the tissue, component near T 2 = 12.7 ms attributed to water in the myelin layers, component near T 2 = 89 ms attributed to intra and extra-cellular water, and a component not attributed to anything at T 2 ~ 340 ms. The short component, T 2 = 12.7 ms, contributed 6.8% to the total water in the T 2 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 T 2 component was by MacKay et al. [15, 16]. In vivo T 2 relaxation measurements [17, 18, 19] distinguish three water compartments in the brain: water trapped between the myelin bilayers (10 < T 2 < 50 ms), intra/extra axonal water (T 2 ss 80 ms) and water associated with cerebrospinal fluid (CSF) (T 2 > 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 T 2 component. Beaulieu et al. [20] compared T 2 distributions from myelinated and non-myelinated nerves from the garfish and found the short T 2 component only in the myelinated axon. More recently, many studies have compared the short T 2 component to other M R I phe- nomenon such as magnetization transfer imaging [21, 22, 23] and diffusion imaging [24]. Each water compartment is "visible" to M R I [17, 18, 19, 25], therefore, M R I 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 T 2 Relaxation The T 2 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 T 2 = 20 ms (circles), T 2 = 80 ms (crosses), and T 2 = 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 C H A P T E R 1. I N T R O D U C T I O N 1.5. T 2 DECAY CURVE 1000r 150 TE (ms) Figure 1.3: Simulated decay curve with signal intensity S(t) = pexp (—£/T 2) with p = 1000, t = 10,20,. . . , 320 ms and T 2 = 20 ms (circles), T 2 = 80 ms (crosses), and T 2 = 2000 ms (diamonds). T R was considered infinite. 1.5.2 Mult i -Exponential T 2 Relaxation Many tissues in the human body have more than one water compartment visible to M R I (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 T 2 = 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 C H A P T E R 1. I N T R O D U C T I O N 1.6. T 2 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 T 2 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 M R I pulse sequence is described by gradients along three orthogonal axes Gx, Gy, and Gz; the R F 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 wil l 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 M R I 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 C H A P T E R 1. I N T R O D U C T I O N 1.6. T 2 DECAY CURVE ACQUISITION Figure 1.4: Standard M R I 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 M R I 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 A n optimized, multiple echo spin-echo pulse sequence (Figure 1.5) was used to acquire multi-echo images from which multiple component T 2 measurements were calculated. The implementation of the pulse sequence used for this thesis had a train of 32 composite refocusing pulses (902-180,,- 10 C H A P T E R 1. I N T R O D U C T I O N 1.6. T 2 DECAY CURVE ACQUISITION 90 x) 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 B i [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 wil l be annotated as MESE„ in this thesis. RF G D A C Q / v 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 ( M E S E 0 ) . 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 G 2 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 C H A P T E R 1. I N T R O D U C T I O N 1.6. T 2 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 G z were added to the constant gradient crushers along G 2 . 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 wil l be denoted by M E S E C in this thesis. , J v y v y v J v y v n (1 - 4 - J l -Jl (1 n (1 n • 1 n f 1 - 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 (MESE C ) . The excitation pulse and refocusing pulses are slice selective. 12 C H A P T E R 1. I N T R O D U C T I O N 1.6. T 2 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 ( M E S E C ) . 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 C H A P T E R 1. I N T R O D U C T I O N 1.6. T 2 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 M R I image, then the noise is Rayleigh distributed. When ///cr is large, for example in the object of an M R I image, then the noise approximates a Gaussian distribution. As /i/cr —>• 0, the Rician distribution becomes a Rayleigh distribution. The Bessel function 14 C H A P T E R 1. I N T R O D U C T I O N 1.6. T 2 DECAY CURVE ACQUISITION 2O(A«) = 1 when JJL = 0 and exp ( - ^ f 2 ) = exp ^ — w h e n |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 P R i c e ( £ ) around a; = fi: x P(x) — exp x exp (J,2 + x 2 2a 2 x 2 + ^ i 2 a exp 2<72 7 v ^ I x 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. A l l simulations, in the following chapters, used Rician noise. The Rician noise was created as ye{U) = yj[y(ti) + e i ] 2 + e?,, where y is the true signal, and e\ and e 2 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 ( i i ) / S N R . 15 C H A P T E R 1. I N T R O D U C T I O N 1.6. T 2 DECAY CURVE ACQUISITION 1.6.5 Ar t i fac ts Inf luencing the T 2 Decay Cu rve Artifacts in the T 2 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 T 2 decay curve are: 1) signal from out-of-slice, 2) magnetization transfer, 3) stimulated echoes, and 4) T 2 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 R F pulse is non-rectangular. The goal of shaped R F 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 R F pulses is the opposite, the desire is to flip all protons across a sample. Therefore, whether shaped R F pulses or hard R F 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 M R I 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 R F 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 R F 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 C H A P T E R 1. I N T R O D U C T I O N 1.6. T 2 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 T 2 decay curves [31]. In the presence of an R F pulse [37], dephasing magnetization will be rotated such that a portion wil l continue dephasing, a portion wil l 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, wil l 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 T 2 Dependence on Inter-Echo Spacing The standard method to measure T 2 in vivo is to use a multi-echo spin-echo pulse sequence where the constant echo spacing is T (TE = 2r). The T 2 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 R 2 (R 2 = 1/T 2 ) as R 2 = R 2 a + A R 2 . The relaxation rate R 2 a was defined independent of microscopic field inhomogeneities, and A R 2 was derived to be a function which was approximated as being proportional to re 3/ 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 A R 2 is a function of T, the measured R 2 would also be a function of r . R 2 is several orders of magnitude smaller than the measurement error for the typical r of 5 ms. Therefore, comparing T 2 from different techniques, requires knowledge of r and the underlying tissue composition. 17 C H A P T E R 1. I N T R O D U C T I O N 1.7. T 2 DECAY CURVE ANALYSIS 1.7 T 2 Decay Curve Analysis Decay curve data acquired by either technique described in Section 1.6 must be fit to determine the underlying T 2 . 1.7.1 Decay Curve Inversion (Decay Curve to T 2 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 T 2 decay curve of y short(U = 200exp•(—t/20), the medium T 2 decay curve of y m ediumW = 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 T 2 times and amplitudes, s, for a given decay curve y. The bottom plot of Figure 1.9 shows the amplitudes, Sj, and T 2 times of the T 2 decay curves in the top plot of Figure 1.9. In this example, s\ = 200 at T 2 ) i = 20 ms and s 2 = 800 at T 2 2 = 80 ms. The plot of s versus T 2 is called the T 2 distribution. Each data inversion technique may assume positions for the spikes in the T 2 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 M i s f i t 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 a n d L2 which measures | |y - y | | 2 = 2~2iLi \Vi ~ Vi?• Another misfit measure, x 2 , 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 C H A P T E R 1. I N T R O D U C T I O N 1.7. T 2 DECAY CURVE ANALYSIS 900 TE (ms) 800r ••. 700 k 600H & 500 - CO c CD ° 400 - o o £ 300 - 200h 1 0 0 " ; ; : 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 T 2 (ms) Figure 1.9: Decay curves (top) and T 2 distribution (bottom) of a bi-exponential white matter model. The short T 2 peak at 20 ms in the T 2 distribution (bottom) corresponds to the fast decaying exponential (exp (—i/20), dashed-dotted line in the top plot). The medium T 2 peak at 80 ms in the T 2 distribution (bottom) corresponds to the slower decaying exponential (exp (—1/80), dash line in the top plot). The two peaks in the T 2 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 C H A P T E R 1. I N T R O D U C T I O N 1.7. T 2 DECAY CURVE ANALYSIS The expected value E(x2) = N and the standard deviation ax2 = y/2N. Solutions with x 2 ^ N reproduce the noise and solutions with x 2 3> N misrepresent the measurement. Intuitively, if .the only difference between the predicted and measured data is noise, then Equation 1.33 is E i L i = £ £ 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 T 2 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 N N L S 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 T 2 . 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 N N L S is used to solve this modified system, it minimizes the x 2 misfit in Equation 1.33. N N L S is guaranteed to minimize Equation 1.34 [39], or the x 2; 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 T 2 = T 2 j at echo times t. Therefore, N N L S minimizes the misfit of a linear combination of mono-exponential decay curves to create multi-exponential decay curve with the minimum x 2 misfit. Regularization The standard N N L S technique can be modified to apply different types of smoothing to the T 2 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 p s - 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 N N L S solution. 20 C H A P T E R 1. I N T R O D U C T I O N 1.7. T 2 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 = 1 0 0 0 1 0 0 0 : 0 0 ••• . . . o . . . o 1 0 0 1 (1.37) and f = 0. 800 700 600 ^ 5 0 0 </> c CD D c o o 400 CL 300 200 100 20 40 60 80 100 T 2 (ms) 120 140 160 180 200 Figure 1.10: T 2 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 C H A P T E R 1. I N T R O D U C T I O N 1.7. T 2 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 N N L S 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 a 1 2 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 T 2 was defined to be 120 logarithmically spaced T 2 = 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 N N L S solution was used and consisted of finding the minimum x2 solution, Xmim then regularizing the T 2 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 x 2 is described in Algorithm 1.1. 1.7.1.2.1 T 2 Distribution Bins The T 2 distribution calculated using N N L S 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 T 2 relaxation time or amplitude, the T 2 distribution is typically binned. A bin is defined with a start T 2 and end T 2 . There are several measures defined for a bin. The most important measure, in particular for the short T 2 , is the fraction of signal in the bin relative to the total signal: fr = (1.39) 22 C H A P T E R 1. I N T R O D U C T I O N 1.7. T 2 DECAY CURVE ANALYSIS Algorithm 1.1 Small norm N N L S regularization. Regularized N N L S 1: Define flow = 1.02 {X2low = fiowXmin} 2: Define f h i g h = 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\ J S between x20W and Xhigh } 6: while x 2 < X20W O R x 2 > xligh d o 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 x 2 vs fj, 12: From linear fit, calculate fi at midpoint between xfow a n d x\ig 13; Solve s = nnls (A, y, Hi) 14: Calculate x 2 from A, y, s 15: if X ? < X L t h e n 16: Set H to be linear spaced between Hi and M(1) 17: else if Xi > xligh t h e n 18: Set H to be linear spaced between /i(end) and Hi 19: else 20: Do Nothing 21: end if 22: end while 23 C H A P T E R 1. I N T R O D U C T I O N 1.7. T 2 DECAY CURVE ANALYSIS The arithmetic mean T 2 of a bin is defined as: T = E T-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 T 2 compartment (related to water trapped between the myelin bilayers, see Section 1.4) was defined as a bin with 10 < T 2 < 50 ms. The medium T 2 bin was defined as 50 < T 2 < 120 ms and is related to the inter and intra-cellular water. 1.7.2 T 2 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 T 2 compo- nents of one order of magnitude apart in a T 2 distribution. The standard method to attain high SNR, in the M E S E G 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 T R = 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 T 2 Decay Curve Sampling The T 2 decay curve must be sampled sufficiently to accurately estimate the T 2 components and amplitudes. Whit tal l 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 T 2 component, underestimated the T 2 , and overestimated the proton density. Multi-exponential fits of the same simulated data sampled at assessed. 24 C H A P T E R 1. I N T R O D U C T I O N 1.8. PHANTOMS 32 echoes (TE = 10, 20 , . . . , 320 ms) accurately reproduced the two T 2 components. Whit tal l et al. compared mono-exponential T 2 solutions from simulated decay curves sampled with different sets of echo times. They found the T 2 and proton density calculated from two different sampling schemes to be statistically different even though the underlying T 2 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 T 2 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 T i and T 2 times (see Table 1.1) similar to the brain [46]. Nickel-chloride (NiCl 2 ) was chosen, as the paramagnetic ion N i is an effective T i 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 T i were measured from a saturation recovery experiment with T R = 100, 300, 800, 2000, 3000 ms. Other parameters were T E = 11.0 ms, F O V = 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 T 2 measurements were done with a 48 echo M E S E 0 pulse sequence with T E = 10,20,. . . , 480 ms. Other parameters were T R = 3000 ms, 2 averages, and F O V = 16 cm. Phantom T 2 (ms) T i (ms) m M N i C l 2 % 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: T i and T 2 of mono-exponential nickel/agarose phantoms. 25 C H A P T E R 1. I N T R O D U C T I O N 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 T 2 decay curves. The thesis focused on three aspects of the noise and artifacts as they relate to T 2 decay curve analysis. First, the curve fitting algorithms used to fit T 2 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 T 2 decay curve is to use a curve fitting algorithm. The standard algorithm is slow and problematic for low SNR multi-echo data. The short T 2 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 T 2 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 T 2 estimation. Chapter 4 discusses a method to estimate the T 2 and refocusing pulse flip angle from a decay curve that has stimulated echoes. The accuracy and precision of the estimation of T 2 and refocusing pulse flip angle were assessed using noise simulations. The phantoms were scanned with the M E S E 0 sequence to validate the theory. Myelin water fractions and maps estimated from the M E S E C decay curves were compared to those estimated from the M E S E 0 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 < T 2 < .50 ms), intra/extra axonal water (T 2 ~ 80 ms) and water associated with cerebrospinal fluid (CSF) (T 2 > 2000 ms). The multiple water compartments in normal white matter result in a multi-exponential M R signal decay as a function of echo time (TE) and the shortest T 2 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 T 2 (T 2 distribution), was calculated voxel-by- voxel from the T 2 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 T 2 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 M R I by Gerig [51]. It is a locally adaptive smoothing filter that incorporates the local image intensity gradient 27 C H A P T E R 2. M W F NOISE R E D U C T I O N 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 M R I volumes and M R I volumes that contain multiple channels (e.g., multi-echo data). In M R I , the anisotropic diffusion filter has been applied to work in quantification of MS lesion volumes [52], intra-cranial boundary detection and R F correction [53], vessel visualization [54] and detection of brain contour [55]. M y hypothesis is that myelin water fractions (MWFs) calculated from spatially filtered spin- echo data should yield better qualitative and quantitative results than M W F s calculated from unfiltered M R I 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 M R I data with different numbers of averages and with and without filtering and is shown in Section 2.6. 2.2 S N R Effect on M y e l i n W a t e r F r a c t i o n 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 wil l 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 T E =10,20, . . . , 320 ms, a short water compartment with p s = 200 and T 2 S = 20 ms, and a medium water compartment with T 2 m = 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 N N L S . 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 C H A P T E R 2. M W F NOISE R E D U C T I O N 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 M R I to obtain higher SNR images, but the cost is a longer scan time. The scan time of the M E S E G 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 i f post-processing of the M R I 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 M R I 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 C H A P T E R 2. M W F NOISE R E D U C T I O N 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, S U S A N 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 M R I 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 < - i 2^ N (2.3) 30 C H A P T E R 2. M W F NOISE R E D U C T I O N 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 + ... + V 7 M ( ^ ) 2 ' 2 (2.4) where J i , . . . ,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) e x P \ ^ — — : . 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 C H A P T E R 2. M W F NOISE R E D U C T I O N 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 log 2(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. A l l 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 M R I data set was created which modeled important aspects of white matter in an M R I 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 C H A P T E R 2. M W F NOISE R E D U C T I O N 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 C H A P T E R 2. M W F NOISE R E D U C T I O N 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 C H A P T E R 2. M W F NOISE R E D U C T I O N 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 C H A P T E R 2. M W F NOISE R E D U C T I O N 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 C H A P T E R 2. M W F NOISE R E D U C T I O N 2.4. SPATIAL FILTERS ON SIMULATED DATA Filter Global Error Measure none 1.964 (0.005) a n i s o l 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 R M S 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) a n i s o l 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 C H A P T E R 2. M W F NOISE R E D U C T I O N 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) a n i s o l 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 C H A P T E R 2. M W F NOISE R E D U C T I O N 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), a n i s o l (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 C H A P T E R 2. M W F NOISE R E D U C T I O N 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 R M S 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 M R I 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 M R I 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 A n optimized multi-echo (ME), spin-echo pulse sequence (descending, alternating crushers [31] and composite 180° pulses) was used to collect 32 echoes at T E = 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 M E 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). A n 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 M R I parameters were F O V = 24 cm, T R = 3,000 ms, slice thickness = 5 mm, 256 x 128 frequency/phase encoding, and pixel size of 0.94 x 0.94 mm. A l l data was collected on a 1.5T G E Signa Horizon system with EchoSpeed gradients (General Electric Medical Systems, Milwaukee WI). Regions were drawn on the T E =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 i l t ) , Median filter applied to M W M from 40 C H A P T E R 2. M W F NOISE R E D U C T I O N 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 f W M 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) a n i s o l 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) medf i l t 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 C H A P T E R 2. M W F NOISE R E D U C T I O N 2.5. SPATIAL FILTERS ON IN VIVO DATA Filter A l l Regions W M Regions none 9 0 a n i s o l 0 0 aniso3 1 0 aniso5 50 35 susan 0 0 medf i l t 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 C H A P T E R 2. M W F NOISE R E D U C T I O N 2.5. SPATIAL FILTERS ON IN VIVO DATA Figure 2.5: Myelin water maps of the unfiltered (top left), susan (top right), a n i s o l (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 C H A P T E R 2. M W F NOISE R E D U C T I O N 2.6. S P A T I A L F I L T E R I N G VS A V E R A G I N G 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 N N L S 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 M R I 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 T E = 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 C H A P T E R 2. M W F NOISE R E D U C T I O N 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, T E = 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 T 2 distribution was calculated by the regularized N N L S algorithm [40]. The M W F was defined as: where s (T 2 ) was the amplitude of the T 2 distribution calculated by N N L S and the maximum T 2 was 2000 ms. Two sets of analysis were performed: 1) M W F s calculated per ROI and 2) M W F s 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 N N L S algorithm to the averaged decay curve then calculate the M W F (average-invert) . 2. Apply the N N L S algorithm to each decay curve in the ROI, calculate the M W F from each result then average the M W F s ( invert-average). 3. Apply the bootstrap algorithm as defined below (bootstrap). The standard M W F measure was defined to be the average- inver t , 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 T 2 distributions calculated voxel-by-voxel by the N N L S 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 M W F less than 0.03. A l l data was processed using in-house software and Matlab (The Mathworks, Natick, Mass.) (2.11) 45 C H A P T E R 2. M W F NOISE R E D U C T I O N 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 M W F s 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 M W F found. The bootstrap resampling is repeated 1,000 times and the mean and standard • deviation of the corresponding M W F s 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 M W F 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 M E data, I compared the mean and variance of the intensities of the anisotropically filtered M E data in the 40 ROIs to the original unfiltered M E data. The maximum SI change in the M E data on the T E = 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 M W F from the bootstrap method, /(,, to the "gold 46 C H A P T E R 2. M W F NOISE R E D U C T I O N 2.6. SPATIAL FILTERING VS AVERAGING standard" average- inver t , 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 M W F s calculated from the bootstrap method were found to be statistically the same as those calculated by the average- inver t method (minimum p=0.11). For the rest of the work the mean M W F within an ROI was calculated using the boots t rap method unless otherwise stated. 2.6.2.3 Myelin Water Fraction 2.6.2.3.1 Four Average Data The standard M W F measure ( M W F calculated from the 4-average, unfiltered M E data) was compared to the M W F calculated from the 4-average, filtered M E data to determine if the M W F measurement had less variability after the application of the filter to the 4-average M E data. The average relative change in mean M W F for the ROIs, calculated with the boots t rap 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 M W F 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 M W F 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 M W F . The mean of the magnitude of the relative decrease in the standard deviation of the M W F 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 M E data can be collected in half the scan time of the M E data required for the standard (4-averages); therefore, if the M W F s calculated from the 2-average, filtered M E data are as robust as those calculated from the standard, then the 2-average, filtered M E data could be used for future work (which would be a significant savings in scan time). The mean relative difference between the M W F calculated from the 2-average (with and without filtering) multi-echo data and the standard M W F over the 40 ROIs was 26%. Only one 2-average, 5 iteration ROI had a mean M W F statistically the same as the standard. The myelin water map was created by applying N N L S voxel-by-voxel to calculate the myelin 47 C H A P T E R 2. M W F NOISE R E D U C T I O N 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 M W F between the invert-average method (based on the 4-average M W F data) and the standard M W F 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 M W F (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 M R I 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 M W F 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 C H A P T E R 2. M W F NOISE R E D U C T I O N 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 C H A P T E R 2. M W F NOISE R E D U C T I O N 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 M E SI was introduced by the anisotropic diffusion filter (a bias could result in artificial changes in the amplitudes or T 2 times of the peaks in the calculated T 2 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 M E 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 M E 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 M W F results. The M W F calculated from the 4-average, filtered M E data was expected to be more homo- geneous than the M W F calculated from the 4-average M E data (standard) because the local variability in the SI of the M E data decreased as more iterations of the filter were applied. As expected, there was only a small change (< 5%) in the mean M W F for all the ROIs. There was no discernible bias in the M W F with filtering because half the ROI's M W F decreased and half increased. The anisotropic filtering of the M E data was shown to improve the M W F ROI measurement compared to the standard. I expected that the mean M W F 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 M W F calculated from the 4-average, unfiltered spin-echo data. The results did not corroborate this. To understand this, I analyzed M W F histograms in each ROI (for example, see Figure 2.7). The M W F in the 4-average, unfiltered ROI data had a distinct peak (top right plot in Figure 2.7) corresponding 50 C H A P T E R 2. M W F NOISE R E D U C T I O N 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 C H A P T E R 2. M W F NOISE R E D U C T I O N 2.6. SPATIAL FILTERING VS AVERAGING to a true M W F that did not vary spatially in the ROI. The M W F 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 N N L S 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 M R I 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 i n v e r t - a v e r a g e 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 M W F (16% for 5 iterations for the 4-average data) between the i n v e r t - a v e r a g e and b o o t s t r a p methods, the relative difference decreased from over 29% for the unfiltered i n v e r t - a v e r a g e to 16%. The M W F calculated from the i n v e r t - a v e r a g e method tended to be closer to the standard as more iterations of the anisotropic filter were applied to the original M E data. Therefore, the myelin water maps calculated from filtered M E data do tend to be representative of the local M W F . The myelin water fractions were greater in 72% of the 40 ROIs as calculated by the i n v e r t - a v e r a g e method compared to the b o o t s t r a p method (on the 4-average data). This apparent increase in the myelin water fraction was due to the non-linear N N L S algorithm applied to lower SNR data in the i n v e r t - a v e r a g e method. The anisotropic diffusion filter applied to the 4-average, M E 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 M R I 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 N N L S algorithm 52 C H A P T E R 2. M W F NOISE R E D U C T I O N 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 M W F 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 M W F estimated from multi-echo data of different number of averages and with the application of no filter, a n i s o l , aniso3, and aniso5 filter were then compared. The expectation 53 C H A P T E R 2. M W F NOISE R E D U C T I O N 2.7. CONCLUSIONS was the M W F 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 N N L S algorithm. The noise reduction filter was good for reducing the variability of the M W F over assumed homogeneous regions. 54 Chapter 3 M y e l i n W a t e r F r a c t i o n f r o m L i n e a r C o m b i n a t i o n 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 T 2 ~ 20 ms and medium compartment with T 2 ~ 80 ms [18]. The ratio of the areas of the short to medium compartments is approximately 1:6. The short T 2 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 T 2 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 N N L S algorithm minimized the total signal of the T 2 distribution, such that l.02xmin < X2 < ^-•^^Xmin where Xmin w a s the minimum misfit of the unregularized solution. A n alternative method to calculate the short T 2 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 ( T 2 ) is the filter of interest, c is a vector of coefficients. The image corresponding to the (3.1) 55 C H A P T E R 3. L I N E A R C O M B I N A T I O N 3.1. INTRODUCTION filter is defined as: N 7L(x) = ^c,/ i(x) (3.2) i=i where 7i(x) = 2~ ĵ=i P j ( x ) e x P ( — f /T 2 ;j(x)). The goal, then, is to define coefficients, c, such that A ( T 2 ) includes T 2 between 10 and 50 ms and excludes all other T 2 . The method to calculate the coefficients c is the novel work in this chapter. 3.1.1 Linear Combination for T 2 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 T 2 distribution. This method is described further in Section 3.6.1. 3.1.1.2 Calculation of Coefficients by Backus-Gilbert Method In 1991 Whittal l [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 ( T 2 ) , and therefore allows a better selection. 56 C H A P T E R 3. L I N E A R C O M B I N A T I O N 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 ) A 5 ( 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 ) ^ f f l ( 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 T 2 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 T 2 < 50 ms relative to the total water, I calculated two sets of coefficients, one set that selects short T 2 , c s , and a second set that selects all T 2 , c a . To verify the linear combination was selecting the appropriate T 2 components, a set of nickel/agarose phan- toms of known T i / T 2 times were scanned using a multi-echo sequence and images of the short 57 C H A P T E R 3. L I N E A R C O M B I N A T I O N 3.2. C A L C U L A T I O N O F C O E F F I C I E N T S 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 N N L S 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 N N L S 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 T 2 of interest and zero elsewhere. Coefficients for both filters, J 4 S ( T 2 ) and A a ( T 2 ) , 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 1 0 < T2 < 5 0 ms and the second set, ca, that selects T2 > 10 ms. The selection filter for the short T 2 coefficients was defined as As (T2) = ^2iLi ct e x P ( ~ ^ i / T 2 ) (and similarly for the total T 2 ) . The desired filter T ! s ( T 2 ) is a boxcar: 1 T lower < To < T o u p p e r AS(T2) = { _ 2 _ 2 ( 3 7 ) 0 elsewhere and for j4 a (T 2 ): , 1 T 2 > 10 ms Aa (T 2 ) = { ~ (3.8) 0 T 2 < 10 ms though neither definition is able to be constructed from a finite set of exponential decays. Similarly, the desired filter for ^4 a (T2) is 1 for T 2 > 10 ms (a step function). Therefore, the goal was to create a filter, J 4 5 ( T 2 ) , 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 ^S j - e x p [-<i/T 2 j-(x)] + a, (x) (3.9) i=i where p (x) is the proton density, k represents R F inhomogeneity and other factors, s, is the frac- 58 CHAPTER 3. LINEAR COMBINATION 3.2. CALCULATION OF COEFFICIENTS tional contribution of component T 2 j 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 : S N R / L = ^ W ^ i ^ x P ( - ^ Z l M ) . ( 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 T 2 Coefficients, cs, for short T 2 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: S N R < = 8 N B . E E • • • * « P ( - V T . * ) ( 3 . 1 2 ) T 2 0 was defined to be 25 ms (as opposed to 20 ms) to force a broader filter in the short T 2 range of interest. Four constraints were necessary in the minimization procedure for calculating c s : 1. AS{T2) > 0 for 10 < T 2 < 70 ms Force the filter A S (T 2 ) to be positive in the range of interest for the short T 2 . Previous to adding this constraint, the algorithm tended to find a minimum where A S (T 2 ) fluctuated positive and negative in the range 10 < T 2 < 70 ms. 2. A S (T 2 ) < 5 for T 2 > 80 ms where 6 = max [A5(T2)] /100 This is the key constraint to suppress signal from T 2 outside of the range of interest. Intu- 59 C H A P T E R 3. L I N E A R C O M B I N A T I O N 3.2. CALCULATION OF COEFFICIENTS itively, the magnitude of the filter response for T 2 > 80 ms, relative to the filter response for T 2 = 15 ms, must be much less than a factor of six, as the area of the medium T 2 component in white matter is approximately six times larger than the area of the short T 2 component. The factor of 100 was chosen to provide greater suppression of the medium T 2 component. 3. £ c s = 0 The filter A S ( T 2 ) must go to 0 at long T 2 . 4. fAs.dT2 = 1 ± 0 . 0 1 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 ' " C M t 0.02 < 0.01 0 - 0 . 0 1 0 50 100 150 200 T 2 (ms) Figure 3.1: Constraints on the short T 2 filter. The physical interpretation of the constraints are shown in Figure 3.1. 3.2.2 Coefficients for A l l T 2 Similarly, coefficients to select for T 2 > 10 ms, c a , were calculated by minimizing the inverse of the SNR in Equation 3.11 based on a set of assumptions and constraints. 5 A s(T2)>0 t 60 C H A P T E R 3. L I N E A R C O M B I N A T I O N 3.2. CALCULATION OF COEFFICIENTS The function to maximize was: S N R o E l = i < e x p ( - t l / T 2 i 0 ) S N R a = (3.13) where T 2 0 = 25 ms. Two constraints were included in the minimization procedure for calculating c a : 1. A a ( T 2 ) > 0 for T 2 > 10 ms Force the filter A a ( T 2 ) to be positive for T 2 greater than 10 ms. 2. | A a ( T 2 ) - M\ < S for T 2 > 10 ms where M = max L4 a(T 2)], 6 = M / 4 0 This wil l constrain the fluctuations in the filter response to be less than 1/40 for T 2 > 10 ms. 0.01 0.009 0.008 0.007 ~ 0.006 a < 0.005 0.004 0.003 0.002 0.001 0 A a (T 2 )>0 50 100 T 2 (ms) 150 200 Figure 3.2: Constraints on the "all" T 2 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, M A ) . The short T 2 fraction was the signal calculated as 61 C H A P T E R 3. L I N E A R C O M B I N A T I O N 3.2. CALCULATION OF COEFFICIENTS the linear combination of the short T 2 coefficients with the images divided by the signal calculated from linear combination of the total T 2 coefficients with the images. The constraints, incorporated into the minimization algorithm, required a T 2 filter be created, at each iteration, based on the current set of coefficients. The requirements placed on the T 2 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 T 2 distribution. Small S were tried but the minimization algorithm was not able to find a feasible solution. 3.2.4 Results The coefficients, c s and c a calculated based on the minimization procedure above, are shown i i i Table 3.1. The corresponding filters, As (T 2 ) and Aa (T 2 ) , are shown in Figure 3.3. t (ms) c s 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 T 2 filter and c a were for selecting all T 2 . 62 C H A P T E R 3. L I N E A R C O M B I N A T I O N 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: T 2 filters for short T 2 components ( A S ( T 2 ) , solid) and all T 2 components ( A a ( T 2 ) , 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, T 2 sa 15 ms, relative to the noise. The variable T 2 0 was defined slightly longer than that of myelin water to increase the width of the filter ^4 S(T 2). Larger values of T 2 0 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 T 2 (ms) 63 C H A P T E R 3. L I N E A R C O M B I N A T I O N 3.3. LINEAR COMBINATION SIMULATIONS to N N L S 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 M W F calculated by linear combination were similar to the mean and standard deviation of the M W F calculated by N N L S . 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 N N L S algorithms. Other parameters included N = 32 echoes, t = 10, 20 , . . . , 320 ms and T R 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 N N L S solution 64 C H A P T E R 3. L I N E A R C O M B I N A T I O N 3.3. LINEAR COMBINATION SIMULATIONS with the true fs — 0.20. The mean short T 2 fraction was within 0.02 of the true fraction. The short T 2 fractions calculated by the linear combination method were normally distributed for each of the true fractions, but only when fs > 0.10 were the N N L S results normally distributed. On average the linear combination method calculated approximately 10,000 solutions per second and the regularized N N L S 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), N N L S (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 N N L S method. The N N L S 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 C H A P T E R 3. L I N E A R C O M B I N A T I O N 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 N N L S 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 T 2 coefficients, in Table 3.1, was analyzed. The short T 2 image was calculated as the linear combination of the multi-echo data to show the selectivity of the short T 2 phantoms and the suppression of the longer T 2 phantoms. 3.4.1 Methods The phantoms (Section 1.8) were scanned using the 32 echo M E S E 0 pulse sequence (Section 1.6.2). The images were combined linearly using the short T 2 filter weights defined in Table 3.1. 3.4.2 Results and Discussion Excellent suppression of the longer T 2 phantoms was found as shown in Figure 3.5. The T 2 of the nine phantoms surrounding the water bottle were (clock-wise starting from bottom left) T 2 = 264,101,25,23,97,363,157,77,23 ms. The short T 2 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 T 2 phantoms in the short T 2 map was at least fifteen times higher than the mean signal of the longer T 2 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 T 2 of approximately 100 ms compared to the phantoms with T 2 greater than 200 ms. This was not unexpected as the T 2 filter is close to zero near 100 ms, but not zero. Overall the quantitative selectivity of the short T 2 phantoms was excellent. 66 C H A P T E R 3. L I N E A R C O M B I N A T I O N 3.4. MONO-EXPONENTIAL PHANTOMS Figure 3.5: M R I image (TE =10 ms, left) of nine nickel/agarose phantoms (and the large water bottle) and the short T 2 image (right) based on the linear combination of the original 32 spin-echo images. Region Mean Signal T 2 (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 T 2 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 C H A P T E R 3. L I N E A R C O M B I N A T I O N 3.5. NNLS VS. T 2 F I L T E R 1.2r T 2 (ms) Figure 3.6: The mean signal from each phantom plotted against the T 2 of the phantom overlayed on the short T 2 filter. 3.4.3 Conclusions The selectivity of the short T 2 phantoms, from the linear combination method was excellent as was the suppression of the long T 2 phantoms. For mono-exponential phantoms, the filter behaved as well as expected. 3.5 In Vivo Myelin Water Fraction: NNLS vs. T 2 Filter The selectivity and suppression of mono-exponential phantoms was excellent. In this section, the linear combination method is applied to a mult i-T 2 component system. The short T 2 fraction calculated by the T 2 filter method was compared to those calculated by N N L S (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 C H A P T E R 3. L I N E A R C O M B I N A T I O N 3.5. NNLS VS. T 2 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 N N L S method. Six regions (each left and right) were selected (on the T E = 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 T 2 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 N N L S . 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 T 2 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 T 2 in the short compartment. A n increase in the short T 2 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 N N L S . Gray matter had 90% ofthe differences within 1.96CT and a slight bias toward larger fractions calculated by N N L S . 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 C H A P T E R 3. L I N E A R C O M B I N A T I O N 3 . 5 . NNLS VS. T 2 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 T 2 nickel/agarose phantoms as shown in Figure 3.5 in 17.5% of the regions for the linear combination method. 70 C H A P T E R 3. L I N E A R C O M B I N A T I O N 3.5. NNLS VS. T 2 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 N N L S (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 . 0 2 c 0) CD - 0 . 0 4 b - 0 . 0 6 - 0 . 0 8 -oj ! 1.96'SD • • • ; • • • • Mean > * i» ! ! . - 1 96 S O * 0.06 0 0 8 0.1 0.12 0.14 0 16 0.18 0 i Mean of Pairs 0.06, 0.05 0.04| 0.03 0.02 0.01 0 - 0 . 0 1 - 0 . 0 2 -0.0.3 ! > ' 1 i • 1.96 S D • • • - • i • -1 .96 S D -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 C H A P T E R 3. L I N E A R C O M B I N A T I O N 3.6. COEFFICIENTS B Y MATRIX INVERSION 3.5.3 Conclusions Myelin water fractions calculated from linear combinations of T 2 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 T 2 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 wil l be very broad compared to the filters defined in Section 3.2. In this section, I calculated coefficients to select for short T 2 (cs, 10 < T 2 < 50 ms), medium T 2 (cm, 70 < T 2 < 120 ms), and long T 2 (c(, T 2 > 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 (—£,/T 2 j ) ] , and L was the amplitudes which correspond to the T 2 . The coefficients were the rows of F~l where Fij — [exp (—i j /T 2 j ) ] for t = 10, 2 0 , . . . , 320 ms and T 2 = 20, 80,3000 ms. The coefficients calculated are shown in Table 3.3. The coefficients for the short T 2 , cs, follow a similar trend of positive in the early echoes, then negative, and then positive. The T 2 filters, / ( T 2 ) , shown in Figure 3.3, were constructed by calculating the filter / ' ( T 2 ) = (F~1)t exp ( - i / T 2 ) . The short T 2 filter, top left filter of Figure 3.10, should be compared with the short T 2 filter in Figure 3.3 (solid line). 72 C H A P T E R 3. L I N E A R C O M B I N A T I O N 3.6. COEFFICIENTS B Y MATRIX INVERSION 1000 T2 (ms; T 2 (ms) Figure 3.10: The filters calculated from the matrix inversion technique to calculate coefficients. Short T 2 filter (top left), medium T 2 filter (top right) and the long T 2 filter (bottom left). The short T 2 filter calculated by the inversion method {AS'M(T2), dashed line) compared to the short T 2 filter (^4S(T 2), solid line) calculated by the non-linear method (bottom right, plotted on a log scale to visualize the whole filter). 73 C H A P T E R 3. L I N E A R C O M B I N A T I O N 3.6. COEFFICIENTS B Y 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 T 2 selection based on the matrix inversion method of calculation. 3.6.2 Results and Discussion The short T 2 filter (Figure 3.10, top left) is a very broad peak centered on T 2 = 22 ms and crosses 0 at T 2 = 80 ms. The primary problem with the short T 2 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 T 2 > 80 ms. If the system has only mono-exponential contributions for all voxels, then this is not a problem, but if multiple T 2 compartments are possible per voxel, then this could lead to a mis-interpretation of water compartments 10 < T 2 < 80 ms, as the myelin water will appear lower even though the short T 2 was not lower. The short T 2 filter, calculated from this method, is compared to the short T 2 filter, A S ( T 2 ) 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 T 2 filter is very effective if objects in the image are all mono-exponential. In this case, objects with a short T 2 near 20 ms will be selected and objects with a T 2 > 80 ms will be negative. The medium T 2 filter (Figure 3.10, top right), and long T 2 filter (Figure 3.10, bottom), both have a very broad selection. 74 C H A P T E R 3. L I N E A R C O M B I N A T I O N 3.6. COEFFICIENTS B Y MATRIX INVERSION T 2 Selected Images The coefficients calculated in Section 3.6.1, were applied to a multi-echo M R I dataset. The short T2 (top left), medium T 2 (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 T 2 (top left), medium T 2 (top right), long T 2 (bottom left) and the myelin water fraction (ratio of the short T 2 to the sum of the three, bottom right). 75 C H A P T E R 3. L I N E A R C O M B I N A T I O N 3 . 7 . BACKUS-GILBERT COEFFICIENTS 3.6.3 Conclus ions 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 T 2 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 Ca lcu la t ion of Coeff icients 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 ( T 2 , T 2 ) 0 ) - n ( T 2 , T 2 , m m , T 2 ! m a x ) ] 2 d T 2 + S m(0) V c 2 c r 2 . (3.15) where 6 is the trade-off parameter, A ( T 2 , T 2 > 0 ) = ^ e x p ( - T E / T 2 ] 0 ) , and n(T 2 , T 2 ) m i n , T 2 , m a x ) is the boxcar function defined as: T T f T T \—) 1 ^ ̂ 2<min - T 2 - T 2 , m a x (o T C \ J-2,mm> i-2,max) ~ \ • ( 3 . 1 0 J I 0 otherwise The coefficients which minimize Equation 3.15 are found by calculating the partial derivative of (j> with respect to each coefficient C j , and setting the derivative to zero: _ 0 = O i = l...N (3.17) 76 C H A P T E R 3. L I N E A R C O M B I N A T I O N 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 T 2 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. A l l four filters had a similar selectivity in the T 2 of interest. The differences in the filter were around T 2 = 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 < T 2 < 1000 ms. The signal from T 2 > 100 ms is undesirable as the filter was designed for short T 2 . 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 C H A P T E R 3. L I N E A R C O M B I N A T I O N 3.7. BACKUS-GILBERT COEFFICIENTS 1 1 I 1 i 1 c b g 1 cbg2 - - c b 9 3 ' . . . . c s • , , , , l t «* \ 1 \ > ^ y - 50 100 150 200 T E (ms) 250 300 350 Figure 3.12: Plot of the coefficients calculated from the Backus-Gilbert technique. Coefficients, c s , calculated in Section 3.2 are shown for reference. 78 C H A P T E R 3. L I N E A R C O M B I N A T I O N 3.7. BACKUS-GILBERT COEFFICIENTS Figure 3.13: The filters calculated from the Backus-Gilbert technique. Filter, A S ( T 2 ) , calculated based on the coefficients cs, is shown for reference. 79 C H A P T E R 3. L I N E A R C O M B I N A T I O N 3.8. O P T I M A L E C H O T I M E S 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, A S ( T 2 ) , resulted in a smaller contribution for T 2 > 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 T 2 . Two sets of experiments were done to determine the quality of the short T 2 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 T 2 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 < T 2 < 50 ms and the T 2 at which the filter reaches 1/100 of the filter maximum. The filter is defined as A ( T 2 ) = E i l i °i e x P ( _ * i / T 2 ) , where N is the number of echoes. The goal of the filter is to be as wide as possible for 10 < T 2 < 50 ms and near zero for T 2 > 80 ms. Figure 3.14 shows a short T 2 filter. Therefore, the measure of the goodness of a filter is based on the full width at half the maximum ( F W H M ) , which should be as wide as possible and the T 2 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 T 2 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 C H A P T E R 3. L I N E A R C O M B I N A T I O N 3 . 8 . OPTIMAL ECHO TIMES 0.045 r -0.005" 1 1 1 1 1 : ' 0 20 40 60 80 100 120 T 2 (ms) Figure 3.14: The measures of the effectiveness of a short T 2 filter. The full width at half the maximum should be as large as possible and the T 2 where the filter reaches 1/100 of the maximum should be less than T 2 = 80 ms and greater than T 2 = 50 ms. 81 C H A P T E R 3. L I N E A R C O M B I N A T I O N 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 T 2 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, A i , 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 M E S E 0 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 T 2 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 T E « 80 ms, and the following coefficients oscillated positive and negative. The top plot of Figure 3.16 shows the width of the F W H M of A S ( T 2 ) . 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 C H A P T E R 3. L I N E A R C O M B I N A T I O N 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 T E . The line segment is the coefficients for two echoes, and appears out of place only because of the perspective of the axes. 83 C H A P T E R 3. L I N E A R C O M B I N A T I O N 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 F W H M the better. 84 C H A P T E R 3. L I N E A R C O M B I N A T I O N 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 T E , 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 T 2 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 T R and T i weighting. max t,c Z i = i ^ e x p ( - ^ / 2 5 ) (3.19) defined to constrain the noise amplification 85 C H A P T E R 3. L I N E A R C O M B I N A T I O N 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 F W H M 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 T E . 86 C H A P T E R 3. L I N E A R C O M B I N A T I O N 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 100 t h max (bottom). 87 C H A P T E R 3. L I N E A R C O M B I N A T I O N 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 T 2 space. 3.8.3.3 Conclusions The short T 2 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 T 2 filter, shorter echo times would be required, which would be difficult on an M R I machine, but would be possible on an N M R spectrometer. 3.8.4 Discussion The selectivity of the short T 2 filter, A 5 ( T 2 ) , for 10 < T 2 < 50 ms is a function of the echo and corresponding coefficients. Figure 3.19 shows the T 2 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 T 2 < t. Therefore, a sharp rise on the left side of A S ( T 2 ) 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 T 2 range of interest. The signal from each echo contributes different information depending on the echo time of the acquisition. The 32 coefficients, c s , 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 A S ( T 2 ) 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 C H A P T E R 3. L I N E A R C O M B I N A T I O N 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 T 2 , based on echo time T E = 10 ms (upper solid line), T E = 80 ms (dash-dotted line), T E = 120 ms (dashed line), T E = 200 ms (dotted line), and T E = 800 ms (lower solid line). Each filter is zero at T E = 0 ms and goes to one as T 2 —> 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 T 2 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 C H A P T E R 3. L I N E A R C O M B I N A T I O N 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 M T 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 T E (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 T 2 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 T 2 filter. The first echo time, T E = 10 ms, defines the rapid increase in the portion of the curve for 0 < T 2 < 18 ms. The second and third echoes, at T E = 30 ms and T E = 120 ms, bracket a portion of the T 2 filter where the filter changes from decreasing to flat. The fourth echo, at T E = 270 ms forces the tail of the filter to be near zero at long T 2 . The short T 2 filter (Figure 3.21) calculated from the four echoes followed the T 2 filter cal- culated from 32 echoes. The ascending portion of the filter (0 ;$ T 2 ;$ 15 ms) was exactly the same in each as that portion of the curve is determined by the T E = 10 ms contribution only. The descending portion of the filter (15 ~ T 2 ;$ 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 T 2 filter, from 4 echoes, was not as close to zero as the short T 2 filter from 32 echoes for T 2 > 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 T 2 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 C H A P T E R 3. L I N E A R C O M B I N A T I O N 3.8. OPTIMAL ECHO TIMES Figure 3.21: Short T 2 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 T 2 . several echoes at small T E , one echo near the T 2 where the filter dropped to zero and one echo at long T E . 91 C H A P T E R 3. L I N E A R C O M B I N A T I O N 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 T 2 ~ 20 ms. 3.9 C o n c l u s i o n s 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 T 2 voxel-by-voxel. The myelin water fraction is calculated as the ratio of the signal of water with 10 < T 2 < 50 ms to the total signal in the distribution. The calculation of the T 2 distribution for each voxel in an image is a slow process. Another method to estimate the signal of water within a range of T 2 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 < T 2 < 50 ms (short T 2 ) and a second set to select water with T 2 > 10 ms (all T 2 ) . The coefficients to select for short T 2 were applied to multi-echo images of a set of phantoms. The resulting short T 2 image had very good selection of the phantoms with T 2 near 20 ms and good suppression of phantoms with T 2 > 70 ms. The N N L S 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 C H A P T E R 3. L I N E A R C O M B I N A T I O N 3.9. CONCLUSIONS the in vivo myelin water fraction, from five volunteer M R I 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 T 2 . It was 20,000 times faster than the regularized N N L S 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 N o n - 1 8 0 R e f o c u s i n g P u l s e s 4.1 Introduction 4.1.1 Motivation Accurate T 2 estimation traditionally depends on accurate refocusing pulses. Many studies [31, 35, 71, 72] showed large deviations in the measured T 2 when the refocusing pulse flip angle deviates from 180°. For example, Majumdar [71] measured biases in R 2 (R 2 = 1/T 2) 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 T 2 decay curves. Levitt [30] proposed using composite refocusing pulses (QO^-lSO^-QOx) to correctly rephase spins in the presence of moderate B i inhomogeneity. Zur and Stokar [73] analyzed artifacts due to inhomogeneous B i 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 R F pulses and of composite R F pulses. In 2000, Sled and Pike [75] described a method to correct inaccurate mono-exponential T 2 measurements due to imperfect refocusing pulses. A multi-echo pulse sequence was used to calculate the T 2 and a separate acquisition was required to calculate the B i map. The calculated T 2 was divided by f (a) to obtain.a "corrected" mono-exponential T 2 - Another approach by Lepage [76] corrected 94 C H A P T E R 4. NON-180 R E F O C U S I N G P U L S E S 4.1. INTRODUCTION the experimentally determined T 2 of a gel phantom based on a B j map calculated from a dif- ferent homogeneous gel phantom. The techniques by Sled and Lepage each required a separate acquisition to calculate the B i 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 e x p ( - 2 r / T 2 ) for F(t) and F*(t), e x p ( - 2 r / T i ) for Mz(t), and a generating term MQ (1 — exp(—2r/Ti)) (where 2r is the echo spacing). For accurate T 2 quantification, Hennig 95 C H A P T E R 4. NON-180 R E F O C U S I N G P U L S E S 4.1. INTRODUCTION suggested all magnetization from stimulated and indirect echoes must be removed by spoilers. On the contrary, accurate and precise T 2 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 C H A P T E R 4. NON-180 R E F O C U S I N G P U L S E S 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) = _ ^ L Z p ) + W l V W (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 R F pulse): U{t + tw) = U(t) V(t + tw) = V(t) cos (0) - Mz sin (0) Mz(t + tw) = V(<)sin(0) + M 2cos(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 C H A P T E R 4! NON-180 R E F O C U S I N G P U L S E S 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) + sin 2 (9/2)] + iV(t) [cos2 (9/2) - sin 2 (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) s in 2 (9/2) - iMz (t) sin (9) (4.7) = U(t) cos2 (0/2) + iV(t) cos2 (0/2) + r/(i) s in 2 (9/2) - iV(t) sin2 (9/2)-iMz(t) sin (9) • (4.8) = [U(t) + iV(t)} cos 2 (9/2) + [C/(t) - iV(t)] s in 2 (0/2) - i M 2 (t) sin (0) (4.9) = F( i ) cos 2 (0 /2 ) + F * ( t ) s i n 2 ( 0 / 2 ) - i M z ( i ) s i n ( 0 ) (4.10) by using the cosine double angle formula: cos (0) = cos2 (0/2) - sin 2 (0/2). (4.11) and the identity: cos2 (0/2) + sin 2 (0/2) = 1. (4.12) Then, F(t + tw) and Mz(t + tw), written as functions of F(t) and Mz(t), wil l be: F(t + tw) = F ( i ) cos 2 (0 /2 ) + F * ( i ) s i n 2 ( 0 / 2 ) - i M 2 ( i ) s i n ( 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 Implementat ion of A l g o r i t h m The algorithm toxcalculate the signal of each echo from a train of R F 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 C H A P T E R 4. NON-180 R E F O C U S I N G P U L S E S 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 T i (Algorithm 4.2) and T 2 (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 T 2 decay, but there is also a weak signal dependence on the T i relaxation time. For the primary work in this chapter, involving the brain, the T i « 650 ms [80]. The T i 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 / T i is not zero. The decay curve dependence on T i could possibly be exploited as a means of calculating the T i along with T 2 , a, and p. The signal intensity, y, was calculated (Algorithm 4.1) given T 2 , T i , 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 R A M . 99 C H A P T E R 4. NON-180 R E F O C U S I N G P U L S E S 4.1. INTRODUCTION A l g o r i t h m 4.1 Signal calculation from a train of refocusing pulses f u n c t i o n [decay] = s t i m u l a t e ( a l p h a , num_echoes, T l , T 2 , t a u ) % 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) j a y = s q r t ( - l ) ; % Preallocate the signal intensity vector d e c a y = z e r o s G , num_echoes) ; % Preallocate the array of num^echoes substates [F F* Z Z*] and initialize B = zeros (num_echoes*4 , 1 ) ; B ( l ) = 0 + j a y * l ; % Precalculate rotation matrix Tone for one substate block Tone = [ [ ( c o s ( a l p h a / 2 ) ) " 2 , ( s i n ( a l p h a / 2 ) ) " 2 , - j a y * s i n ( a l p h a ) , 0 ] ; . . . [ ( s i n ( a l p h a / 2 ) ) ~ 2 , ( c o s ( a l p h a / 2 ) ) " 2 , 0 , j a y * s i n ( a l p h a ) ] ; . . . [ ( - j a y ) * l / 2 * s i n ( a l p h a ) , ( j a y ) * l / 2 * s i n ( a l p h a ) , c o s ( a l p h a ) , 0] ; . . . [ ( - j a y ) * l / 2 * s i n ( a l p h a ) , ( j a y ) * l / 2 * s i n ( a l p h a ) , 0 , c o s ( a l p h a ) ] ] ; % Calculate the sparse rotation matrix of "num-echoes" blocks Tp = k r o n ( s p a r s e ( d i a g ( o n e s ( l . n u n u e c h o e s ) ) ) , c r e a t e . T p ( a l p h a ) ) ; % Apply initial T2 decay and Tl regrowth (between 90-180 pulses). B = t 2 r e l a x ( B , t a u , T 2 ) ; B = t l r e l a x ( B , t a u , T l ) ; f o r 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* d e c a y ( e c h o ) = a b s ( B ( 2 ) ) * e x p ( - t a u / T 2 ) ; % Apply transverse decay for time 2*tau. B = t 2 r e l a x ( B , 2 * t a u , T 2 ) ; % Apply longitudinal regrowth for time 2*tau. B = t l r e l a x ( B , 2 * t a u , T l ) ; % Let the spins dephase for time 2*tau B = e v o l u t i o n ( B ) ; end 100 C H A P T E R 4. NON-180 R E F O C U S I N G P U L S E S 4.1. INTRODUCTION Algorithm 4.2 T i r e l a x a t i o n f u n c t i o n [B] = t l r e l a x ( B , t a u , T l ) % The Z and Z* substates will be modified i n d i c i e s = [ 3 : 4 : l e n g t h ( B ) 4 : 4 : l e n g t h ( B ) ] ; % Apply the Tl relaxation to the Z and Z* substates B ( i n d i c i e s ) = B ( i n d i c i e s ) * e x p ( - t a u / T l ) ; Algorithm 4.3 T 2 r e l a x a t i o n f u n c t i o n [B] = t 2 r e l a x ( B , t a u , T2) % The F and F* substates will be modified i n d i c i e s = [ 1 : 4 : l e n g t h ( B ) 2 : 4 : l e n g t h ( B ) ] ; % Apply the T2 relaxation to the Z and Z* substates B ( i n d i c i e s ) = B ( i n d i c i e s ) * e x p ( - t a u / T 2 ) ; Algorithm 4.4 Phase evolution of magnetization f u n c t i o n [Bnew] = e v o l u t i o n ( 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 = l e n g t h ( B ) ; % All F(n) go to F(n+1). Bnew([ 5 : 4 : L ] ) = B ( [ l : 4 : L - 4 ] ) ; % F*(l) goes to F(l). B n e w ( l ) = B ( 2 ) ; % All F*(n) go to F*(n-1). Bnew([ 2 : 4 : L - 4 ] ) = B ( [ 6 : 4 : L ] ) ; 101 C H A P T E R 4. NON-180 R E F O C U S I N G P U L S E S 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 M E S E C 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 M E S E C sequence were compared to myelin water maps calculated from the M E S E 0 . 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 T 2 Decay Curve Fit A n 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 T 2 values, were scanned with different sets of refocusing pulses. The estimated T 2 was compared to the T 2 estimated from the M E S E 0 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 T 2 , baseline offset d, and refocusing pulse flip angle a. A non- linear least squares algorithm ( lsqcurvef i t in Matlab [The Mathworks, Natick, MA]) minimized the misfit x 2 = 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, T 2 = 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]. T R was considered infinite for all fitting. 102 C H A P T E R 4. NON-180 R E F O C U S I N G P U L S E S 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. M e t h o d s A noiseless decay curve was created based on the model: y = p e x p ( - t / T 2 ) (4.17) where p = 1000 and T 2 = 20,80,140 ms. The decay curve was created based on flip angles of a = 90°, 1 0 0 ° , . . . , 180° and echo times t = 10, 20, . . . , 320 ms. One thousand realizations of quadrature noise (SNR =200) were added. T R was set to infinity. The T 2 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. Resu l t s and Discuss ion Each parameter estimated from the decay curve, for T 2 = 20,80,140 ms was within 5% of the true parameter. Table 4.1 shows the results for T 2 = 80 ms. The parameters estimated from the decay curve were within 2% of the true parameter for T 2 = 80,140 ms. The increased accuracy, as a function of T 2 , was expected as more echoes were able to contribute information to the solution as T 2 increased. The variability in the estimated parameters was higher for the lower flip angles (for example, variability in T 2 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, sin 2 (a/2), on the relative attenuation due to the flip angle. For example, I multiplied the random noise by sin 2 (a/2) = 1.0 for a = 180° and by sin 2 (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 C H A P T E R 4. NON-180 R E F O C U S I N G P U L S E S 4.2. MONO-EXPONENTIAL FIT "true ( ° ) A(° ) P T 2 ( 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 Oi l ) 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 T 2 (standard deviation in brackets) given the true p = 1000 and true T 2 = 80 ms. "true ( ° ) a{°) P T 2 ( 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 T 2 from a simulation with SNR scaled as a factor of sin 2 (a/2) to determine whether the increased variability in the estimated parameters was due to the decrease in a. 104 C H A P T E R 4. NON-180 R E F O C U S I N G P U L S E S 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 T 2 remain relatively constant as a decreased to 90°. The constant variability was true for T 2 = 140 ms as well. When T 2 = 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 T 2 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 T 2 was affected more by the decreased flip angle than the longer T 2 components. 4.2.3 Phantom Verification A set of mono-exponential phantoms were used to confirm the accuracy of T 2 from a set of scans. A set of six nickel/agarose phantoms of known T 2 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 M E S E C pulse sequence with parameters: echo spacing = 11.1 ms, T R = 3000 ms and slice selective refocusing pulses. Regions were drawn on the T E = 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 (T 2 = 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 C H A P T E R 4. NON-180 R E F O C U S I N G P U L S E S 4.2. M O N O - E X P O N E N T I A L F I T alternated positive and negative for successive echoes. The curve fitting algorithm worked well over a range of refocusing pulse flip angles and T 2 (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 T 2 was typically within ±10% of the T 2 estimated from the 48 echo M E S E 0 data. OLp "best P180 PP Pbest T 2 ,180 t 2 , p t2,best X l 8 0 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 (T 2 = 90.8 ms). Fl ip angles are in degrees and T 2 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 R F 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 T 2 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 T 2 , 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 C H A P T E R 4. NON-180 R E F O C U S I N G P U L S E S 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. T 2 = 90.8 ms for the phantom. 107 C H A P T E R 4. NON-180 R E F O C U S I N G P U L S E S 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. T 2 = 90.8 ms for the phantom. 108 C H A P T E R 4. NON-180 R E F O C U S I N G P U L S E S 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. T 2 = 90.8 ms for the phantom. 109 C H A P T E R 4. NON-180 R E F O C U S I N G P U L S E S 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. T 2 = 90.8 ms for the phantom. 110 C H A P T E R 4. NON-180 R E F O C U S I N G P U L S E S 4.3. MULTI-EXPONENTIAL FIT Conclus ions 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 T 2 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 M u l t i - E x p o n e n t i a l F i t A l g o r i t h m 4.3.1.1 T 2 Decay C u r v e F i t U s i n g N N L S Multi-exponential fits of the in vivo data were done with the non-negative least squares (NNLS) algorithm [18, 39, 40]. N N L S solves 4̂s = y where Aij = exp ( - T E i / T 2 j ) , y is the measured decay curve and Sj are the unknown amplitudes that correspond to T 2 j . The solution minimizes the x2 misfit. In the standard N N L S algorithm, each column Aj, of the matrix A, is a mono- exponential decay curve with T 2 = T 2 j . 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 T 2 = T 2 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, N N L S a , to find the minimum x2 (a) by solving Aas = y for a between 40° and 180°. The multi-exponential fit algorithm I introduced, N N L S " , minimized x2 a s 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 C H A P T E R 4. NON-180 R E F O C U S I N G P U L S E S 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 / 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 N N L S " algorithm is to find the minimum x2 over a. 112 C H A P T E R 4. NON-180 R E F O C U S I N G P U L S E S 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 T 2 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, T E = 10 ms, T i = 800 ms. The second set was based on 200 echoes, T E = 5 ms, T i = 800 ms. The T i was chosen to be 800 ms as that is approximately the T i 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 C H A P T E R 4. NON-180 R E F O C U S I N G P U L S E S 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 T E = 10 ms and T R = 800 ms. ratio (smallest to largest singular values) was 105.2686 x 10~ 1 7 at a = 158° and the lowest was 0.4996 x 10~ 1 7 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 1 0 - 1 7 at a = 150° and the lowest was 0.0021 x 10~ 1 7 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 A 2 0 8 , but the difference between The mono-exponential decay curves in A 2 0 8 and A 1 5 2 was less than I O - 1 3 . Conclusions The matrix A 1 8 0 has the lowest rank and the highest condition number and therefore is about the worst matrix with which to calculate T 2 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 C H A P T E R 4. NON-180 R E F O C U S I N G P U L S E S 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 T E = 5 ms and T R = 800 ms. number (two orders of magnitude) should allow for better T 2 quantification. Another benefit of a multi-echo pulse sequence using 159° non-composite pulses is that the S A R 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 M R I 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 T 2 compartment of p s = 200 and T 2 S = 20 ms and p m = 800 and T 2 m = 80 ms. A true decay curve was created using the bi-exponential model parameters and refocusing pulse flip angles of a t r U e = 90°, 1 0 0 ° , . . . , 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 C H A P T E R 4. NON-180 R E F O C U S I N G P U L S E S 4.3. MULTI-EXPONENTIAL FIT The myelin water fraction, proton densities and T 2 s had less than 10% error for a > 170° (Table 4.4). A l l parameters were incorrectly determined when a < 170°. Therefore, accurate and precise quantification of T 2 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, i y = 20, fm = 0.80, dnm = 800, T 2 m = 80. 4.3.3 N N L S Q M u l t i - E x p o n e n t i a l S imulat ions A set of noise simulations were used to assess the accuracy and precision of the parameters of the T 2 distribution and refocusing pulse flip angle. A bi-exponential T 2 distribution as in Section 4.3.2 was used. The true parameters were fs = 0.20, T 2 S = 20 ms, fm = 0.80, and T 2 m = 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/sin 2 (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 N N L S Q algorithm. A l l 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 C H A P T E R 4. NON-180 R E F O C U S I N G P U L S E S 4.3. M U L T I - E X P O N E N T I A L F I T 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 T 2 ) , fm T 2 ,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°, 1 0 0 ° , . . . , 180°. The true bi-exponential distribution was / s = 0.2, T2 ] S = 20 ms, fm = 0.8 and T2,m = 80 ms. Fl ip angles are in degrees and T2 in milliseconds. 4.3.4 Phantoms The multi-exponential algorithm, N N L S " , to calculate the multi-component T2 was validated based on a set of M R I scans of the nickel/agarose phantoms. The phantoms were scanned with the M E S E C sequence, from which, the T2 were compared to T2 calculated from the M E S E G sequence. The refocusing pulse flip angles were calculated from the M E S E C sequence, as well, and compared to flip angles calculated from a technique by Stollberger [81]. The T i 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 M E S E 0 pulse sequence with T E = 10, 20, . . . , 320 ms. Other parameters were T R = 3000 ms, 2 averages. The T 2 measurements were done with a 32 echo M E S E C pulse sequence with and echo spacing of 12.1 ms. Other parameters were T R = 3000 ms, 2 averages. The B i measurements were calculated from a pair of spin echo images [81], the 117 C H A P T E R 4. NON-180 R E F O C U S I N G P U L S E S 4.3. MULTI-EXPONENTIAL FIT first image was a standard spin-echo image, the second had an excitation pulse of 180°. Other parameters included T E = 11.1 ms, T R = 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 N N L S . The B i , 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 - e x p ( - T R / T i ) ] . 4.3.4.2 Results T 2 Calculation The T2 relaxation times, shown in Table 4.6, were calculated from the M E S E 0 sequence (for reference) and from the M E S E C sequence. The T2 relaxation times calculated from the decay curve collected by the M E S E C sequence were within 8% of the relaxation times calculated from the decay curve collected by the optimized M E S E 0 sequence. The T i 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 MESE Q T 2 (ms) MESE C T 2 (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: T 2 time calculated from the M E S E 0 sequence and the MESE C (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 MESE C 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 MESE C sequence (a). Phantom T i (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 T i from curve fits of five points (TR = 100, 300, 800, 2000, 3000 ms) from a pulse saturation experiment. 119 C H A P T E R 4. NON-180 R E F O C U S I N G P U L S E S 4.3. MULTI-EXPONENTIAL FIT 4.3.4.3 Conclusions The nickel/agarose phantoms were scanned using the M E S E C 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 T 2 from the M E S E 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 a n o m = 180° and for the second scan, the flip angle was set to a n o m = 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 M E S E C with T R — 3 s, T E = 11.096 ms, 256 x 128 freqxphase, 24 cm F O V and 4 averages. The first scan was acquired as per normal with a n o m = 180°. The second scan was acquired with a n o m = 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 M E S E 0 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 Fl ip Angle") in the image and for the second analysis the flip angle was assumed a = 180° at every voxel ("Fixed Fl ip Angle"). 4.3.5.2 Results Variable Flip Angle Maps of the myelin water were estimated from the multi-echo data with a n o m = 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 M E S E 0 . The map from the 120 C H A P T E R 4. NON-180 R E F O C U S I N G P U L S E S 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 a n o m = 150° (Figure 4.9, right). The histogram of the myelin water fraction, within the brain, was narrower for a n o m = 150° than for a n o m = 180°. As well, there was more myelin water structure visible in the middle portion of the myelin water map acquired with a n 0 m = 150°. The bright phantoms on the right side of each myelin water map are the short T 2 nickel/agarose phantoms. Region M W F ( a n o m = 150°) M W F ( a n o m = 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 a n o m = 150° and a n o m = 180°. 121 C H A P T E R 4. NON-180 R E F O C U S I N G P U L S E S 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 a n o m = 150°, and are the top pair of images in Figure 4.10. The mean refocusing pulse flip angle from M E S E C with « n o m = 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 C H A P T E R 4. NON-180 R E F O C U S I N G P U L S E S 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 a n o m = 150° and a n o m = 180°. 120 100 Figure 4.11: The arithmetic mean T 2 maps calculated from the pulse sequence with a refocusing pulse flip angle prescribed as 180° (left) and 150° (right). 123 C H A P T E R 4. NON-180 R E F O C U S I N G P U L S E S 4.3. MULTI-EXPONENTIAL FIT The geometric mean T 2 of the medium T 2 component estimated from M E S E C with ttn0m = 180° and a M m = 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 T 2 m from the two scans were within 4%. Region T 2 m (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 T 2 of the medium component from select regions throughout the image. The x2) calculated during the estimation of the T 2 distribution, is shown in Figure 4.12. The regions of high \ 2 • , a n d 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 T 2 distribution of CSF. 124 C H A P T E R 4. NON-180 R E F O C U S I N G P U L S E S 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 T 2 distribution was calculated voxel-by-voxel from the two scans assuming a n o m = 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 a n 0 m = 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 C H A P T E R 4. NON-180 R E F O C U S I N G 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 N N L S algorithm was applied assuming, incorrectly, the refocusing pulse flip angle was 180°. 4.3.5.3 Conc lus ions The myelin water map estimated from both scans was similar to myelin water maps estimated from an optimized M E S E 0 sequence. Quantitatively the estimated percent myelin water, refocusing pulse flip angle, and geometric mean T 2 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 T 2 decay curve. 4.4 Mul t i -Sl ice , Mul t i -Echo 4.4.1 Introduction The current multi-echo acquisition method, M E S E 0 , is able to collect only a single slice within the given scan time, due to the hard, composite pulses. The M E S E C 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 C H A P T E R 4. NON-180 R E F O C U S I N G P U L S E S 4.4. MULTI-SLICE, MULTI-ECHO use of the myelin water imaging. I scanned a volunteer using the 3D M E S E C sequence. Six slices were acquired, within the brain, as well as a slice acquired with the M E S E 0 for comparison. Myelin water maps were calculated for each slice of the M E S E C sequence and the single slice of the M E S E 0 . 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 M E S E 0 with 4 averages and T E = 10 ms; and second, a 6 slice slab, 42 echo 3 D - M E S E C sequence with 1 average and T E = 10.1 ms. Other parameters for both sequences were: F O V = 24 cm, T R = 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 T 2 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 N N L S model described in Section 1.7.1.2. A regional T 2 analysis compared decay curve parameters estimated by a mono-exponential and N N L S from the M E S E 0 dataset and the corresponding slice of the 6 slice slab M E S E C 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 T E = 10 ms image (top) and the T E = 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 F O V 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 M E S E 0 sequence and each slice of the M E S E C 127 C H A P T E R 4. NON-180 R E F O C U S I N G P U L S E S 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 C H A P T E R 4. NON-180 R E F O C U S I N G P U L S E S 4 . 4 . MULTI-SLICE, MULTI-ECHO Figure 4.14: The T E = 10 ms image from the four average, 32 echo M E S E 0 sequence (top row). The six slices of the six slice slab M E S E C 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 C H A P T E R 4. NON-180 R E F O C U S I N G P U L S E S 4.4. MULTI-SLICE, MULTI-ECHO Figure 4.15: Myelin water maps calculated using the small norm N N L S solution based on data from the 6 slice, 42 echo pulse sequence M E S E C . 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 a n d global density) which is why the two phantoms on the right side were thresholded. 130 C H A P T E R 4. NON-180 R E F O C U S I N G P U L S E S 4.4. MULTI-SLICE, MULTI-ECHO C H A P T E R 4. NON-180 R E F O C U S I N G P U L S E S 4.4. MULTI-SLICE, MULTI-ECHO 4.4.3.2 R e g i o n a l A n a l y s i s 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 T 2 decay curve, the N N L S Q algorithm was applied to the decay curve. The short T 2 and medium T 2 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 M E S E C sequence, was always higher than the myelin water fraction calculated from the data acquired with the M E S E G sequence (Table 4.12). The decay curves were further assessed by comparing curves acquired from the two pulse sequences. One example T 2 decay curve is shown in Figure 4.17. In general, the T 2 decay curve acquired with the slab selected, multi-slice M E S E C sequence decayed faster than the T 2 decay curve acquired with the M E S E 0 sequence. The M E S E C sequence may have subtle differences in the acquisition method to accommodate the acquisition of multiple slices within the slab. 132 C H A P T E R 4. NON-180 R E F O C U S I N G P U L S E S 4.4. MULTI-SLICE, MULTI-ECHO N N L S Structure Sequence Short p Short T 2 Med p Med T 2 frontal_wm_l M E S E 0 (180) 47.258 ( 4.6%) 15.000 931.387 (91.5%) 79.337 frontaLwmJ M E S E C (180) 78.021 ( 8.0%) 36.109 863.296 (88.1%) 80.566 frontal_wm_l M E S E C (174) 83.473 ( 8.5%) 34.430 873.974 (88.8%) 81.122 frontal_wm_r M E S E 0 (180) 36.290 ( 3.5%) 15.000 977.200 (93.7%) 79.755 frontaLwm_r M E S E C (180) 138.347 (13.5%) 24.656 885.655 (86.5%) 83.284 frontal_wm_r M E S E C (180) 138.347 (13.5%) 24.656 885.655 (86.5%) 83.284 genu_wmJ M E S E 0 (180) 94.612 ( 9.2%) 22.819 908.342 (88.0%) 77.340 genu-wmJ M E S E C (180) 150.983 (15.2%) 37.115 837.326 (84.4%) 78.994 genu_wmJ M E S E C (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 M E S E C (180) 179.500 (17.8%) 35.845 798.541 (79.3%) 81.258 genu_wm_r M E S E C (173) 179.886 (17.8%) 35.065 799.187 (78.9%) 80.791 post_w'mJ M E S E 0 (180) 88.007 ( 8.5%) 17.971 942.529 (91.5%)- 92.041 post_wm_l M E S E C (180) 244.997 (24.4%) 35.921 686.287 (68.4%). 96.742 post_wm_l M E S E C (180) 244.997 (24.4%) 35.921 686.287 (68.4%) 96.742 post_wm_r M E S E 0 (180) 70.068 ( 6.7%) 15.291 808.996 (77.5%) 87.281 , post_wm_r M E S E C (180) 224.328 (22.3%) 40.148 589.941 (58.5%) 86.406 post_wm_r M E S E C (180) 224.328 (22.3%) 40.148 589.941 (58.5%) 86.406 mid_wm_l M E S E 0 (180) 177.964 (17.2%) 29.219 607.018 (58.8%) 82.458 mid_wmJ M E S E C (180) 215.597 (22.4%) 39.133 500.218 (52.0%) 79.342 mid_wmJ M E S E C (180) 215.597 (22.4%) 39.133 500.218 (52.0%) 79.342 mid_wm_r M E S E 0 (180) 168.491 (16.7%) 29.417 606.233 (60.0%) 85.955 mid_wm_r M E S E C (180) 228.307 (23.1%) 30.727 567.696 (57.3%) 86.775 mid_wm_r M E S E C (178) 229.005 (23.1%) 30.709 567.653 (57.3%) 86.831 perivent_gm_l M E S E 0 (180) 0.000 ( 0.0%) — 1161.671 (100.0%) 76.568 perivent_gm_l M E S E C (180) 0.000 ( 0.0%) — 1196.713 (100.0%) 68.768 perivent-gmJ M E S E C (172) 0.000 ( 0.0%) — 1200.621 (100.0%) 68.533 perivent_gm_r M E S E 0 (180) 17.381 ( 1.5%) 15.000 1173.212 (98.5%) 81.559 perivent_gm_r M E S E C (180) 0.000 ( 0.0%) — 1163.000 (100.0%) 73.643 perivent_gm_r M E S E C (180) 0.000 ( 0.0%) — 1163.000 (100.0%) 73.643 post_gm_l M E S E 0 (180) 0.000 ( 0.0%) — 1018.847 (100.0%) 86.499 post_gm_l M E S E C (180) 197.091 (19.1%) 42.066 745.337 (72.2%) 79.251 post_gm_l M E S E C (159) 139.822 (13.1%) 22.560 926.926 (86.9%) 79.282 post_gm_r M E S E 0 (180) 0.000 ( 0.0%) — 873.843 (83.6%) 84.274 post_gm_r M E S E C (180) 199.367 (19.3%) 42.229 607.816 (59.0%) 79.542 post_gm_r M E S E C (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 N N L S solution was calculated from the M E S E 0 and M E S E C 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 N N L S short and medium parameters were estimated from bins 10 < T 2 < 50 ms and 50 < T 2 < 120 ms, respectively. 133 C H A P T E R 4. NON-180 R E F O C U S I N G P U L S E S 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 , 4 t / « > ^ ^ . w - . , v . . J , - 1 1 1 50 100 150 200 250 300 350 400 450 200 m u a 1 1000 800 O Measured _ . MESE c Fitted (a=180) . . . . MESE Fitted (a=176.0) 600 400 20 40 8 » » p « 8 8 t e t a a a 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, T 2 decay curves acquired by the M E S E 0 and M E S E C (middle figure). NNLS fits to the M E S E C data, one assuming the re- focusing pulse flip angle was 180° and the other where the optimal flip angle was calculated (bottom figure). 134 C H A P T E R 4. NON-180 R E F O C U S I N G P U L S E S 4.4. MULTI-SLICE, MULTI-ECHO 4.4.3.3 Flip Angle Estimation The slice from the slab M E S E C that corresponded to the slice from the single-slice M E S E 0 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 M E S E C . 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 wil l 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 T 2 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 C H A P T E R 4. NON-180 R E F O C U S I N G P U L S E S 4.5. NON 90° EXCITATION PULSE maps calculated from M E S E 0 . The T 2 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 M R I pulse sequence contains two types of R F 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 (a e ) , the refocusing pulse (a r ) and the tissue parameters, p, T 2 , T i . 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 T 2 , T i , p and the pulse sequence parameters T E , a e and a r , 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) sin 2 (a/2) sin (a) 0 sin 2 (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 C H A P T E R 4. NON-180 R E F O C U S I N G P U L S E S 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 wil l assess the effect of having an inaccurate excitation pulse but a perfect refo- cusing pulse ar = 180°. M e t h o d s A set of true decay curves were calculated with the excitation pulse flip angle set' to ae — 9 0 ° , 8 5 ° , . . . , 5 0 ° . Other parameters were: ar = 180, 32 echoes, T i = 1000 ms, T 2 = 80 ms, T R = 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, T 2 and x2 (the misfit). Resu l t s and Discuss ion 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 s in(a e ) . 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 C H A P T E R 4. NON-180 R E F O C U S I N G P U L S E S 4.5. NON 90° EXCITATION PULSE True ae (°) ar (°) P T 2 ( ms) X 2 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 S imula t ion 2: Inaccurate Exc i t a t i on and Inaccurate Refocus ing, Assume Accura te Exc i ta t i on This simulation wil l 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. M e t h o d s The parameters for the simulation are: 32 echoes, T j = 1000 ms, T 2 = 80 ms, T R = oo, echo spacing = 1 0 ms and p — 1000. The excitation pulse was varied: ae = 90°, 8 5 ° , . . . , 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, T 2 and x2 (the misfit). The excitation pulse flip angle was assumed to be ae = 90°. Resu l t s and Discuss ion 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, T 2 decreased as ae decreased. The x2 misfit increased as ae decreased, therefore the curve fits got progressively worse as the excitation and 138 C H A P T E R 4. NON-180 R E F O C U S I N G P U L S E S 4.5. NON 90° EXCITATION PULSE refocusing pulses decreased. True ae/ar (°) ar (°) P T 2 (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°, 8 5 ° , . . . , 50° and ar = 2ae. The curve fitting calculated ar, T 2 , p and x2 assuming ae = 90°. The true T 2 = 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 T 2 estimation. This simulation wil l assess the effect of having ae = ar/2 and calculating both ae and ar along with the other parameters. M e t h o d s The parameters for the simulation were: 32 echoes, T i = 1000 ms, T 2 = 80 ms, T R = oo, echo spacing = 10 ms and p = 1000. The excitation pulse was varied: ae = 90°, 8 5 ° , . . . , 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 a r , ae, p, T 2 and %2 (the misfit) were calculated. The decay curve was fit assuming the ratio ae/ar = 1/2. Resu l t s and Discuss ion A l l parameters (shown in Table 4.15) were estimated accurately, including ae and ar. The stan- dard deviation of the estimated T 2 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 C H A P T E R 4. NON-180 R E F O C U S I N G P U L S E S 4.5. NON 90° EXCITATION PULSE True ae/ar (°) a S ( ° ) ar (°) P T 2 (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°, 8 5 ° , . . . , 50° and ar = 2ae. The curve fitting calculated ae, ar, T 2 , 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 T 2 , when the excitation pulse flip angle was assumed to be ae = 90°, was as accurate and precise (Table 4.14) as the T 2 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 C H A P T E R 4. NON-180 R E F O C U S I N G P U L S E S 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 N o n - C o n s t a n t T r a i n of F l i p A n g l e s The work previously showed that T 2 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 T 2 is of primary interest and is determined from the signal in the decay curve with T E < 50 ms. Therefore, the refocusing pulses with T E < 50 ms should be near 180° to retain the signal to noise, but the refocusing pulses where T E > 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, s a , 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 • sL i } =' i .o , . r . , i.o 16 • s i 2 ) = 0.5,. , 0.5 13 • s^3) = 1.0,1.0,1.0,0.5,0.5* . . ,0.5 141 C H A P T E R 4. NON-180 R E F O C U S I N G P U L S E S 4 . 6 . N O N - C O N S T A N T T R A I N • 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 . 1 9 ) , is an amalgamation of the first two, it retains 1 8 0 ° refocusing pulses for the first part of the train to retain the higher SNR and the last part of the train is 9 0 ° to decrease the power deposition. The fourth set of scalings (Figure 4 . 2 0 ) , S Q \ was defined as Sa\i) = e x p ( — i / 1 0 ) for i = 0 , 1 , . . . , 1 5 . 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 — - j x — / 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 C H A P T E R 4. NON-180 R E F O C U S I N G P U L S E S 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 I M I M I M l 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, T 2 and flip angle from a decay curve that was "collected" with non-constant refocusing pulses through the train. 4.6.1.1 M e t h o d s I created a set of simulations with p = 1000, T 2 = 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 Resul t s The parameters a, p, and T 2 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 C H A P T E R 4. NON-180 R E F O C U S I N G P U L S E S 4.6. NON-CONSTANT TRAIN S Q ^ compared to the constant refocusing pulse train of ar = 1 8 0 ° . As expected, p and T 2 had a (2) higher standard deviation in simulation s a . 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 M E S E C 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: F O V = 16 cm, T E = 12.12 ms, T R = 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 R F 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 N N L S 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 ŝ 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 C H A P T E R 4. NON-180 R E F O C U S I N G P U L S E S 4.6. NON-CONSTANT TRAIN a: 180 16 a: 90 16 800 ^ 6 0 0 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 ^ 6 0 0 CO c 5 400 CO c 200 0 800 ^ 6 0 0 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°, 9 0 ° , . . . , 90° (bottom, right). Measured data shown with whiskers and fitted data shown as circles. 145 C H A P T E R 4. NON-180 R E F O C U S I N G P U L S E S 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°, 9 0 ° , . . . , 90° (bottom, right). Measured data shown with whiskers and fitted data shown as circles. 146 C H A P T E R 4. NON-180 R E F O C U S I N G P U L S E S 4.6. NON-CONSTANT TRAIN The T 2 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 T 2 by approxi- mately 10%. Phantom Sequence T 2 (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 T 2 and refocusing pulse flip angle (for the first echo) calculated from each train of refocusing pulses. The T 2 , from the M E S E 0 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 a c a i c = 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 R F 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 R F pulses results is linearly correlated with the actual applied R F pulse flip angle. 147 C H A P T E R 4. NON-180 R E F O C U S I N G PULSES- 4.6. NON-CONSTANT TRAIN S90 x2 T 2 (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 R F pulses was 1.0,1.0,1.0,590 from ŝ 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 C H A P T E R 4. NON-180 R E F O C U S I N G P U L S E S 4.7. CONCLUSIONS 4.6.3 Conc lus ion The amplitude of the R F pulses in the train can be scaled to have any amplitude and the T 2 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 T 2 correction due to B i inhomogeneities have relied on the acquisition of a B i map or based on uniform phantoms. Other methods required two acquisitions to calculate the corrected mono-exponential T 2 , one acquisition to calculate the multi-echo decay curve and the second acquisition to estimate the B i field. The method I propose requires only one acquisition, the multi-echo decay curve is sufficient to calculate T 2 and a. Therefore, the T 2 and B i are inherently registered and the B i 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 B i fields, due to the dielectric effect [83], and require high power deposition for 180° refocusing pulses [84]. Gareau et al. [85] compared T 2 distributions of guinea pig brain at 1.5T and 4T and found both the short and medium T 2 peaks had a faster decay at the higher field strength. The decrease in T 2 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 T 2 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 T 2 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 M R I 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 M R I today. Image processing is an evolving field and there wil l 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 T 2 filter was to select for T 2 near 20 ms only, and suppressed all T 2 above 70 ms. The myelin water maps were comparable to myelin water maps calculated by the N N L S algorithm. The linear combination method to analyze T 2 from a multi-echo M R I dataset is restricted primarily to 150 C H A P T E R 5. C O N C L U S I O N S 5.2. LINEAR COMBINATION the short T 2 range. The restriction is due to the broad selection if longer T 2 are of interest. The broad selection would negate the use of linear combination for T 2 longer than 50 ms. Therefore, the curve fitting algorithms which provide information about the T 2 distribution are still required if the T 2 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 T 2 . Therefore, it may not be possible to create a better filter than what is presented. The short T 2 filtered associated with the coefficients calculated in this chapter is slightly asymmetric in that the filter for 20 < T 2 < 70 ms is less sharp than the filter for 10 < T 2 < 20 ms. This provides a difficulty that a myelin water compartment which has an increase in the T 2 is indistinguishable from a myelin water compartment that has only a decrease in the amplitude of the T 2 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 T 2 toward higher T 2 and is corroborated by results based on the N N L S algorithm. The shift toward higher T 2 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 N N L S technique, though, using fewer echoes may prove to be an important step. Shorter trains of R F 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 C H A P T E R 5. C O N C L U S I O N S 5.3. Bi/T 2 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 T 2 curve fitting algorithm should be able to estimate the T 2 distribution, as well as the refocusing pulse flip angle, from a multi-echo decay curve. I presented an algorithm to estimate T 2 and a from a mono-exponential decay curve and a modified N N L S algorithm, N N L S a , to estimate the T 2 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 T 2 . 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 T 2 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, T 2 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 T 2 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 T 2 decay curves, which have resulted in the Meiboom-Gill [86] extension to the original Carr-Purcell [87] pulses, and M L E V phase cycling [30]. The effectiveness of adjusting the phase of the R F pulse, rather than the amplitude, is that adjusting the phase can be done more precisely. The relationship between the amplitude of an R F 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 C H A P T E R 5. C O N C L U S I O N S the cause of some discrepancies seen in Section 4.6, tailored R F 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 T 2 is the pulse sequence must have a long T R 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 P re fe r red T2 Quant i f ica t ion P r o t o c o l The quantification of myelin water is important for continued basic science research. The current protocol, using M R I , 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 T 2 using the N N L S algorithm. The multi-echo sequence and N N L S 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 M R I dataset with a long T R . 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 B i 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 N M R 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, Apr i l 2000. [12] V . Vasilescu, V . Simplaceanu E . Katona and, and D. Demco. Water compartments in the myelinated nerve. III. Pulsed N M R results. Experientia, 34:1443-1444, 1978. [13] R M Kroeker and Henkelman R M . Analysis of biological N M R 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 . L i , 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. R S N A , 1991. [16] A . L . MacKay, K . P. Whittall , K . S. Cover, D. K . B . L i , and D. W . Paty. In vivo T 2 relaxation measurements of brain may provide myelin concentration. In Proceedings of the SMRM, page 917. S M R M , 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 C P M G 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. L i , 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 . L i , and D. W . Paty. In vivo measurement of T 2 distributions and water contents in normal human brain. Magn. Reson. Med., 37:34-43, 1997. [20] C Beaulieu, F R 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] I M Vavasour, K P Whittall , A L MacKay, D K L i , G Vorobeychik, and D W 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] G J Stanisz, A Kecojevic, M J Bronskill, and R M 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 N J , 1974. [29] D B Twieg. The k-trajectory formulation of the N M R 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 N M R spin-echo experiments. J. Magn. Reson., 43:65-80, 1981. [31] C S . Poon and R . M . Henkelman. Practical T 2 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 M R I 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 T 2 estimation using multislice multiple-echo imaging. Magn. Reson. Med., 4:34-47, 1987. [36] I. M . Vavasour, K . P. Whittal l , D. K . L i , and A . L . MacKay. Different magnetization transfer effects exhibited by the short and long T 2 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 C P M G relaxation rate in iron-rich gray matter. Magn. Reson. Med., 46:159-165, 2001. [39] C L . Lawson and R . J . Hanson. Solving least square problems. Prentice Hall , Englewood Cliffs N J , 1974. [40] K . P . Whit tal l and A . L . MacKay. Quantitative interpretation of N M R 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 N M R relaxation data. J. Magn. Reson., 95:221-234, 1991. [42] R. M . Kroeker and R. M . Henkelman. Analysis of biological N M R relaxation data with continuous distributions of relaxation times. J. Magn. Reson., 69:218-235, 1986. [43] K . P. Whittal l and A . L . MacKay. Quantitative analysis of N M R 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 T 2 relaxation data. Magn. Reson. Med., 35(3):370-378, 1996. 157 Bibliography [45] K . P. Whittall , A . L . MacKay, and D. K . B . L i . 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. A n M R I 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 N M R imaging. Magn. Reson. Imag., 4:263-266, 1986. [49] A . L . MacKay, K . P . Whittal l , J . Adler, D . K . B . L i , 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 M R I 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 M R 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 M R angiography by nonlinear anisotropic filtering. Magn Reson Med, 37(6):914-9, 1997. [55] R. Zingale and A . Zingale. Detection of M R I brain contour using isotropic and anisotropic diffusion filter. A comparative study. J Neurosurg Sci, 42(2) : l l l -4 , 1998. 158 Bibliography [56] S.M. Smith and J . M . Brady. S U S A N - 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 M R I 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 N M R . 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 N M R 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 T 2 using multiple-echo M R I 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 B I inhomogeneities. In Proceedings of the ISMRM, page 2071, 1997. [73] Zur Y . and Stokar S. A phase-cycling technique for cancellling spurious echoes in N M R 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 B I and BO variations in quantitative T2 measure- ments using M R I . 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 T 2 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 N M R fourier spectroscopy. J. Chem. Phys., 60(8):2966-2979, 15 Apr i l 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 T l of human brain tissue in vivo. Magn Reson Imag, 15(10):1133- 1143, 1997. [81] R Stollberger and P Wach. Imaging of the active B l field in vivo. Magn Reson Med, 35(2):246-51, 1996. [82] E . T. Lebsack and S. M . Wright. Iterative R F 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, C V . 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. Gi l l . 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., Whittal l , 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, Apr i l 21-27, (2001). Page 820. • Chapter 3: Myelin Water Fraction from Linear Combination - Paper: Jones, C.K., Xiang, Q.S., Whittal l , K .P . , Mackay, A . L . , "Short T 2 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 T 2 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 T 2 and B\ from Decay Curves Collected with non-1800 Refocusing Pulses", Proceedings of the International Society of Magnetic Resonance in Medicine. Toronto, O N , July 10-16, (2003). 162

Cite

Citation Scheme:

    

Usage Statistics

Country Views Downloads
France 23 1
United States 18 1
Canada 9 1
China 7 0
United Kingdom 5 0
Germany 2 2
Sweden 2 32
Singapore 2 0
Italy 2 0
Japan 1 0
City Views Downloads
Unknown 26 2
Champigneulles 11 1
Mountain View 8 0
Edmonton 7 0
Hangzhou 5 0
Bethesda 4 0
Stockholm 2 0
Vancouver 1 1
Toronto 1 0
Denver 1 0
Redmond 1 0
Tucson 1 1
Beijing 1 0

{[{ mDataHeader[type] }]} {[{ month[type] }]} {[{ tData[type] }]}

Share

Share to:

Comment

Related Items