T2 Decay Curve Acquisition and Analysis i n M R I Noise Considerations, Short T , and B 2 x Field Encoding by Craig K . Jones M . S c , The University of Western Ontario, 1997 B . S c , Simon Fraser University, 1992 A T H E S I S 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 O F THE REQUIREMENTS FOR T H E D E G R E E OF 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 UNIVERSITY 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 ( M W F ) 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 distribution and refocusing pulse 2 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. ii Table of Contents Abstract ii Table of Contents iii List of Tables vii List of Figures ix List of Algorithms xi Acknowledgements xii Glossary xiv 1 Introduction 1.1 1.2 1.3 1.4; 1.5 ' 1 Magnetic Resonance Imaging T and T i Relaxation !' 1.2.1 Spectral Density 1.2.2 Bloch Equations Central Nervous System Tissue Relaxation T Decay Curve 1.5.1 Mono-Exponential T Relaxation .1.5.2 Multi-Exponential T Relaxation . . T Decay Curve Acquisition 1.6.1 K-Space . . . 1.6.2 Optimized Multi-Echo Pulse Sequence 1.6.3 Standard Multi-Echo Pulse Sequence 1.6.4 Noise in Decay Curve 1.6.5 Artifacts Influencing the T Decay Curve T Decay Curve Analysis 1.7.1 Decay Curve Inversion (Decay Curve to T Distribution) 1.7.2 T Decay Curve Requirements Phantoms Overview of Thesis . . 2 2 2 2 1.6 2 2 1.7 2 2 2 1.8 1.9 in 1 2 2 3 4 67 7 8 9 9 10 11 14 16 18 18 24 25 26 Table of Contents 2 Myelin Water Fraction Noise Reduction 2.1 Introduction 2.2 SNR Effect on Myelin Water Fraction 2.3 Description of Spatial Filters 2.3.1 Anisotropic Diffusion Filter 2.3.2 Susan Filter 2.3.3 Median Filter 2.3.4 Wavelet Filter 2.4 Spatial Filters on Simulated Data 2.4.1 Methods 2.4.2 Results 2.4.3 Conclusions 2.5 Spatial Filters Applied to Zn Vivo Data 2.5.1 Data Acquisition and Analysis 2.5.2 Results 2.5.3 Conclusions 2.6 Spatial Filtering vs Averaging 2.6.1 Methods 2.6.2 Results 2.6.3 Discussion 2.6.4 Conclusions 2.7 Conclusions 3 . , v.. • • • ••'44. 46 :. . . . 50 53 53 Myelin Water Fraction from Linear Combination 3.1 Introduction 3.1.1 Linear Combination for T2 Selectivity 3.1.2 Properties of a Linear Combination 3.1.3 Overview 3.2 Calculation of Coefficients 3.2.1 Coefficients for Short T 3.2.2 Coefficients for All T 3.2.3 Description of the Algorithm 3.2.4 Results 3.2.5 Discussion 3.3 Linear Combination Simulations 3.3.1 Model 3.3.2 Results and Discussion 3.3.3 Conclusions 3.4 Selectivity of Mono-Exponential Phantoms 3.4.1 Methods 3.4.2 Results and Discussion 3.4.3 Conclusions 3.5 In Vivo Myelin Water Fraction: NNLS vs. T Filter 3.5.1 Methods 3.5.2 Scan Results 3.5.3 Conclusions 2 2 2 iv 27 27 28 29 30 31 31 32 . . 32 32 36 40 .40 40 41 44v .44 . 55 55 56 . . 57 57 58 59 60 '. . . . 61 62 63 64 64 64 65 66 66 66 68 68 69 69 72 Table of Contents 3.6 3.7 3.8 3.9 Calculation of Coefficients by Matrix Inversion 3.6.1 Calculation of Coefficients 3.6.2 Results and Discussion 3.6.3 Conclusions Calculation of Coefficients by Backus-Gilbert 3.7.1 Calculation of Coefficients 3.7.2 Results and Discussion 3.7.3 Conclusions Optimal Echo Times . 3.8.1 Filter Measures 3.8.2 Equally Spaced Echoes 3.8.3 Non-Equidistant 3.8.4 Discussion. . . 3.8.5 Four Echo Non-equidistant Sampling Conclusions 72 72 74 76 76 76 77 80 80 80 82 . . . . : 85 88 90 • • • • .-92 : . 4 Non-180 Refocusing 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 Decay Curve F i t 102 4.2.1 Mono-Exponential Fit Algorithm 102 4.2.2 Monoexponential F i t - Noise Simulation 103 4.2.3 Phantom Verification 105 4.3 Multi-Exponential T Decay Curve F i 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 2 2 : v Table of Contents 4.6 4.7 4.5.7 Conclusions Non-Constant Train of F l i p Angles 4.6.1 Simulations 4.6.2 Scanning 4.6.3 Conclusion Conclusions 141 141 143 144 149 149 5 Conclusions 5.1 5.2 5.3 5.4 150 Filtering Linear Combination B!/T Preferred T Quantification Protocol 150 150 152 153 2 2 Bibliography A Related Publications . . . . . . . . . 154 • 162 List of Tables 1.1 T i and T of mono-exponential nickel/agarose phantoms 2.1 2.4 2.5 2.6 The mean global error in myelin water fraction over the 10 noisy data sets for each filter The mean error in myelin water fraction over the 10 noisy data sets for each filter for the area of 15% myelin water The mean error in myelin water fraction over the 10 noisy data sets for each filter for the area of 30% myelin water In vivo myelin water fraction as a function of filter The number of times each filter had the smallest myelin water fraction . . . . . . . . Number of holes per R O I as a function of iterations of the anisotropic diffusion filter ; 38 ,41 42 48 3.1 3.2 3.3 3.4 3.5 The coefficients for the short T and all T linear combination method Mean and standard deviation of the signal in the short T phantom image . . . . . . Coefficients for the short, medium and long T selection based on matrix inversion. Coefficients calculated by the Backus-Gilbert method Coefficients for short and all T filter based on four unevenly spaced echoes. . . . . 62 67 74 77 90 4.1 4.2 Mono-exponential noise simulation as a function of a 104 Mono-exponential noise simulation as a function of a and S N R scaled by a function of the refocusing pulse . . . 104 Parameters of phantom derived from mono-exponential fit ,. 106 Bi-exponential noise simulation to assess accuracy and precision of parameter esti-.... mation from non-180 refocusing pulse train ;. 116 Noise simulation to estimate variability of multi-exponential decay curve parameters. 117 T time calculated from the M E S E sequence and the M E S E (with the flip angle calculated) .119 Refocusing pulse flip angle calculated from multi-echo decay curve •'. 119 Estimated T i from curve fits of five points ( T R = 100,300,800,2000,3000 ms) from a pulse saturation experiment .119 Percent myelin water estimated from scans with (x = 150° and a = 180°. . 121 Refocusing pulse flip angle calculated over the regions from the scans with <x — 150° and a = 180° .123 Geometric mean T of the medium component from select regions throughout the image 124 Decay curve parameters estimated from in vivo scan as a function of pulse sequence. 133 The excitation pulse was varied from 90° to 50° and the parameters were calculated assuming only the refocusing pulse was incorrect 138 2.2 2.3 4.3 4.4 25 2 2 2 2 2 2 v 0 4.5 4.6 4.7 4.8 4.9 4.10 2 Q C nom nom n o m 4.11 4.12 4.13 n o m 2 vii 37 37 List of Tables 4.14 The excitation and refocusing pulses were varied with a e = 90°, 8 5 ° , . . . ,50° and a = 2a 139 e r 4.15 The excitation and refocusing pulses were varied with a e = 90°, 8 5 ° , . . . , 50° and a = 2a 140 e r 4.16 Simulation to assess decay curve parameter estimation from non-constant refocusing pulse train (16 echoes) 4.17 The T and refocusing pulse flip angle (for the first echo) calculated from each train of refocusing pulses 4.18 The resulting x misfit given the scaling of the R F pulses was 1.0,1.0,1.0, S90 from 144 2 147 2 s<) 148 4 viii List of Figures 1.1 1.2 1.3 1.4 1.5 1.6 1.7 1.8 1.9 1.10 Myelinated axon •'• 5 Electron-micrograph of a myelinated cell 6 Simulated mono-exponential decay curves 8 Standard M R I image and corresponding k-space ' 10 Optimized single-slice, multi-echo pulse sequence ( M E S E ) 11 Conventional multiple echo 2D spin-echo pulse sequence ( M E S E ) . . . . . . . . . . 12 Conventional multiple echo 3D spin-echo pulse sequence ( M E S E ) . . . . . . . . . . 13 .: Rice distribution as a function of the mean . .14 Decay and T distribution of a bi-exponential white matter model . . . 19 Least squares and regularized T2 distributions of a white matter model . . . . . . 21 2.1 2.2 2.3 Histograms of myelin water fraction from simulated decay curves 29 Simulated myelin water fraction image annotated with six regions 33 Simulated myelin water fraction image (truth) from which the 32-echo data set was constructed 35 Myelin water maps calculated from spatially filtered, simulated multi-echo data. . 39 Myelin water maps calculated from spatially filtered M R I data 43 Myelin water maps from unfiltered and filtered (anisotropic diffusion filter) multiecho data .49 Histograms of unfiltered and filtered myelin water fractions . 51 . 2.4 2.5 2.6 2.7 3.1 3.2 3.3 3.4 3.5 3.6 3.7 3.8 3.9 3.10 3.11 3.12 3.13 3.14 3.15 3.16 0 C C 2 Constraints on the short T2 filter. Constraints on the "all" T filter T2 filters for short T components and all T components Myelin water fraction noise simulation to compare linear combination and N N L S . . Linear combination short T selectivity of phantoms Mean signal from phantoms in selected image versus T filter In vivo myelin water maps calculated by linear combination and N N L S Myelin water fractions calculated per region from the linear combination method and N N L S Bland-Altman plots for white matter regions and gray matter regions T2 filters based on coefficients calculated by matrix inversion Brain images T filtered from coefficients calculated by matrix inversion Coefficients calculated by Backus-Gilbert method T2 filters based on coefficients calculated by Backus-Gilbert method The measures of the effectiveness of a short T2 filter The coefficients as a function of the number of echoes and T E W i d t h of T2 filter for equally spaced echo times 2 2 2 2 2 2 ix 60 61 63 65 67 68 70 71 71 73 75 78 79 81 83 84 List of Figures 3.17 3.18 3.19 3.20 The coefficients as a function of the number of echoes and T E W i d t h of T filter for unequally spaced echo times The T filter, exp ( - T E / T ) as a function of T E The short T filter as a function of a consecutive subset of echoes based on the thirty-two coefficients 3.21 Short T filter calculated from four echoes compared to filter calculated from 32 echoes 3.22 Myelin water map calculated by linear combination from four echoes and from 32 echoes 2 2 2 86 87 89 2 89 2 4.1 4.2 4.3 4.4 4.5 4.6 4.7 4.8 4.9 91 92 Phase coherence diagram 96 Decay curves, collected with a — 180°, fit to obtain T and a 107 Decay curves, collected with a = 150°, fit to obtain T and a 108 Decay curves, collected with a = 110°, fit to obtain T and a 109 Decay curves, collected with a = 90°, fit to obtain T and a 110 The x misfit plotted as a function of refocusing pulse flip angle ' 112 The rank and inverse of the condition number of A based on a 32 echo decay . . 114 The rank and inverse of the condition number of A based on a 200 echo decay . . 115 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 prescribed as 180° and 150° 122 4.11 The arithmetic mean T maps calculated from the pulse sequence with a refocusing pulse flip angle prescribed as 180° and 150° . . 123 4.12 The x 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 sequence and six slices of the six slice slab M E S E acquisition 129 4.15, Myelin water maps calculated from a 6-slice slab M E S E . . . . . . . 130 4.16 Difference between the myelin water map calculated from the M E S E sequence and the myelin water map calculated from the M E S E sequence 131 4.17 In vivo decay curve from M E S E and M E S E 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 . . . . 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 = 85 ms 145 4.22 Decay curve fit of decay curves collected with non-constant refocusing pulse train, phantom with T = 23 ms 146 4.23 Measured and fitted data from non-constant train of refocusing pulses 148 p 2 p 2 p 2 p 2 2 Q a 2 2 G C C C C C 0 C 2 2 x List of Algorithms 1.1 4.1 4.2. 4.3 4.4 Small norm N N L S regularization Signal calculation from a train of refocusing pulses T i relaxation T relaxation . Phase evolution of magnetization . 2 xi 23 . 100 101 . . . . 101 101 : 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 P h D . M y 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." W i t h 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 P h D . The financial, academic and other support (Ken Bigelow and Vilmos Soti, in particular), during the past four years, has allowed me to complete this P h D . 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 P h D . 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 P h D 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 P h D 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. O n 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. K e n 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, K e n Nixon, who I met with Wednesday mornings for several years. The one question K e n 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. K e n Whittal. K e n 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 K e n was... xiii Glossary Abbreviation Meaning BG Backus-Gilbert CSF Cerebrospinal Fluid FOV Field of View (cm) FWHM Full-Width-at-Half-Max GM Gray Matter 'ME" " • - • Multi-echo MESE Multi-echo Spin-Echo Pulse Sequence MRI Magnetic Resonance Imaging MS Multiple Sclerosis MT Magnetization Transfer MWF Myelin Water Fraction NNLS Non-Negative Least Squares RF Radiofrequency ROI Region-of-Interest RMS Root-Mean-Square SAR Specific Absorption Rate SI Signal Intensity SNR Signal to Noise Ratio T Echo spacing (ms) TE Echo Time (ms) Ti Spin-Lattice Relaxation Time (ms) TR Repetition Time (ms) T Spin-Spin Relaxation Time (ms) 2 xiv Glossary WM Glossary White Matter xv Chapter 1 Introduction In magnetic resonance imaging, information about the local water environment is quantified by T i and T proton relaxation times and the density of the protons, p. The measurement process 2 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 to assess tissue and pathology (Section 1.4). The T decay curve 2 2 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 ) and the 2 spin-lattice relaxation time ( T i ) . 1 C H A P T E R 1. I N T R O D U C T I O N 1.2 1.2. T AND T i R E L A X A T I O N 2 T and T i Relaxation 2 The T relaxation time calculated from a multi-echo decay curve acquired using the M R I scanner 2 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 , therefore, a modified set of magnetization' 2 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 , ofthe spins through the spectral density function [2]: c J where (1.1) M = ^iry~T^ in) 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]: where 1/T, = k\jW(cj ) + J (2u )\ 1/T = k \jM(0) + 2 0 (1.2) {2) Jo -jU("o) + 5 (w) is defined in Equation 1.1 and k = (^f 2 \jW(2u ) 0 7 / i § J (I + 1) [2]. 4 2 (1.3) : 1.2. T AND Tj R E L A X A T I O N C H A P T E R 1. I N T R O D U C T I O N 1.2.2 2 Bloch Equations The magnetization, as a function of time, was defined in a set of phenomenological equations by Felix Bloch [1]: c? .. - M ( t ) = r 7 ,,/N -i M,.(/.) „ MM) „ [ M ( t ) x B ] - ^ - ^ i y + Mo T l MM) „ * where M ( i ) is the magnetization vector; B is the applied magnetic field; T ( 2 ,L 4 x ) 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 i n each orthogonal direction: j MM) = jM {t) = t y M , 5 7 l M B x 0 0 - ^ (1.5) - ^ (1.6) , ± ( MM 7) The solution to the coupled differential equations i n Equations 1.5,1.6,1.7 is: MM) =[M (0)cos(w <) M (t) = [-M (0)sia(u t) + M {Q)cos(u t)}exp(-t/T ) (1.9) MM) = (1.10) y x + Af (0)sin(woO]exR(-t/T2). 0 y x Q y 0 M + [M (0)-M ]exp(-i/T ). 0 z 0 1 2 . (1.8). where OJQ = jBo and Mo is the equilibrium magnetization. B y defining the transverse magnetization, M (t), xy as M (t) xy = M (t) + iM (t), x y 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: M {t) = M e ^ M (t) = M (l-e-'/ xy z 0 0 2 e (1.11) 1Wot T l ). (1.12) C H A P T E R 1. where M (0) xy INTRODUCTION 1.3. C E N T R A L NERVOUS S Y S T E M = 1 and M ( 0 ) = 0. The frequency of precession OJQ is typically demodulated out 2 of the transverse magnetization to obtain: M (t) xy M (t) z = M e-*/ (1.13) T 2 0 = M (l -e-*/ i) . (1.14) T 0 The measured signal intensity, y(ti) is a function of the magnitude of the transverse magnetization y(t) oc K\M (t)\ xy = K Mn e-*/ 2 T 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 relaxation, 2 therefore, the acquisition and analysis of T 1.3 2 decay curves will be discussed further. 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 between 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. INTRODUCTION 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 . Simple 2 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 exponentials (multi-exponential decay). 2 5 C H A P T E R 1. I N T R O D U C T I O N 1.4. TISSUE R E L A X A T I O N 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 analysis of frog sciatic nerve. They found three distinct T 2 2 decay curve 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. component, T 2 A n intermediate relaxing ~ 80 ms ascribed to axoplasmic water. A n d 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 relaxation 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 peak ( T ~ 15 ms) was attributed to water trapped between the myelin bilayers. 2 Menon fit T 2 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 DECAY CURVE 2 four distinct components: component near T = 1 ms attributed to phospholipid protons in the 2 tissue, component near T = 12.7 ms attributed to water i n the myelin layers, component near 2 T 2 = 89 ms attributed to intra and extra-cellular water, and a component not attributed to anything at T ~ 340 ms. The short component, T = 12.7 ms, contributed 6.8% to the total 2 2 water in the T distribution. Menon also noted the myelin bound water would be in slow exchange 2 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 component was by MacKay et al. [15, 16]. In vivo 2 T 2 relaxation measurements [17, 18, 19] distinguish three water compartments in the brain: water trapped between the myelin bilayers (10 < T < 50 ms), intra/extra axonal water ( T ss 80 ms) 2 2 and water associated with cerebrospinal fluid (CSF) ( T > 2000 ms). The fraction of the water 2 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 component. Beaulieu et al. [20] compared T distributions from myelinated and 2 2 non-myelinated nerves from the garfish and found the short T component only i n the myelinated 2 axon. More recently, many studies have compared the short T component to other M R I phe2 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 Relaxation 2 The T decay curve is the signal intensity measured from the radiofrequency coils as a function of 2 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) = \M (t)\ as defined in Equation 1.13. A xy 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 = 20 ms (circles), T = 80 ms (crosses), 2 2 and T = 2000 ms (diamonds). For each curve, y = pexp (—t/T ), the proton density p = 1000 2 2 and the time of the echoes was t = 10,20,... , 320 ms. 7 C H A P T E R 1. INTRODUCTION 1.5. T DECAY CURVE 2 1000r 150 TE (ms) Figure 1.3: Simulated decay curve with signal intensity S(t) = pexp (—£/T ) with p = 1000, t = 10,20,... , 320 ms and T = 20 ms (circles), T = 80 ms (crosses), and T = 2000 ms (diamonds). T R was considered infinite. 2 2 1.5.2 Multi-Exponential T 2 2 2 Relaxation Many tissues in the human body have more than one water compartment visible to M R I (Section 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 Ja (a < x < b) (1.16) 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 D E C A Y CURVE ACQUISITION 2 To solve integral equations on the computer, Equation 1.15 is discretized: M Y, 8(T-T ) s(T) = Sj j 3=1 where M is the number of T components and 6 (x) is the Dirac delta function defined as: 2 (1.18) y 0, otherwise. Substituting Equation 1.17 into Equation 1.15 gives: y(U) = ^2s exp{-U/T ) j j 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 G , G , and G ; the R F pulse train; and the acquired data. The x y z acquisition ofthe data and the relationship to the image is discussed (Section 1.6.1). Each pulse sequence will be discussed (Sections 1.6.2 and 1.6.3), as well as noise considerations (Section 1.6.4) and artifacts (Section 1.6.5). 1.6.1 K-Space Standard 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 (k , k ), is related by: x S(k ,k )= x y y / / M'(x,y)exp(i x k )exp(i y k )dy dx x 9 y (1.20) C H A P T E R 1. I N T R O D U C T I O N 1.6. T D E C A Y CURVE ACQUISITION 2 Figure 1.4: Standard M R I image of a brain (left) and the magnitude ofthe corresponding k-space image (right). where S (k , k ) is the signal at k , k and M (x, y) is the image space signal at x, y. The variables x y x y k and k describe the "distance" moved in k-space and are defined as: x y fc* = 7 / G (t')dt' x Jo k = i I G {t')dt' Jo y y (1.21) (1.22) where 7 is the gyromagnetic ratio, Gx(t') and Gy(t') are the gradient shapes, t is the duration of the gradient. Each point i n 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 measurements were calculated. The implementation 2 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. INTRODUCTION 1.6. T D E C A Y CURVE ACQUISITION 2 90 ) where the duration of the 90° hard pulse was 300 fis and the duration of the 180° hard x 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 will be annotated as MESE„ in this thesis. / RF G / 11 11 A A A A v V I1 i - A A A •' A- A ... A V V t\z v~ DACQ 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 ) . 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 were used to remove signal from stimulated echoes and signal out-of-slice. 0 2 1.6.3 Standard Multi-Echo Pulse Sequence The standard multi-echo pulse sequence (Figure 1.6) consists of a selective excitation pulse followed 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 G and G were used to encode the data acquisition in the 3D y z 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. INTRODUCTION 1.6. T D E C A Y CURVE ACQUISITION 2 number of slices and s is the thickness of each slice. In the 3D implementation used in this thesis, the phase encodes along G were added to the constant gradient crushers along G . The large, z 2 constant gradient crushers along G were applied to reduce signal from out-of-slice spins, but did z not reduce the stimulated echo artifact. This pulse sequence will be denoted by M E S E in this C thesis. , v J n • 1 yv nuf 1 yv (1 -4 (1 D n a-- n j -5 0 5 10 15 yv Jv -Jl -Jl II (L (1 n u n n n .j V.. 20 25 Time (ms) II' 30 35 40 Figure 1.6: Conventional multiple echo 2D spin-echo pulse sequence ( M E S E ) . pulse and refocusing pulses are slice selective. C 12 45 50 The excitation C H A P T E R 1. I N T R O D U C T I O N 1.6. T D E C A Y CURVE ACQUISITION 2 -A- RF — -4 Balanced Z Phase Encode: Slab Selective A Balanced Z Phase Encode: Balanced Z Phase Encode: DACQ 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 ) . The excitation pulse is slab selective (gray bar centered around 0 ms) and signal is phase encoded along G and G . The slice phase encodings are added to the large crushers in a way analogous to the phase encoding already on G ^ . C y z 13 C H A P T E R 1. 1.6.4 INTRODUCTION 1.6. T D E C A Y CURVE ACQUISITION 2 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: (1.23) -PRice(z) = where a, is the "standard deviation" of the function, / i is the mean, x is the parameter, and I (x) is a modified Bessel function of the first kind. The two interesting implications of the Rice 0 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. 2O(A«) = INTRODUCTION 1.6. T D E C A Y CURVE ACQUISITION 2 1 when JJL = 0 and exp ( - ^ f ) = exp ^ — w h e n |t/ = 0. Therefore, 2 -PRice(z) = -^exp = (1.24) '2^2 (1.25) ^fcayleigh (^) As p/a —> oo, the Rician distribution becomes similar to a Gaussian distribution. The asymptotic expansion [34] of the Bessel function IQ(X) when x is large is: I6{x) V2 7TX then for sufficiently large x, P(x) + = 77!= — J_ fe 1•9 2!(8x) + + 1-9-25 3!(8x) 3 (1.26) + • Therefore, considering P R i c e ( £ ) around a; = fi: x exp (J, + x x x + ^i 2 2a 2 (1.27) 2 2 exp 2 2 2<7 a 7 2 exp v^Ix exp (S) (1.28) (1.29) 2a 2 (1.30) then, because x ^ fi the scaling factor in front of the exponential can be simplified leaving: P(x) = — V2ira 2 •. exp -^Gaussian (1.31) 2a 2 (1.32) (^0 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 y {U) = yj[y(ti) + e i ] + e?,, where y is the true signal, and e\ and e are random numbers from a 2 e 2 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. 1.6.5 INTRODUCTION 1.6. T D E C A Y CURVE ACQUISITION 2 A r t i f a c t s Influencing the T D e c a y C u r v e 2 Artifacts in the T decay curve are non-random changes to the decay curve during data acquisition. 2 The artifacts could result from mis-tuned hardware or limitations of the underlying assumptions. Four effects that arise during acquisition of the T decay curve are: 1) signal from out-of-slice, 2) 2 magnetization transfer, 3) stimulated echoes, and 4) T dependence on the echo spacing r. 2 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 G z 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 carefully. For slice-selective imaging, a magnetic field gradient is applied along G which spatially z 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. INTRODUCTION 1.6. T D E C A Y CURVE ACQUISITION 2 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 will continue dephasing, a portion will begin to rephase and a portion will be stored along the longitudinal axis. The magnetization along the longitudinal axis does not accrue phase, therefore, when the magnetization along the longitudinal axis is rotated back into the transverse plane, it, will not have the same phase as the magnetization that was in the transverse plane. The lack of phase coherence results in magnetization refocusing at a different echo time. The stimulated echo artifact is described further in Chapter 4 T Dependence on Inter-Echo Spacing 2 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 ( T E = 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 1/T ) as R 2 2 = R 2 a + A R . The relaxation rate R 2 2 (R = 2 was defined independent of microscopic 2 a field inhomogeneities, and A R was derived to be a function which was approximated as being 2 proportional to re / where x = ADT/R , 3 2 2 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 is a function of T, the measured R would 2 2 also be a function of r . R is several orders of magnitude smaller than the measurement error for 2 the typical r of 5 ms. Therefore, comparing T from different techniques, requires knowledge of 2 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 1.7. T D E C A Y CURVE ANALYSIS 2 Decay Curve Analysis Decay curve data acquired by either technique described i n Section 1.6 must befitto 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 decay curve of y hort(U = 200exp•(—t/20), the 2 s medium T decay curve of y ediumW = 800 exp (—1/80), and the total bi-exponential decay curve 2 ybiexp(i) = 200 exp m {-t/20) + 800exp (-t/80). Echo times at t = 1 0 , 2 0 , . . . , 320 ms. The goal of any data inversion technique is to calculate the single or multiple T times and 2 amplitudes, s, for a given decay curve y. The bottom plot of Figure 1.9 shows the amplitudes, Sj, and T times of the T decay curves in the top plot of Figure 1.9. In this example, s\ = 200 at 2 2 T i = 20 ms and s = 800 at T 2 ) 2 22 = 80 ms. The plot of s versus T is called the T distribution. 2 2 Each data inversion technique may assume positions for the spikes in the T distribution, upper 2 and lower bound for the positions, or non-negativity of the amplitudes. Most inversion techniques minimize a misfit based on the measured data. 1.7.1.1 Misfit Measures The misfit is a scalar measure of the difference between the predicted data y and the measured data y . Two common definitions for the misfit between the predicted and measured data are Li, which measures | | y - y | | i = Y,iLi \Vi ~ Vii a n d 2 which measures ||y - y | | = 2~2iLi \Vi ~ Vi?• L 2 Another misfit measure, x , accounts for the difference between the predicted and measured 2 data, as well as the noise in the measurement, is defined as: (1.33) where N is the number of data points, yi is a measured datum, ai is the standard deviation of the i th is a reconstructed datum, and datum. The x misfit is a logical misfit to minimize in 2 curve fitting algorithms, as it accounts for noise in the measurements. 18 C H A P T E R 1. INTRODUCTION 1.7. T DECAY CURVE ANALYSIS 2 900 TE (ms) 800r ••. 700 k 600H & 500 CO c CD ° o o 400 - £ 300 200h 1 0 0 " ; ; : I j i Q\ 1 1 1 20 40 60 i i 80 ! : i i 100 120 T (ms) j • i i 160 140 j i 180 J 200 2 Figure 1.9: Decay curves (top) and T distribution (bottom) of a bi-exponential white matter model. The short T peak at 20 ms in the T distribution (bottom) corresponds to the fast decaying exponential (exp (—i/20), dashed-dotted line in the top plot). The medium T peak at 80 ms in the T distribution (bottom) corresponds to the slower decaying exponential (exp (—1/80), dash line in the top plot). The two peaks in the T distribution (bottom) correspond to the sum of the two decaying exponentials and is the solid line in the top plot (200exp ( - i / 2 0 ) + 800exp (-t/80)). The echo times are t = 10, 2 0 , . . . , 320 ms. 2 2 2 2 2 2 19 C H A P T E R 1. INTRODUCTION 1.7. T D E C A Y CURVE ANALYSIS 2 The expected value E(x ) = N and the standard deviation a 2 = y/2N. Solutions with x 2 x 2 ^ N reproduce the noise and solutions with x 3> N misrepresent the measurement. Intuitively, 2 if .the only difference between the predicted and measured data is noise, then Equation 1.33 is EiLi = £ £ i i = N. 1.7.1.2 Non-Negative Least Squares The primary curve fit algorithm used i n this thesis was a non-negative least squares algorithm. The T distribution is calculated from a decay curve by applying the non-negative least squares 2 (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/T j ) , yi is the datum at ti, and s is the set of amplitudes over T . Therefore, 2 2 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 misfit in Equation 1.33. 2 N N L S is guaranteed to minimize Equation 1.34 [39], or the x ; if the standard deviation is 2 divided through both sides As = y. It is interesting to note that each column of A is a monoexponential decay curve with T = T 2 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 (Equation .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. INTRODUCTION 1.7. T D E C A Y CURVE ANALYSIS 2 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 0 1 -1 0 0 0 : -1 (1.36) and f — 0. Or the matrix could be designed to minimize the amplitudes as in: 1 0 H= 0 ... o 0 1 0 ... o 0 0 : 1 0 0 0 ••• 0 1 (1.37) and f = 0. 800 700 600 ^500 c D 400 co o </> CD CL 300 200 100 20 40 60 80 100 120 T (ms) 140 160 180 200 2 Figure 1.10: T 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 distribution. 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). 2 21 C H A P T E R 1. INTRODUCTION 1.7. T D E C A Y CURVE ANALYSIS 2 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: a an 1 1 2 a.2i aN 1 <Z3iV 1 2 032 Si -. r 2/1 S2 — V2 (1.38) SN C-M2 AMI •• 1 O-MN VM B The set of T was defined to be 120 logarithmically spaced T = 10.00,10.52,11.06,... , 4000.0 ms. 2 2 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 x solution, Xmim 2 then regularizing the T distribution until a solution was found with 1.02 Xmin — X < 1-025 Xmin2 2 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 1.7.1.2.1 2 is described in Algorithm 1.1. T Distribution Bins 2 The T distribution calculated using N N L S from a tissue sample is expected to have broad distri2 butions due to the underlying heterogeneity of the tissue [42] and the regularization method [42, 43]. To quantify the T relaxation time or amplitude, the T 2 2 distribution is typically binned. A bin is defined with a start T and end T . 2 2 There are several measures defined for a bin. The most important measure, in particular for the short T , is the fraction of signal in the bin relative to the total signal: 2 (1.39) fr = 22 C H A P T E R 1. INTRODUCTION 1.7. T D E C A Y CURVE ANALYSIS 2 Algorithm 1.1 Small norm N N L S regularization. Regularized N N L S 1: Define f = 1.02 { = fiowXmin} 2 2: Define f = 1.025 {x high = fhighX in\ {Calculate the least squares solution} 3: Solve s = nnls (A, y) 2 low X low 2 h i g h m 4: Calculate Xmin fr° A y> 5: Define a vector, /u, logarithmically spaced between 10~ and 1 0 {Find m such that x\ between x20W and Xhigh } 6: while x < X O R x > xligh d o 7: for all m in \x do 8: Solve s = nnls (A, y, Hi) 9: Calculate xj from A, y, s m s 5 _ J S 2 2 2 0W 10: end for 11: 13; 14: 15: Do linear fit of x vs fj, From linear fit, calculate fi at midpoint between Solve s = nnls (A, y, Hi) Calculate x from A, y, s if X ? < X L 16: Set H to be linear spaced between Hi and M(1) 12: 17: 18: 19: 20: 2 xf 2 t h e n else if Xi > xligh t h e n Set H to be linear spaced between /i(end) and Hi else Do Nothing end if 22: end while 21: 23 a n ow d x\i g C H A P T E R 1. I N T R O D U C T I O N 1.7. T DECAY CURVE ANALYSIS 2 The arithmetic mean T of a bin is defined as: 2 s{T) T E T-T T=T tart T—T i E T=T„tart , s(T) end T = (1.40) s ent and the energy for the bin is £ start For this thesis, two compartments were defined. The short T compartment (related to water 2 trapped between the myelin bilayers, see Section 1.4) was defined as a bin with 10 < T < 50 ms. 2 The medium T bin was defined as 50 < T < 120 ms and is related to the inter and intra-cellular 2 2 water. 1.7.2 T Decay Curve Requirements 2 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 compo2 nents of one order of magnitude apart in a T distribution. The standard method to attain high 2 SNR, in the M E S E G sequence, is to acquire four averages for each line acquired i n k-space [18]. The S N R is proportional to VN where N is the number of averages - therefore, an image acquired from four averages has twice the S N R 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 i n the scanner for one sequence. Therefore, other methods of attaining the required: S N R must be assessed. 1.7.2.2 T Decay Curve Sampling 2 The T decay curve must be sampled sufficiently to accurately estimate the T components and 2 amplitudes. Whittall 2 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 ( T E = 30, 60, 90,120 ms) missed the short T component, underestimated the T , and 2 2 overestimated the proton density. Multi-exponential fits of the same simulated data sampled at 24 C H A P T E R 1. INTRODUCTION 1.8. PHANTOMS 32 echoes ( T E = 10, 2 0 , . . . , 320 ms) accurately reproduced the two T components. Whittall et 2 al. compared mono-exponential T solutions from simulated decay curves sampled with different 2 sets of echo times. They found the T and proton density calculated from two different sampling 2 schemes to be statistically different even though the underlying T and proton density were the 2 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 and amplitudes. 2 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 ( N i C l ) was chosen, as the paramagnetic ion N i is an 2 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 measurements were done with a 48 echo M E S E pulse sequence with T E = 10,20,... , 480 2 0 ms. Other parameters were T R = 3000 ms, 2 averages, and F O V = 16 cm. Phantom 1 2 3 4 5 6 T (ms) 2 237 91 98 24 84 26 T i (ms) mM NiCl 712 712 i858 333 342 690 2 2 0.5 5 5 2 2 % Agarose 0.2 1 1 4 1 4 Table 1.1: T i and T of mono-exponential nickel/agarose phantoms. 2 25 C H A P T E R 1. 1.9 INTRODUCTION 1.9. OVERVIEW OF THESIS Overview of Thesis Noise and artifacts which result from decay curve acquisition can confound the analysis and interpretation of T decay curves. The thesis focused on three aspects of the noise and artifacts 2 as they relate to T decay curve analysis. 2 First, the curve fitting algorithms used to fit T S N R of the acquired data. decay curves are strongly affected by the 2 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 decay curve is 2 to use a curve fitting algorithm. The standard algorithm is slow and problematic for low S N R multi-echo data. The short T is of particular interest as it relates to the myelin water fraction. 2 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 curvefittechnique. Third, it was previously shown that stimulated echoes encoded over the primary echoes will lead to inaccurate T estimation. Chapter 4 discusses a method to estimate the T and refocusing 2 2 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 sequence to validate the theory. Myelin water fractions 0 and maps estimated from the M E S E decay curves were compared to those estimated from the C M E S E decay curves. Several extensions are discussed. 0 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 i n 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 ~ 80 ms) and water associated with cerebrospinal fluid 2 (CSF) ( T > 2000 ms). The multiple water compartments in normal white matter result in a 2 multi-exponential M R signal decay as a function of echo time (TE) and the shortest T component 2 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 ( T distribution), was calculated voxel-by2 2 voxel from the T decay curves by a regularized non-negative least-squares (NNLS [28]) algorithm. 2 A high signal-to-noise ratio (SNR) is required to distinguish components in the T distribution [44] 2 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 N O I S E R E D U C T I O N 2.2. SNR AND M W F 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. 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 M S 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 spinecho 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 o n 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 (Section 1.7.2:1). The effect of S N R on myelin water fraction will be shown by a noise simulation study. The variability of the myelin water fractions, as a function of S N R , was calculated from simulated decay curves. A 32-echo decay curve was created with T E = 1 0 , 2 0 , . . . , 320 ms, a short water compartment with p = 200 and T s T = 80 ms and p m 2 m S 2 = 20 ms, and a medium water compartment with = 800. One thousand realizations of noise were added in quadrature for each of S N R = 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 N O I S E 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 S N R < 100 and the shape of the histogram was non-Gaussian. When the S N R 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 SNR 50 - - • SNR 100 - - SNR 200 SNR400 200 150 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 S N R 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 S N R images:; averaging during data acquisition and noise reduction image processing. Data acquisition averaging is a standard technique in M R I to obtain higher S N R 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 N O I S E 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 i n 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): D = I (x + l,y)-I (x,y) D w = I (x, y) - I (x - 1, y) D N = I (x, y + 1) - I (x, y) D = I (x,y) - I {x,y - 1) E e e s (2.1) e e e e e e where I (x, y) was the signal intensity at spatial location (x,y) i n the e th e echo image. Then the signal intensity at (x,y) was updated: I'{x,y) = I{x,y) + L\{c D -cwD +c D -csD ) E E w N N (2.2) s 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: i ^ 2 N Cd(x,y) = exp < - 30 (2.3) C H A P T E R 2. M W F N O I S E 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: >/V/iC?) + V / ( ^ ) + ... + V 7 ( ^ ) ' 2 2 exp < - 2 2 2 (2.4) M 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( ) - Ui,J)^0,0)^ X r V + ^y ^ ^{ Z.(ij)^(o,o) where r = \fi 2 2.3.3 + e x e P \ ^ > 2^ \ ——: . + j , a controls the scale of the smoothing, t is the brightness threshold. 2 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 N O I S E R E D U C T I O N 2.3.4 2.4. SPATIAL FILTERS ON SIMULATED DATA 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 (A') levels 2 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 N O I S E R E D U C T I O N 2.4. SPATIAL FILTERS ON SIMULATED DATA A C B 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 N O I S E 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 0 (2.6) 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) 0.15[c - s(r)]/[117 - s(r)] ~ ~ s{r) < c < 118 (2-7) 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 (2.8) s(r) < c < 256 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 N O I S E 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 S N R =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 ( a n i s o l ) , 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 N O I S E 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 compared to the true myelin water fraction. The measure was defined as: (2.10) 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 such that it did not come within 5 pixels of the border, so the boundary does not create a bias 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 36 C H A P T E R 2. M W F N O I S E R E D U C T I O N Filter 2.4. SPATIAL FILTERS ON SIMULATED DATA Global Error Measure none anisol aniso3 aniso5 median susan postmedian wavelet46 wavelet48 wavelet50 wavelet52 1.964 1.208 0.758 0.637 1.175 0.741 1.340 0.723 2.504 0.818 1.048 (0.005) (0.009) (0.012) (0.011) (0.014) (0.013) (0.010) (0.016) (5.809) (0.084) (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 none anisol aniso3 aniso5 median susan postmedian wavelet46 wavelet48 wavelet50 wavelet52 Error Measure 2.623 1.773 1.038 0.813 1.584 0.928 1.719 1.021 2.340 1.208 1.570 (0.081) (0.081) (0.061) (0.056) (0.075) (0.054) (0.086) (0.052) (4.329) (0.162) (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 N O I S E 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 none anisol aniso3 aniso5 median susan postmedian wavelet46 wavelet48 wavelet50 wavelet52 Error Measure 2 891 1 927 1 247 0 957 1 733 1 124 1 326 1 270 4 030 1 499 1 819 (0 (0 (0 (0 (0 (0 (0 (0 (8 (0 (0 065) 079) 113) 123) 087) 091) 082) 079) 903) 147) 080) Table 2.3: The mean error i n 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 N O I S E 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 N O I S E R E D U C T I O N 2.4.3 2.5. SPATIAL FILTERS ON IN VIVO DATA 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 ( M E ) , 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 C S F 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 ( a n i s o l ) , 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 N O I S E 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 none anisol aniso3 aniso5 susan medfilt postmed wavelet20 wavelet25 wavelet30 fWM 9.6 8.9 8.1 7.8 8.5 8.7 8.8 8.5 8.6 8.8 (4.8) (3.6) (2.2) (2.0) (3.9) (3.0) (2.2) (3.4) (3.4) (3.6) genu 15.7 15.5 15.3 15.2 15.5 15.2 15.8 15.1 15.3 15.4 postWM (4.1) (3,6) (3.5) (3.0) (3.6) (3.9) (2.9) (4.3) (4.2) (4.2) 12.1 11.7. 12.0 11.7 12.3 12.2 12.2 12.0 12.0 12.2 (4.1.) (3.3) (2.6) (2.4) (3.1) (3.7) (2.6) (3.4) (3.2) (3.6) postIC 12.1 12.2 12.9 13.0 12.1 12.9 11.5 12.4 12.6 12.3 (5.7) (6.2) (6.8) (6.8) (6.0) (6.6) (4.6) (5.0) (5.7) (5.8) spl en 15.0 14.9 14.7 14.4 15.1 16.0 15.2 14.8 15.0 15.1 (5.2) (4.8) (4.0) (3.4) (5.0) (3.8) (3.0) (5.8) (5.5) (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 N O I S E 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 anisol aniso3 aniso5 susan medfilt postmed wavelet20 wavelet25 wavelet30 9 0 0 0 0 35 0 1 50 0 4 33 5 o 0 1 17 1 0 0 Table 2.5: The number of times each filter had the smallest myelin water fraction standard deviation over all the regions (middle column) and over the white matter ( W M ) regions (last column). 42 C H A P T E R 2. M W F N O I S E 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 N O I S E R E D U C T I O N 2.5.3 2.6. SPATIAL FILTERING VS A V E R A G I N G 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 ( M E ) 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 N O I S E R E D U C T I O N 2.6.1.2 2.6. SPATIAL FILTERING VS AVERAGING 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 S N R 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: (2.11) where s ( T ) was the amplitude of the T distribution calculated by N N L S and the maximum T 2 2 2 was 2000 ms. Two sets of analysis were performed: 1) M W F s calculated per R O I 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 R O I : 1. Average the decay curves in the R O I and apply the N N L S algorithm to the averaged decay curve then calculate the M W F ( a v e r a g e - i n v e r t ) . 2. Apply the N N L S algorithm to each decay curve in the R O I , 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 a v e r a g e - i n v e r t , R O I 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 distributions calculated voxel-by-voxel by the N N L S algorithm (no averaging 2 of neighboring decay curves). The number of "holes" in an R O I of the myelin water map was defined as the number of voxels in a white matter R O I 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.) 45 C H A P T E R 2. M W F N O I S E R E D U C T I O N 2.6.1.3 2.6. SPATIAL FILTERING VS AVERAGING Bootstrap Theory The bootstrap algorithm [62] is a means of estimating confidence intervals and other statistical measures on a set of data. Consider an R O I 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 R O I . 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 R O I . The standard deviation is an estimate of the standard error of the M W F mean. 2.6.2 2.6.2.1 Results 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 R O I 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 R O I . 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 N O I S E R E D U C T I O N 2.6. SPATIAL FILTERING VS AVERAGING standard" a v e r a g e - i n v e r t , f i. The statistic D — h , f ai a 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 R O I the M W F s calculated from the b o o t s t r a p method were found to be statistically the same as those calculated by the a v e r a g e - i n v e r t method (minimum p=0.11). For the rest of the work the mean M W F within an R O I was calculated using the b o o t s t r a p method unless otherwise stated. 2.6.2.3 2.6.2.3.1 Myelin Water Fraction 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 b o o t s t r a p 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 i n 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 R O I 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 N O I S E R E D U C T I O N 2.6. SPATIAL FILTERING VS AVERAGING water fraction therefore it was important to compare i n v e r t - a v e r a g e 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 iteration), 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 2 4 20 (8) 15 (9) 8(6) 19 (9) 14 (9) 6(5) 17 (9) 13 (9) 4(3) 14 (9) 12 (10) 3(3) Table 2.6: Number of holes per R O I 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 N O I S E 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 N O I S E R E D U C T I O N 2.6.3 2.6. SPATIAL FILTERING VS AVERAGING 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 times of the peaks in the calculated T 2 2 distribution) and the standard deviation decreased for each R O I 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 homogeneous 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 R O I measurement compared to the standard. I expected that the mean M W F calculated from the 2-average, multi-echo data with 5 iterations 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 R O I (for example, see Figure 2.7). The M W F in the 4-average, unfiltered R O I 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 N O I S E R E D U C T I O N 2.6. SPATIAL FILTERING VS AVERAGING 0.2 0.4 I . 1 0.2 1 Average 0.2 2 Averages 0.2 4 Averages 0.4 Figure 2.7: Histograms of myelin water fraction calculated from the 1-average (first column), 2average (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 N O I S E 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 R O I 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 S N R =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 R O I but did not shift it lower. I believe that this was a result of using the non-linear curve fitting technique on low S N R data. The invert-average method has been used to create the myelin water maps and there- fore is important to validate compared to the standard. Even though there was a large rel- ative difference in the mean M W F (16% for 5 iterations for the 4-average data) between the invert-average the unfiltered and bootstrap invert-average methods, the relative difference decreased from over 29% for to 16%. The M W F calculated from the invert-average method tended to be closer to the standard as more iterations of the anisotropic filter were applied to the original 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 invert-average method compared to the bootstrap method (on the 4-average data). This apparent increase in the myelin water fraction was due to the non-linear N N L S algorithm applied to lower S N R data in the invert-average 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 N O I S E 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 S N R of the multiecho data. The current method to acquire sufficient S N R is to collect multiple averages of the multi-echo data, but this results in a long scan time. A second approach to have higher S N R 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 echo data. For the in vivo multi- 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 N O I S E 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 S N R to obtain similar results to the 4-average data. The minimum lower bound on the S N R 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 Myelin Water Fraction from Linear Combination 3.1 Introduction Many tissues in the human body consist of several water compartments (or components) within the same voxel. White matter, for example, has two compartments, a short compartment with T 2 ~ 20 ms and medium compartment with T ~ 80 ms [18]. The ratio of the areas of the short 2 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 fraction, in vivo, is to collect a decay curve with 2 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 stimulated echoes. A regularized N N L S algorithm minimized the total signal of the T distribution, 2 such that l.02x i m n < X 2 < ^-•^^Xmin where Xmin w a s the minimum misfit of the unregularized solution. An 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: (3.1) where A ( T ) is the filter of interest, c is a vector of coefficients. The image corresponding to the 2 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 7 (x) = ^ c , / ( x ) L (3.2) i i=i where 7i(x) = 2~^j=i P j ( ) x e x P ( f /T j(x)). The goal, then, is to define coefficients, c, such that — 2; A ( T ) includes T between 10 and 50 ms and excludes all other T . The method to calculate the 2 2 2 coefficients c is the novel work in this chapter. 3.1.1 Linear Combination for T Selectivity 2 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;/T (x) 2j and L the amplitudes of the T distribution. This method is described further in Section 3.6.1. 2 3.1.1.2 Calculation of Coefficients by Backus-Gilbert Method In 1991 Whittall [65] used the Backus-Gilbert technique [66, 67] to calculate the coefficients. T o 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 ) , and therefore allows a better selection. 2 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.2 3.1. INTRODUCTION 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) -* ( x ) (3.3) s A/(x) A (x) (3.4) 5 • 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) + / (x) ^ 2 A system that satisfies ( x ) + <? (x) f f l (3.5) 2 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 component. The nonlinear algorithm allows a larger variety of constraints on the filter 2 that were not possible using previous methods. As my primary focus for the work was to calculate 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 , c , and a second set that selects all T , c . To verify the s a 2 2 linear combination was selecting the appropriate T toms of known T i / T 2 2 components, a set of nickel/agarose phan- 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. CALCULATION OFCOEFFICIENTS 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 T 2 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 of interest and zero elsewhere. Coefficients for both filters, 2 J 4 ( T ) and A ( T ) , were calculated in a similar manner, but each with slightly different sets of S a 2 2 constraints applied to the desired filters. The algorithm described below was used to calculate the two sets of coefficients: one set, c , that selects for 1 0 < T 2 < 5 0 ms and the second set, c , s a that selects T 2 > 1 0 ms. The selection filter for the short T coefficients was defined as 2 A s (T2) = ^2iLi t c e x P (~^i/T ) 2 (and similarly for the total T ) . The desired filter T ! ( T ) is a boxcar: s 2 2 1 T lower < To < T o A (T ) = { _ S 2 0 u p p e r 2_ 2 ( 3 7 ) elsewhere and for j4 (T ): a 2 , A (T ) = { a 1 T 2 0 T 2 2 > 1 0 ms ~ < 1 0 ms (3.8) though neither definition is able to be constructed from a finite set of exponential decays. Similarly, the desired filter for ^4 (T2) is 1 for T a 2 > 1 0 ms (a step function). Therefore, the goal was to create a filter, J 4 ( T ) , as broad as possible and near 0 for T 2 > 8 0 ms (to suppress the signal 5 2 from the large peak at T 2 ~ 8 0 ms). Let a T 2 weighted image Ii be defined as: M 7i(x) = kp(x) J ^ S j - e x p [-<i/T -(x)] + a, (x) i=i 2j (3.9) where p (x) is the proton density, k represents R F inhomogeneity and other factors, s, is the frac58 CHAPTER 3. LINEAR COMBINATION tional contribution of component T 2 j 3.2. CALCULATION OF COEFFICIENTS 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 l ic a An effective criterion for calculating the coefficients was to maximize the S N R at a specified T2,0 : SNR/ = ^ W ^ i ^ L 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, c , for short T selection were calculated by minimizing the inverse of the S N R s 2 (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 < = 8NB.EE•••*«P(-VT.*) ( 3. 1 2 ) T 0 was defined to be 25 ms (as opposed to 20 ms) to force a broader filter in the short T range 2 2 of interest. Four constraints were necessary in the minimization procedure for calculating c : s 1. A {T ) > 0 for 10 < T < 70 ms S 2 2 Force the filter A ( T ) to be positive in the range of interest for the short T . Previous to S 2 2 adding this constraint, the algorithm tended to find a minimum where A ( T ) fluctuated S 2 positive and negative in the range 10 < T < 70 ms. 2 2. A ( T ) < 5 for T > 80 ms where 6 = max [A (T )] /100 S 2 2 5 2 This is the key constraint to suppress signal from T outside of the range of interest. Intu2 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. itively, the magnitude of the filter response for T for T 2 C A L C U L A T I O N OF COEFFICIENTS > 80 ms, relative to the filter response = 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 2 component. 3. £ c s = 0 The filter A ( T ) must go to 0 at long T . S 2 2 4. fA .dT = 1 ± 0 . 0 1 s 2 The area under the filter must be near 1 so there is no increase or decrease in I S L from the filter. 0.05 0.04 0.03 '"CM t < 0.02 0.01 5 0 A (T )>0 s 2 t -0.01 0 50 T 2 100 (ms) Figure 3.1: Constraints on the short T 150 2 200 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 > 10 ms, c , were calculated by minimizing the inverse of a 2 the S N R in Equation 3.11 based on a set of assumptions and constraints. 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: SNR where T a = SNRoE =i<exp(-t /T l l 2 i 0 ) (3.13) = 25 ms. 2 0 Two constraints were included in the minimization procedure for calculating c : a 1. A ( T ) > 0 for T > 10 ms a 2 2 Force the filter A ( T ) to be positive for T greater than 10 ms. a 2 2 2. | A ( T ) - M\ < S for T > 10 ms where M = max L4 (T )], 6 = M / 4 0 a a 2 2 2 This will constrain the fluctuations in the filter response to be less than 1/40 for T > 10 ms. 2 0.01 0.009 0.008 0.007 a ~ 0.006 < 0.005 0.004 A (T )>0 a 2 0.003 0.002 0.001 0 50 T 2 100 (ms) 150 Figure 3.2: Constraints on the "all" T 2 200 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 LevenburgMarquardt routine in Matlab (Natick, M A ) . The short T fraction was the signal calculated as 2 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. C A L C U L A T I O N OF COEFFICIENTS the linear combination of the short T coefficients with the images divided by the signal calculated 2 from linear combination of the total T coefficients with the images. 2 The constraints, incorporated into the minimization algorithm, required a T filter be created, 2 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 distribution. Small S were tried but the minimization 2 algorithm was not able to find a feasible solution. 3.2.4 Results The coefficients, c and c calculated based on the minimization procedure above, are shown i i i s a Table 3.1. The corresponding filters, A ( T ) and A ( T ) , are shown in Figure 3.3. s a 2 t (ms) 10 20 30 40 50 60 70 80 90 100 110 120 130 140 150 160 c s 3.0820 -0.5698 -1.8381 -1.9002 -1.4213 -0.7650 -0.1330 0.3858 0.7583 0.9756 1.0576 1.0266 0.9113 0.7384 0.5188 0.2860 2 c t (ms) c c 4.8727 -5.9451 -1.3997 1.9780 2.5442 1.2543 0.6837 -0.4731 -1.3316 -1.4437 -0.7389 -0.9576 -0.0542 0.1944 0.6518 0.2573 170 180 190 200 210 220 230 240 250 260 270 280 290 300 310 320 0.0576 -0.1574 -0.3481 -0.5033 -0.6275 -0.6940 -0.7184 -0.6962 -0.6253 -0.5100 -0.3503 -0.1441 0.1020 0.3769 0.6918 1.0333 -0.4878 -0.4606 1.1084 1.4641 -0.0698 -0.1675 -0.1988 0.3600 -0.1289 0.4137 0.3065 0.1169 -0.0424 -0.4811 -0.7215 -0.1254 a s a Table 3.1: The coefficients for the linear combination method. Coefficients c were for the short T filter and c were for selecting all T . s a 2 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. C A L C U L A T I O N OF COEFFICIENTS 1 0.9 0.8 0.7 CO c 0.6 1 0.4 0.3 0.2 0.1 0 T (ms) 2 Figure 3.3: T filters for short T components ( A ( T ) , solid) and all T components ( A ( T ) , dashed). The maximum of each curve was normalized to 1 for sake of comparison. S 2 3.2.5 2 a 2 2 2 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 sa 15 ms, relative to the noise. 2 The variable T 2 0 was defined slightly longer than that of myelin water to increase the width of the filter ^4 (T ). Larger values of T S 2 2 0 were tried, but resulted in poor filters A (T ) or large S 2 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 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 NNLS. 3.3.1 Model A white matter model was created based on the equation: I =p[f exp(-t /T ) s + f exp(-t /T )} s i l m 2 where p = 1000 was the proton density, T S 2 and T m used ( / s 2 = 80 ms were the fraction and T (3.14) m l 2 — 20 ms for the short component, f m = 1—f s of the medium component. Seven values for f were s 2 = 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 f was calculated by both s the linear combination and the N N L S algorithms. Other parameters included N = 32 echoes, t = 10, 2 0 , . . . , 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 f — f . The estimated short fraction f was compared to the true s s s fraction f by both a t-test (assumes the data comes from a normally distributed population) and s 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 f were normally s distributed. 3.3.2 Results and Discussion The mean short T fractions calculated by both methods agreed with the true short fraction 2 (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 with the true f s — 0.20. The mean short T 3.3. LINEAR COMBINATION SIMULATIONS 2 fraction was within 0.02 of the true fraction. The short T fractions calculated by the linear combination method were normally distributed for each 2 of the true fractions, but only when f s > 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 image was 2 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 3.4.1 2 phantoms. Methods The phantoms (Section 1.8) were scanned using the 32 echo M E S E pulse sequence (Section 1.6.2). 0 The images were combined linearly using the short T filter weights defined in Table 3.1. 2 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 map was analyzed quantitatively by calculating the mean and standard deviation 2 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 mean signal of the longer T 2 phantoms. 2 map was at least fifteen times higher than the 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 of approximately 100 ms compared 2 to the phantoms with T greater than 200 ms. This was not unexpected as the T filter is close 2 2 to zero near 100 ms, but not zero. Overall the quantitative selectivity of the short T was excellent. 66 2 phantoms 3.4. M O N O - E X P O N E N T I A L PHANTOMS 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 Figure 3.5: M R I image ( T E =10 ms, left) of nine nickel/agarose phantoms (and the large water bottle) and the short T image (right) based on the linear combination of the original 32 spin-echo images. 2 Region Mean Signal 1~ 2 3 4 5 6 7 8 9 Outside 1 Outside 2 Outside 3 Water 1 Water 2 Water 3 -0.3 (17.1) -12.3 (34.6) 470.9 (30.8) 433.5 (19.4) -9.9 (15.1) -2.7 (15.0) -28.9 (22.6) -29.1 (20.0) 547.1 (28.8) 0.7 (10.4) -1.4 (10.5) 0.8 (10.5) 6.8 (15.9) 4.3 (16.6) 8.0 (13.6) T 2 (ms) 264 101 25 23 97 363 157 77 23 2000 2000 2000 Table 3.2: The mean and standard deviation of the signal in the short T 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. 2 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 FILTER 1.2r T 2 (ms) Figure 3.6: The mean signal from each phantom plotted against the T of the phantom overlayed on the short T filter. 2 2 3.4.3 Conclusions The selectivity of the short T phantoms, from the linear combination method was excellent as 2 was the suppression of the long T phantoms. For mono-exponential phantoms, the filter behaved 2 as well as expected. 3.5 In Vivo Myelin Water Fraction: NNLS vs. T Filter 2 The selectivity and suppression of mono-exponential phantoms was excellent. In this section, the linear combination method is applied to a multi-T component system. The short T 2 calculated by the T 2 2 fraction 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.1 3.5. NNLS VS. T 2 FILTER 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 fractions from each technique were compared with a 2 Wilcoxon signed rank test. A Bland and A R m a n [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 d calculated by measurement technique 2, Bland and Altman plotted the difference of the points 2 dj — d , versus the mean of the points (dj + d )/2 for all i. 2 3.5.2 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 in the short compartment. A n increase in the short T 2 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 ( T E = 10 ms), single slice optimized pulse sequence. Bright circles were the same short T nickel/agarose phantoms as shown in Figure 3.5 2 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 FILTER 2 LinCom NNLS 2-2.5% 0.25 0.2r C 1 0.15 g LL I C a> 0.1 0.05 I 0 1 I 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.06, 1.96'SD > ' 1 • • 0.03 • 0.02 -0.02 8 c • 0) -0.04 • Mean -0.06 0 i» ! -0.08 • • 0.01 > * b - •i • -0.01 ! . - 1 96 S O -0.02 * 0.06 1.96 S D 0.04| 0 -oj i • •; o CD ! 0.05 008 0.1 0.12 0.14 M e a n of Pairs 0 16 0.18 0i -0.0.3 -0.01 -1.96 S D 0.01 0.02 0.03 M e a n 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.5.3 3.6. COEFFICIENTS B Y M A T R I X INVERSION Conclusions Myelin water fractions calculated from linear combinations of T decay data is extremely fast and 2 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 will be very broad compared to the filters defined in Section 3.2. In this section, I calculated coefficients to select for short T T (c , 70 < T m 2 to 2 < 120 ms), and long T (c , 10 < T s 2 2 < 50 ms), medium (c , T > 2000 ms). Each set of coefficients was applied ( 2 2 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 j ) ] , and L was the amplitudes which correspond to the 2 T . The coefficients were the rows of F~ l 2 where Fij — [exp ( — i j / T ) ] for t = 10, 2 0 , . . . , 320 ms 2j and T = 20, 80,3000 ms. 2 The coefficients calculated are shown in Table 3.3. The coefficients for the short T , c , follow s 2 a similar trend of positive in the early echoes, then negative, and then positive. The T filters, / ( T ) , shown in Figure 3.3, were constructed by calculating the filter / ' ( T ) = 2 2 2 (F~ ) exp ( - i / T ) . The short T filter, top left filter of Figure 3.10, should be compared with 1 t 2 2 the short T filter in Figure 3.3 (solid line). 2 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 T (ms) 2 T2 (ms; Figure 3.10: The filters calculated from the matrix inversion technique to calculate coefficients. Short T filter (top left), medium T filter (top right) and the long T filter (bottom S M left). The short T filter calculated by the inversion method {A ' (T2), dashed line) compared to the short T filter (^4 (T ), solid line) calculated by the non-linear method (bottom right, plotted on a log scale to visualize the whole filter). 2 2 2 2 S 2 2 73 C H A P T E R 3. L I N E A R COMBINATION t (ms) c 10 20 30 40 50 60 70 80 90 100 110 120 130 140 150 160 1.72393 0.49791 -0.16482 -0.49547 -0.63314 -0.66120 -0.62934 -0.56695 -0.49114 -0.41172 -0.33408 -0.26105 -0.19391 -0.13308 -0.07849 -0.02982 s -0.40229 0.08544 0.32873 0.42999 0.45058 0.42708 0.38113 0.32531 0.26685 0.20972 0.15599 0.10661 0.06188 0.02176 -0.01400 -0.04573 3.6. COEFFICIENTS B Y MATRIX INVERSION c' t (ms) c 0.04356 -0.01328 -0.03897 -0.04681 -0.04474 -0.03749 -0.02781 -0.01728 -0.00681 0.00314 0.01233 0.02068 0.02819 0.03487 0.04080 0.04603 170 180 190 200 210 220 230 240 250 260 270 280 290 300 310 320 0.01337 0.05159 0.08533 0.11505 0.14120 0.16418 0.18434 0.20202 0.21749 0.23102 0.24283 0.25313 0.26209 0.26987 0.27660 0.28242 Table 3.3: The coefficients for the short, medium and long T inversion method of calculation. 3.6.2 c' s -0.07379 -0.09854 -0.12034 -0.13950 -0.15632 -0.17107 -0.18398 -0.19527 -0.20513 -0.21372 -0.22120 -0.22768 -0.23330 -0.23815 -0.24233 -0.24591 0.05064 0.05469 0.05823 0.06134 0.06404 0.06640 0.06845 0.07022 0.07176 0.07308 0.07421 0.07518 0.07601 0.07670 0.07728 0.07776 selection based on the matrix 2 Results and Discussion The short T filter (Figure 3.10, top left) is a very broad peak centered on T 2 crosses 0 at T = 80 ms. The primary problem with the short T 2 2 2 = 22 ms and 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 > 80 ms. If the system has only mono-exponential contributions for all 2 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 was not lower. The short T filter, calculated 2 from this method, is compared to the short T 2 filter, A ( T ) calculated in Section 3.2.4 and is S 2 2 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 > 80 ms will be negative. 2 The medium T 2 filter (Figure 3.10, top right), and long T have a very broad selection. 74 2 filter (Figure 3.10, bottom), both 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 T 2 3.6. COEFFICIENTS B Y MATRIX INVERSION S e l e c t e d 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 (top left), medium T (top right), long T (bottom left) and the myelin water fraction (ratio of the short T to the sum of the three, bottom right). 2 2 2 2 75 3.7. BACKUS-GILBERT COEFFICIENTS 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.3 Conclusions The method of calculating the coefficients by matrix inversion is extremely fast, but lacks the ability to apply constraints on the coefficients or filter. The filter would be appropriate for use to select T when the underlying data are mono-exponential, but would not be appropriate when 2 in vivo brain. the underlying data are multi-exponential as is the case for 3.7 Calculation of Coefficients by Backus-Gilbert The Backus-Gilbert method is a linear algorithm to calculate the coefficients, c , which minimize 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 A (T ) as defined in Section 3.2.4. S corresponding filters are compared to the filter 3.7.1 2 C a l c u l a t i o n of Coefficients The standard misfit function, to minimize, includes two terms: the first term is an error term between the constructed filter ("averaging function") and a target filter and the second term is a noise amplification term. The two terms are linearly, combined with a trade-off parameter, 6, which can be adjusted to arbitrarily weight one term more than another. The coefficients, c, were calculated to.minimize: TV /•oo <t> = cos(0) [A(T ,T 2 2 ) 0 )-n(T ,T , 2 where 6 is the trade-off parameter, A ( T , T 2 2 > 0 2 ) = ,T m m ) ] d T + m(0) 2 2 ! m a x 2 ^ exp(-TE/T 2 ] 0 S Vc 2 cr ) , and n(T , T 2 2 2 ) m i n . (3.15) , T , 2 m a x ) is the boxcar function defined as: TTfT J-2,mm> T \—) 1 ^^< 2 min - i-2,max) ~ \ I 0 otherwise T 2 - 2,max T • (o T C \ (3.10J 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 76 (3.17) 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, c , had a higher trade-off parameter, the second set, bgl c , had a lower trade-off parameter and third set, c bg2 , was designed to have a similar noise ba3 amplification factor as c . The short T filters, corresponding to each set of coefficients, is shown s 2 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 of interest. The differences in the filter were 2 around T = 500 ms. The filters, corresponding to the coefficients c b91 2 and c , deviated from ba3 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 > 100 ms is undesirable as the filter was designed for short T . The filter, corresponding 2 2 to coefficients, c , had a similar response to the filter defined by the coefficients calculated in b92 was 113.4 for, c bg2 his paper, but the noise amplification t (ms) 10 20 30 40 50 60 70 80 90 100 110 120 130 140 150 160 bgi c 3.3504 -1.6519 -1.6415 -1.1697 -0.7032 -0.3204 -0.0287 0.1796 0.3183 0.3985 0.4340 0.4337 0.4088 0.3642 0.3078 0.2455 bg2 c 6.8123 -7.5974 -2.2797 -0.1489 0.5016 0.6398 0.6525 0.6430 0.5675 0.5403 0.4623 0.3956 0.2962 0.2009 0.1028 -0.0174 bg-i t (ms) c 4.0186 -2.3374 -2.0012 -1.2865 -0.6791 -0.2204 0.1119 0.3350 0.4668 0.5312 0.5414 0.5123 0.4564 0.3812 0.2950 0.2066 3.0820 -0.5698 -1.8381 -1.9002 -1.4213 -0.7650 -0.1330 0.3858 0.7583 0.9756 1.0576 1.0266 0.9113 0.7384 0.5188 0.2860 170 180 190 200 210 220 230 240 250 260 270 280 290 300 310 320 bgl c 0.1805 0.1162 0.0558 0.0006 -0.0470 -0.0873 -0.1196 -0.1420 -0.1561 -0.1611 -0.1575 -0.1456 -0.1259 -0.0980 -0.0634 -0.0228 and was 29.8 for bg2 c -0.1123 -0.2007 -0.2456 -0.3335 -0.3556 -0.3552 -0.3735 -0.3452 -0.2870 -0.2300 -0.1456 -0.0473 0.0793 0.2389 0.3758 0.5630 bgi c 0.1191 0.0372 -0.0350 -0.0984 -0.1489 -0.1872 -0^2113 -0.2212 -0.2187 -0.2021 -0.1741 -0.1314 -0.0800 -0.0171 0.0556 0.1385 0.0576 -0.1574 -0.3481 -0.5033 -0:6275 -0.6940 -0.7184 -0.6962 -0.6253 -0.5100 -0.3503 -0.1441 0.1020 0.3769 0.6918 1.0333 Table 3.4: Coefficients calculated by the Backus-Gilbert method with with a different set of tradeoff parameters and, therefore, properties. Coefficients calculated by the method defined in this thesis, c , are shown for reference. s 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 1I 1 3.7. BACKUS-GILBERT COEFFICIENTS 1i 1 c c - - c .... c bg1 bg2 b 9 3 ' s • , , , , «* \ lt >^y 1\ - 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 , calculated in Section 3.2 are shown for reference. s 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 ( T 2 ) , calculated based on the coefficients c , is shown for reference. S s 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.7.3 3.8. OPTIMAL E C H O TIMES 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 ( T ) , resulted in a smaller contribution for T > 100 ms and a small noise amplification. S 2 2 Better filters can be created when non-linear constraints can be used. Optimal Echo Times 3.8 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 as a function of 2 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 filters, were defined to allow comparison as the number 2 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 ) = E i l i °i P ( * i / T ) , where N is the number of e x _ 2 2 echoes. The goal of the filter is to be as wide as possible for 10 < T T 2 > 80 ms. Figure 3.14 shows a short T 2 2 < 50 ms and near zero for 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 F W 1 0 0 M is important as the medium T peak around 2 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 3.8. OPTIMAL ECHO TIMES 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 0.045 r -0.005" 0 1 20 1 40 1 1 60 T (ms) 1 80 100 : ' 120 2 Figure 3.14: The measures of the effectiveness of a short T filter. The full width at half the maximum should be as large as possible and the T where the filter reaches 1/100 of the maximum should be less than T = 80 ms and greater than T = 50 ms. 2 2 2 81 2 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.2 3.8. OPTIMAL ECHO TIMES E q u a l l y 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 filter is 2 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 coefficients. The objective function used was the same as the one defined previously. The inverse of the objective function was defined as: max 2~2j=i Cjexp (-tj/2b) (3.18) c,A4 with constraints - 2 < Cj < 2 for i = 1... N, 10 < At < 1000 ms, and t = At, 2 At,NAt. The constraints on Cj were defined to constrain the noise amplification ^ 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 ( T ) . The width is relatively S 2 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. 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 Number of Echoes 100 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 Number of Echoes 25 30 35 Figure 3.16: The full width at half max ( F W H M ) of the filter A {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. S 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 filter had little change for 5 < N < 19 and 2 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: max Zi=i^ e x p(-^/ 2 5 ) (3.19) t,c 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 defined to constrain the noise amplification 3.8.3.2 Results The echo times were calculated along with the coefficients. The minimum spacing between consecutive 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. 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. O P T I M A L 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 F W 1 0 0 M was approximately 2 ms less, when the echo sampling was allowed to vary compared to echoes equally spaced. 400 40 0 Number of Echoes 100 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 ( F W H M ) of the filter A {T2) as a function of the number of echoes (top) and the full width at 100 max (bottom). S th 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 space. 2 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 filter, shorter echo times would be required, which would be difficult 2 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 filter, A ( T ) , for 10 < T 5 2 2 2 < 50 ms is a function of the echo and corresponding coefficients. Figure 3.19 shows the T information in spin-echo signal acquired at 2 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 < t. Therefore, a sharp rise on the left side of A ( T ) requires signal S 2 2 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 range of interest. 2 The signal from each echo contributes different information depending on the echo time of the acquisition. The 32 coefficients, c , in Table 3.1 vary from positive to negative and back to s 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 ( T ) is shown S 2 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 * /" / / / / / / / TE=10 TE=80 - - TE=120 TE=200 TE=800 s • / / 50 100 150 200 250 T (ms) 300 350 400 2 Figure 3.19: The filter, as a function of T , 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 —> oo. 2 2 0.15 1 I After 1 echo " / / N \ ~. ~~ ~ ~~~ ^^^^.c-. \ After 32 .echoes After 28 echoes ; '"V N. 'V s. -0.05 -0.111 0 I. I. 50 100 =l=l 150 T (ms) 2 Figure 3.20: The short T 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. 2 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.5 3.8. OPTIMAL ECHO TIMES 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) c 10 30 120 270 3.1872 -4.7952 2.3515 -0.7539 Table 3.5: Coefficients for short and all T 2 c s a 2.5642 -2.5133 1.8698 -1.1654 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 filter. The first echo time, T E = 10 ms, defines 2 the rapid increase in the portion of the curve for 0 < T < 18 ms. The second and third echoes, 2 at T E = 30 ms and T E = 120 ms, bracket a portion of the T filter where the filter changes from 2 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 culated from 32 echoes. The ascending portion of the filter (0 ;$ T 2 2 filter cal- ;$ 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 T 2 2 filter, from 4 echoes, was not as close to zero as the short T 2 filter from 32 echoes for > 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 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 2 several echoes at small T E , one echo near the T 2 at long T E . 91 where the filter dropped to zero and one echo 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 ~ 20 ms. 2 3.9 Conclusions Fast accurate estimation of myelin water is important for assessing diseases such as multiple sclerosis. The current method to quantify the myelin water fraction in vivo is to acquire multiecho images and compute the distribution of T voxel-by-voxel. The myelin water fraction is 2 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 distribution for each voxel in an image is a slow process. 2 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 < 50 ms (short T ) and a second set to select water with T > 10 ms 2 2 2 (all T ) . The coefficients to select for short T 2 phantoms. The resulting short T 2 2 were applied to multi-echo images of a set of image had very good selection of the phantoms with T 20 ms and good suppression of phantoms with T 2 2 near > 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 the 3.9. CONCLUSIONS 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 Non-180 Refocusing 4.1 Pulses Introduction 4.1.1 Motivation Accurate T estimation traditionally depends on accurate refocusing pulses. Many studies [31, 35, 2 71, 72] showed large deviations in the measured T when the refocusing pulse flip angle deviates 2 from 180°. For example, Majumdar [71] measured biases in R 2 (R 2 = 1/T ) of greater than 2 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 i th focusing pulses, to be attenuated by the i th echo, measured from a train of N re- 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 separate acquisition was required to calculate the B i map. The calculated T 2 2 and a was divided by f (a) to obtain.a "corrected" mono-exponential T - Another approach by Lepage [76] corrected 2 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 the experimentally determined T 2 4.1. INTRODUCTION 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 M (t z + t) = F(t) + t) = M (t) w w z cos 2 (a/2)+F*(t) cos (a) sin - ^i[F(t) 2 (a/2)-iM (t) - F*(t)} sin (a) z (4.1) sin (a) where F(t) is the dephasing transverse magnetization, F*(t) is the rephasing transverse magnetization, M (t) z is the longitudinal magnetization, t w 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 coherence through a multi-pulse experiment. Magnetization after a non-180 refocusing pulse is dis0 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 accumulation 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 relaxation of e x p ( - 2 r / T ) for F(t) 2 and F*(t), e x p ( - 2 r / T i ) for M (t), z MQ (1 — exp(—2r/Ti)) (where 2r is the echo spacing). For accurate T 95 and a generating term 2 quantification, Hennig 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. O n 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 M ,M if the refocusing pulses were 180°. x 96 y 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.2.1 4.1. INTRODUCTION Mathematical Description of Woessner's Paper Woessner [37] defines the Bloch equations i n terms of U (in phase component), V (out-of-phase component) and M (longitudinal component) as: z jU(t) = jV{t) = ±M.(t) = -( - )V-U(t)/T Uo u 2 {Lj -u3)U-V(t)lT -uiM (t) Q _ ^ L Z p ) 2 + W l V z (4.2) W 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)\t (where t w w is the width of the R F pulse): U{t + t ) = U(t) V(t + t ) = V(t) cos (0) - M sin (0) M (t + t ) = V(<)sin(0) + M cos(6») w w z w z 2 (4.3) and F{t + t ) w = U (t + t ) + iV (t + t ) w w = U (i) + iV (t) cos (9) - iM sin (9) z 97 (4.4) 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 Equation 4.4 is then written as a function of F(t + t ) 4.1. INTRODUCTION F(t), F*(t), M(t), and M*(t): = U(t) +iV(t) cos {9) - iM (t)sin(0) w (4.5) z = U(t) [cos (9/2) + sin (9/2)] + iV(t) [cos (9/2) 2 2 2 sin (9/2)] - iM (t) sin (9) (4.6) 2 z = U(t) cos (9/2) +U(t) sm (9/2)+ iV(t) cos (9/2)2 2 2 iV(t) sin (9/2) - iM (t) sin (9) (4.7) 2 z = U(t) cos (0/2) + iV(t) cos (0/2) + r/(i) sin (9/2) 2 2 iV(t) sin (9/2)-iM (t) sin (9) 2 z - 2 • (4.8) = [U(t) + iV(t)} cos (9/2) + [C/(t) - iV(t)] sin (0/2) - i M (t) sin (0) = F(i)cos (0/2) + F*(t)sin (0/2)-iM (i)sin(0) 2 2 2 2 (4.9) (4.10) 2 z by using the cosine double angle formula: cos (0) = cos (0/2) - sin (0/2). (4.11) cos (0/2) + sin (0/2) = 1. (4.12) 2 2 and the identity: 2 Then, F(t + t ) and M (t + t ), written as functions of F(t) and M (t), will be: w z F(t + t ) w w = z F(i)cos (0/2) + F*(i)sin (0/2)-iM (i)sin(0) 2 2 2 F*(t + t ) = F(t)sin (9/2) + F* (t) cos (912)+iM* (t) sin (9) M (t + t ) = M (t)cos(9)-^[F(t)-F*(t)}sin(9) M*(t + t ) = M*(t) cos (9) - \[F(t) - F*(t)]sin(9) w z w w 4.1.3 2 2 2 z z (4.13) (4.14) (4.15) (4.16) I m p l e m e n t a t i o n 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 Algorithm 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. 98 Then, the algo- 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 substate 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 defocusing 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 e x p ( — r / T i ) . The important aspect, for this work, is e x p ( — r / T i ) « 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 , a, and p. 2 The signal intensity, y, was calculated (Algorithm 4.1) given T , T i , a, r, and the number of 2 refocusing pulses. The algorithm calculated the 32-echo decay curve in approximately 25 ms on a P4 1.5 G H z 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 function [decay] % alpha % num.echoes % Tl = s t i m u l a t e ( a l p h a , num_echoes, = rotation = = Tl T2 of tissue (ms) T2 = tau = time % Preallocate between % the signal Preallocate tau) (ms) intensity of num^echoes rotation matrix Tone (cos(alpha/2))"2, Calculate substates [F F* Apply Z Z*] and initialize 1); jay*l; Precalculate [[ vector for one substate block (sin(alpha/2))"2, -jay*sin(alpha), 0]; (sin(alpha/2))~2, [ (-jay)*l/2*sin(alpha), (jay)*l/2*sin(alpha), cos(alpha), [ (-jay)*l/2*sin(alpha), (jay)*l/2*sin(alpha), 0, the sparse of "num-echoes" rotation (cos(alpha/2))"2, matrix initial B = t2relax(B, T2 decay tau, and Tl regrowth 0, . . . [ 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 ) ) ) , % T l , T2, reference. num_echoes) ; the array = 0 + Tone = % of echoes. 90-180 B = zeros(num_echoes*4, B(l) frame sqrt(-l); decay = z e r o s G , % of (ms) % = number in x'y'z of tissue % jay vector jay*sin(alpha)]; . . . 0] ; cos(alpha)]]; blocks create.Tp(alpha)); (between 90-180 pulses). T2); B = tlrelax(B, tau, T l ) ; for echo = 1:nunuechoes % Apply % Each rotation block matrix. of the rotation matrix acts on [F(k), B = (Tp * B ) ; % Save echo decay(echo) % Apply intensity, transverse B = t2relax(B, % Apply this = abs(B(2)) decay 2*tau, longitudinal Let the spins dephase Fl* for time 2*tau. T2); regrowth B = t l r e l a x ( B , 2*tau, % is * exp(-tau/T2); for time 2*tau. Tl); for time 2*tau B = evolution(B); end 100 F*(k), Z(k), Z*(k)] . . . 4.1. 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 Algorithm 4.2 T i function % [B] The Z and Z* indicies % = Apply substates the Tl % % [B] relaxation be = 2 [ modified 4:4:length(B) * ]; substates exp(-tau/Tl); relaxation substates tau, will 1:4:length(B) the T2 relaxation B(indicies) Tl) to the Z and Z* = t2relax(B, The F and F* Apply will = B(indicies) Algorithm 4.3 T indicies tau, [ 3:4:length(B) B(indicies) function relaxation = tlrelax(B, be T2) modified 2:4:length(B) to the Z and Z* = B(indicies) * ]; substates exp(-tau/T2); Algorithm 4.4 Phase evolution of magnetization function % [Bnew] This % as they = evolution(B) function is a book-keeping are flipped and as they method progress to track through Bnew = B ; % L % Number = of elements All F(n) go to Bnew([ 5 : 4 : L ]) % F*(l) Bnew(l) % in B length(B); All goes = F(n+1). = B([ to l:4:L-4 ]); F(l). B(2); F*(n) go to Bnew([ 2 : 4 : L - 4 ]) F*(n-1). = B([ 6:4:L ]); 101 the phases the inter-pulse of the time. spins INTRODUCTION C H A P T E R 4. NON-180 R E F O C U S I N G 4.1.4 PULSES 4.2. MONO-EXPONENTIAL FIT Overview of Work I created a method tofitdecay 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 sequence C 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 sequence were C compared to myelin water maps calculated from the M E S E . Several extensions to the work will 0 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 tofitmono-exponential decay curves is described. Noise simulations were used to assess the accuracy and precision of thefitalgorithm. Finally, a set of nickel/agarose phantoms, of known T values, were scanned with different sets of refocusing pulses. The estimated T 2 2 was compared to the T estimated from the M E S E sequence and the refocusing pulse were compared 2 0 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 , baseline 2 offset d, and refocusing pulse flip angle a. A non- linear least squares algorithm ( l s q c u r v e f i t in Matlab [The Mathworks, Natick, MA]) minimized the misfit x = YliLiiVi 2 ~ Vi) l 1 2 a 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 = 100 ms, 2 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.2 4.2. MONO-EXPONENTIAL FIT 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-180 refocusing pulses. The mean and standard 0 deviation were calculated for each parameter over the noise realizations. Methods A noiseless decay curve was created based on the model: y = pexp(-t/T ) (4.17) 2 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, 2 0 , . . . , 320 ms. One thousand realizations of quadrature noise (SNR =200) were added. T R was set to infinity. The T decay curve parameters 2 used in the simulation are similar to values for brain tissue. The mean and standard deviation of the estimated parameters were calculated and compared to the true values. Results and Discussion Each parameter estimated from the decay curve, for T 2 = 20,80,140 ms was within 5% of the true parameter. Table 4.1 shows the results for T = 80 ms. The parameters estimated from the 2 decay curve were within 2% of the true parameter for T = 80,140 ms. The increased accuracy, as 2 a function of T , was expected as more echoes were able to contribute information to the solution 2 as T increased. 2 The variability in the estimated parameters was higher for the lower flip angles (for example, variability in T in Table 4.1). To determine if the increased variability is related to the decrease 2 in the flip angle, I re-ran the simulations and scaled the 1,000 realizations of noise by, sin (a/2), 2 on the relative attenuation due to the flip angle. For example, I multiplied the random noise by sin (a/2) = 1.0 for a = 180° and by sin (a/2) = 0.5 for a = 90°. In this way the noise decay 2 2 curves estimated with a = 90° had the same S N R as the curves simulated with a = 180°. Then any discrepancies in the standard deviations of the parameters, were not be due to S N R issues. 103 C H A P T E R 4. NON-180 R E F O C U S I N G PULSES "true ( ° ) A(°) P 090 100 110 120 130 140 150 160 170 180 89.836 (0.964) 99.848 (1.015) 109.849 (1.117) 119.851 (1.148) 129.839 (1.269) 139.848 (1.450) 149.841 (1.669) 159.851 (1.984) 169.992 (2.726) 177.474 (2.925) 1001.220 (12.455) 1000.919 (11.116) 1000.774 (10.058) 1000.622 (8.817) 1000.551 (8.040) 1000.370 (7.292) 1000.269 (6.490) 1000.144 (5.717) 1000.027 (4.823) 1000.391 (4.179) Table 4.1: The estimated p and T true T = 80 ms. 4.2. MONO-EXPONENTIAL FIT T ( ms) 79 296 (1.266) 79 407 (1.139) 79 492 (1.035) 79 557 (0.939) 79 594 (0.867) 79 630 (0.799) 79 650 (0.735) 79 669 (0.686) 79 681 (0.638) 79 636 (0.605) x 1 2 28 28 28 28 28 28 28 28 28 28 413 391 285 333 270 185 137 141 177 642 (8 050) (8 017) (8 025) (8 070) (8 035) (8 058) (7 992) (8 Oil) (8 049) (8 030) (standard deviation in brackets) given the true p = 1000 and 2 2 "true ( ° ) 090 100 110 120 130 140 150 . 160 170 180 a{°) 89.900 (0.464) 99.902 (0.587) 109.912 (0.670) 119.932 (0.887) 129.982 (1.026) 139.988 (1.245) 149.964 (1.592) 160.159 (1.969) 170.412 (2.770) 177.861 (2.579) P 1000.469 (5.860) 1000.287 (6.211) 1000.094 (6.026) 999.904 (6.408) 999.612 (6.367) 999.597 (6.037) 999.713 (5.996) 999.381 (5.327) 999.490 (4.645) 1000.032 (4.348) T ( ms) 532 (0.658) 588 (0.709) 638 (0.711) 678 (0.755) 707 (0.756) 704 (0.728) 689 (0.719) 711 (0.696) 693 (0.672) 634 (0.661) 2 79 79 79 79 79 79 79 79 79 79 X 1 7.132 (1.753) 9.794 (2.384) 12.822 (3.137) 15.928 (3.916) 19.117 (4.719) 22.015 (5.450) 24.446 (6.032) 26.423 (6.618) 27.776 (6.919) 28.757 (6.986) Table 4.2: The estimated p and T from a simulation with S N R scaled as a factor of sin (a/2) to determine whether the increased variability in the estimated parameters was due to the decrease in a. 2 2 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 i n Table 4.2. The standard deviation of p and T remain relatively constant as a decreased to 90°. 2 The constant variability was true for T = 140 ms as well. When T = 20 ms, the variability was 2 2 still a function of a, though decreased by approximately half when the S N R was scaled by the effect of the refocusing pulse. Conclusions Mono-exponential T parameters were accurately and precisely estimated from decay curves sim2 ulated from non-180 refocusing pulses. The accuracy and precision were found to depend on a, 0 and the short T was affected more by the decreased flip angle than the longer T components. 2 4.2.3 2 Phantom Verification A set of mono-exponential phantoms were used to confirm the accuracy of T from a set of scans. 2 A set of six nickel/agarose phantoms of known T were scanned (Section 1.8) with a 16 echo pulse 2 sequence with prescribed refocusing pulse flip angles. Methods The data was acquired with the single slice, 16 echo M E S E pulse sequence with parameters: echo C 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 a , flip angle set to aiso = 180°, and flip angle estimated to minimize p 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 FIT 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 a = 180°. The T was typically within ± 1 0 % of the T estimated from the 48 echo p MESE 0 2 2 data. OLp 90 110 130 150 170 180 "best 85 103 123 137 151 156 P180 PP 498.2 596.8 687.3 750.3 784.3 800.2 792.7 777.2 783.1 786.2 788.5 800.2 Pbest 842.7 816.5 813.0 826.2 819.2 819.9 T ,180 2 136.4 123.7 110.4 101.7 98.3 97.0 t 2,p 88.5 96.4 97.4 97.2 97.7 97.0 2,best , 83.4 87.8 92.4 92.5 94.0 92.0 t Xl80 4 .141.1 81.8 49.9 23.7 7.5 5.4 4.0 3.3 3.4 6.5 5.7 5.4 Xbest 1.9 0.6 0.9 0.8 1.3 1.0 Table 4.3: Phantom decay curve parameters estimated from mono-exponential fits to decay curves collected using a train of 16 a refocusing pulses. Parameters with a subscript of p were estimated with a fixed to a . Parameters with a subscript of best were estimated allowing a to vary. Results are from one phantom ( T = 90.8 ms). F l i p angles are in degrees and T in milliseconds. p p 2 2 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 and a. To assess the effect of 2 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 , p and a from each phantom decay curve were estimated to be the same as the result from 2 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 phantoms (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 a — 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 = 90.8 ms for the phantom. p 2 107 C H A P T E R 4. NON-180 R E F O C U S I N G PULSES 4.2. MONO-EXPONENTIAL FIT 700 TE (ms) Figure 4.3: Decay curves calculated as the mean signal intensity over a region of one of the phantoms (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 a = 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 = 90.8 ms for the phantom. v 2 108 C H A P T E R 4. NON-180 R E F O C U S I N G PULSES 4.2. MONO-EXPONENTIAL FIT Figure 4.4: Decay curves calculated as the mean signal intensity over a region of one of the phantoms (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 a = 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 = 90.8 ms for the phantom. p 2 109 C H A P T E R 4. NON-180 R E F O C U S I N G PULSES 200 4.2. MONO-EXPONENTIAL FIT 1 1 1 n 3 TJ W O tr -200 0 50 100 150 200 Figure 4.5: Decay curves calculated as the mean signal intensity over a region of one ofthe phantoms (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 a — 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 = 90.8 ms for the phantom. p 2 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 Conclusions The decay curves from the mono-exponential phantoms were reliably fit using a mono-exponential curvefitalgorithm designed to calculate the refocusing pulse flip angle during thefit.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 x misfit. 2 Multi-Exponential T 4.3 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 Multi-Exponential Fit Algorithm 4.3.1.1 T 2 Decay Curve F i t Using 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 ^4s = y where Aij = exp ( - T E i / T ) , y is the measured 2 j decay curve and Sj are the unknown amplitudes that correspond to T 2 j . The solution minimizes the x misfit. In the standard N N L S algorithm, each column Aj, of the matrix A, is a mono2 exponential decay curve with T 2 = T . To account for a, I created a matrix A where each a 2 j column was a decay curve, calculated using Hennig's method (Section 4.1.2), with T = T j and 2 flip angle a. Therefore, the objective was to find the minimum x i ) 2 a 2 misfit of ^4 s = y over the Q range of a. To do this, I created a new algorithm, N N L S , to find the minimum x (a) by solving a 2 A s = y for a between 40° and 180°. a The multi-exponentialfitalgorithm I introduced, N N L S " , minimized x 2 a s Figure 4.6 shows three example curves of x 2 a function of a. as a function of a for a = 90° (solid line), a = 150° p p (dashed line) and a = 180° (dash-dotted line). Each of the curves has a distinct minimum at p the true a. Therefore, the goal is to find a fast, robust method to find the refocusing pulse flip angle that minimizes the x misfit. 2 The x misfit increases quickly as the a° refocusing pulse flip angle deviates from the true 2 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 200 180 11 // ! — ' / 1 160 140 : co i 100 O 80 : 1 1 1 1 ; / : / : I : \ : : \ ' True a = 90 - - • Truea=150 ' - - Truea=180 : 1 i 1 i / 1 i / . I / \I:: r * / 1 £ 120 ro 1 / 1 4.3. MULTI-EXPONENTIAL FIT \ ' > 7 \ i •: 40 \ 20 60 80 100 ' 120 Flip Angle / ' 11 140 \ *y / ^ 1 / \ * ' / V \ • / \ *: ' I / 60 1 / • ;: // t \ \ 1 i / 1 160 180 Figure 4.6: The x 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 x over a. 2 2 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 x • Two types of algorithms 2 were assessed for speed and robustness. The first type of algorithm samples several points along the x / 2 a curve then fits the points and determines the minimum x 2 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 x 2 is found by minimizing the function X {a) over a. The second technique was inherently slower, but more robust. 2 4.3.1.2 Analysis ofthe A Matrix a 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 measurements is higher when a < 180° than 2 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 A a 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 A a 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• A decay with T E = 10 ms and T R = 800 ms. a ratio (smallest to largest singular values) was 105.2686 x 1 0 ~ 0.4996 x 1 0 ~ 17 17 based on a 32 echo at a = 158° and the lowest was 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 0.0021 x 1 0 ~ 17 - 1 7 at a = 150° and the lowest was 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 A a the largest inverse condition number was A decay curves in A 2 0 8 and A 1 5 2 2 0 8 matrices. For example, the matrix with , but the difference between The mono-exponential 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. B y 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 Flip Angle (degrees) 300 Flip Angle (degrees) Figure 4.8: The rank (left) and inverse of the condition number (right) of A decay with T E = 5 ms and T R = 800 ms. a based on a 200 echo number (two orders of magnitude) should allow for better T quantification. Another benefit of a 2 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) = 1.28. Even though there is a large 2 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 S N R 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 compartment of p = 200 and T s 2 and p m = 800 and T m 2 S 2 = 20 ms = 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 S N R — 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 s had less than 10% error for a > 170° 2 (Table 4.4). A l l parameters were incorrectly determined when a < 170°. Therefore, accurate and precise quantification of T parameters requires the flip angle be calculated as well. 2 90 100 110 120 130 140 150 160 170 180 s 0.000 0.000 0.000 0.000 0.003 0.150 0.255 0.219 0.209 0.204 (0 (0 (0 (0 (0 (0 (0 (0 (0 (0 T 2m m f Q'true f 000) 000) 000) 000) 029) 147) 076) 048) 048) 044) 0.000 (0.000) 0.000 (0.000) 0.000 (0.000) 0.000 (0.000) 0.716 (5.805) 25.070 (22.278) 34.956 (7.518) 26.921 (5.214) 21.753 (5.071) 19.804 (4.620) 1.000 1.000 1.000 0.984 0.919 0.838 0.744 0.779 0.790 0.794 (0 (0 (0 (0 (0 (0 (0 (0 (0 (0 000) 000) 002) 029) 083) 142) 074) 047) 047) 044) x 1 107.707 (0.956) 98.013 (1.204) 89.886 (1.212) 82.727 (1.783) 75.387 (4.323) 82.402 (9.003) 86.947 (4.419) 82.931 (3.244) 80.270 (3.217) 79.259 (2.962) 1123.617 955.397 767.491 579.248 405.741 256.958 . 142.229 66.923 32.087 26.948 Table 4.4: Measured parameters of a bi-exponential T decay curve when the nominal flip angle is assumed to be a = 180°, but isn't. The true parameters are: f = 0.20, dn = 200, m i y = 20, f = 0.80, dn = 800, T = 80. 2 s m s m 2 4.3.3 NNLS Q M u l t i - E x p o n e n t i a l Simulations A set of noise simulations were used to assess the accuracy and precision of the parameters of the T distribution and refocusing pulse flip angle. A bi-exponential T distribution as in Section 4.3.2 2 2 was used. The true parameters were f s = 0.20, T S 2 = 20 ms, f m = 0.80, and T m 2 = 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 (a/2). When the noise was scaled, 2 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 90 100 110 120 130 140 150 160 170 180 89.87 (1.44) 99.81 (1.47) 109.83 (1.57) 119.81 (1,63) 129.80 (1.68) 139.72 (1.72) 149.88 (1.96) 159.80 (2T0) 170.01 (2.98) 176.76 (2.97) T , 15.90 (10.63) 16.65 (9.67) 17.75 (9.17) 18.06 (8.20) 18.32 (7.70) 18.33 (6.43) 18.61 (6.11) 18.29 (6.29) 17.81 (6.87) 17.97 (6.62) fs 0 20 0 20 0 21 0 20 0 21 0 20 0 20 0 20 0 20 0 21 fm 2 ) (0 13) (0 12) (0 11) (0 10) (0 09) (0 08) (0 07) (0 07) (0 07) (0 07) 0 79 0 79 0 78 0 79 0 78 0 79 0 79 0 79 0 79 0 79 T (0 17) (0 15) (0 12) (0 10) (0 10) (0 08) (0 07) (0 07) (0 07) (0 07) 2,m 73.33 (15.66) 74.50 (13.89) 76.12 (12.15) 76.55 (9.98) 76.77 (9.78) 77.06 (7.98) 77.12 (7.44) 77.16 (7.03) 77.04 (7.21) 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 / = 0.2, T2 = 20 ms, f = 0.8 and T2, 80 ms. F l i p angles are in degrees and T2 in milliseconds. = s 4.3.4 ]S m m 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 sequence, as well, and C 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 pulse sequence with T E = 10, 2 0 , . . . , 320 ms. 0 Other parameters were T R = 3000 ms, 2 averages. The T measurements were done with a 32 echo 2 M E S E pulse sequence with and echo spacing of 12.1 ms. Other parameters were T R = 3000 ms, C 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 PULSES 4.3. MULTI-EXPONENTIAL FIT first image was a standard spin-echo image, the second had an excitation pulse of 180°. Other parameters included 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]: (4.18) where I a 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. For each phantom, the mean signal intensity over the phantom was fit to the equation SI(TRi) 4.3.4.2 T 2 = p[l -exp(-TR/Ti)]. Results Calculation The T2 relaxation times, shown in Table 4.6, were calculated from the M E S E reference) and from the M E S E curve collected by the M E S E C C sequence. 0 sequence (for The T2 relaxation times calculated from the decay 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. 118 C H A P T E R 4. NON-180 REFOCUSING PULSES Phantom bottomJeft bottom-middle bottom_right top-left top-middle topjight 4.3. MULTI-EXPONENTIAL FIT M E S E T (ms) 24.3 86.4 26.1 316.6 99.3 110.6 Q 2 M E S E T (ms) 22.9 80.7 24.2 304.6 91.7 102.8 C 2 Table 4.6: T time calculated from the M E S E sequence and the MESE (with the flip angle calculated). 2 0 C The refocusing pulse flip angle was calculated by the Stollberger method (as) and from the decay curves acquired with the M E S E sequence (a) and are shown in Table 4.7. The refocusing C pulses calculated from the decay curve were within 6% of the refocusing pulses calculated by the Stollberger method. Phantom bottom-left bottom_middle bottom_right topJeft top-middle top_right a (°) 163.0 159.0 163.0 156.0 158.0 160.0 a (°) 165.7 166.3 167.7 164.8 165.1 166.9 s Table 4.7: The refocusing pulse flip angle was calculated by the Stollberger method (as) and from the decay curves acquired with the M E S E sequence (a). C Phantom bottom-left bottom_middle bottom-right top-left top-middle top-right T i (ms) 343.4 352.9 690.7 730.7 712.0 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 4.3. MULTI-EXPONENTIAL FIT 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.4.3 Conclusions The nickel/agarose phantoms were scanned using the M E S E C sequence from which the T and 2 refocusing pulse flip angle a were estimated. The T estimated from each phantom was within 2 10% of the estimated T from the M E S E o sequence. The estimated refocusing pulse flip angle 2 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 flip angle was set to a n o m n o m = 180° and for the second scan, the = 150°. From these two scans, the myelin water map was calculated, as was a map of the refocusing pulse flip angle and the mean T . 2 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 = 180°. The second scan was acquired with a n o m = 150°, by scaling the n o m 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 scan was not acquired. 0 The T distribution was estimated voxel-by-voxel in each scan as well as the refocusing pulse 2 flip angle using the algorithm described i n Section 4.3.1. Maps of the myelin water fraction, refocusing pulse flip angle and mean T calculated from each scan were compared visually. For 2 the first analysis, the T 2 information was calculated along with the flip angle at every voxel ("Variable Flip Angle") i n the image and for the second analysis the flip angle was assumed a = 180° at every voxel ("Fixed Flip Angle"). 4.3.5.2 Results Variable Flip Angle Maps of the myelin water were estimated from the multi-echo data with a n o m = 180° and with anom = 150° and are shown i n 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 120 The map from the 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 nom = 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 nickel/agarose phantoms. Region M W F ( a = 150°) M W F ( a = 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) 2 n o m Table 4.9: Percent myelin water estimated from scans with a 121 n o m n o m = 150° and a n o m = 180°. 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 a nom = 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 with « C the mean refocusing pulse flip angle with a nom n o m = 150° was 138.2° and = 180° was 156.3°. The bottom pair of images i n 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 i n 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 Region caudate (left) caudate (right) putamen (left) putamen (right) frontal (left) frontal (right) genu ic (left) ic (right) post (left) post (right) & (anom 140.0 142.3 139.9 140.8 139.8 141.5 142.2 143.4 143.4 139.1 138.4 = 150°) (1.9) (0.9) (1.3) (1.1) (1.2) (1.1) (4.4) (1.4) (2.1) (2.0) (1.8) 4.3. MULTI-EXPONENTIAL FIT a ( nom 157.5 158.6 157.2 156.6 159.3 161.4 163.6 161.1 158.1 159.1 158.7 tt = 180°) (1.6) (1.1) (1.0) (2.4) (1.5) (1.5) (3.9) (1.3) (2.7) (2.0) (1-7) Table 4.10: Refocusing pulse flip angle calculated over the regions from the scans with 150° and a n o m a n o m = = 180°. 120 100 Figure 4.11: The arithmetic mean T maps calculated from the pulse sequence with a refocusing pulse flip angle prescribed as 180° (left) and 150° (right). 2 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 The geometric mean T 180° and a M m 2 of the medium T 2 4.3. MULTI-EXPONENTIAL FIT component estimated from M E S E C with tt m = n0 = 150° are the left and right images, respectively, in Figure 4.11. The mean T 2 of the medium component was calculated from several regions within the brain and is shown in Table 4.11. The T m 2 from the two scans were within 4%. Region caudate (left) caudate (right) putamen (left) putamen (right) frontal (left) frontal (right) genu ic (left) ic (right) post (left) post (right) Table 4.11: Geometric mean T image. 2 T m 2 81.1 80.5 71.6 71.6 76.9 76.8 74.6 77.3 76.9 78.4 76.9 (ms) (2 (1 (1 (1 (1 (1 1) 5) 6) 7) 3) 1) (5 4) (0 8) (1 0) (1 9) (1 2) T m 2 (ms) 79.0 79.4 70.3 70.5 77.5 75.4 .74.1 79.3 78.4 81.3 77.7 (1.6) (1.1) (0.6) (1.4) (2.4) (1.7) (2.0) (5.2) (4.2) (4.1) (1.2) of the medium component from select regions throughout the The x ) calculated during the estimation of the T distribution, is shown in Figure 4.12. The 2 2 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 C S F is flowing and therefore may cause artifacts on the decay curves that can not be accounted for by echo formation. Second, the T of C S F is very 2 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 distribution of C S F . 2 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 x maps calculated from the pulse sequence with a refocusing pulse flip angle prescribed as 180° (left) and 150° (right). 2 Fixed Flip Angle For comparison, the T a n o m 2 distribution was calculated voxel-by-voxel from the two scans assuming = 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 m = 180°. As expected, the regions of the image with a < 180°, for example, n0 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 P U L S E S 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 Conclusions 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 from the two scans corresponded very well. The voxels 2 in the brain with C S F tended to be poorly fit (high x 2 misfit) which could be due to flow effects or an under-sampled T decay curve. 2 4.4 4.4.1 Multi-Slice, M u l t i - E c h o Introduction The current multi-echo acquisition method, M E S E , is able to collect only a single slice within 0 the given scan time, due to the hard, composite pulses. The M E S E pulse sequence can be used C to acquire multiple slices within one slab with the same S N R 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 M S 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 sequence. Six slices were acquired, within the brain, C as well as a slice acquired with the M E S E for comparison. Myelin water maps were calculated 0 for each slice of the M E S E sequence and the single slice of the M E S E . The myelin water maps C 0 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 with 4 averages and T E = 10 ms; 0 and second, a 6 slice slab, 42 echo 3 D - M E S E sequence with 1 average and T E = 10.1 ms. Other C 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 dataset and the corresponding slice of the 6 slice slab M E S E 0 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 127 0 sequence and each slice of the M E S E C 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 4.4. MULTI-SLICE, MULTI-ECHO PULSES Figure 4.14: The T E = 10 ms image from the four average, 32 echo M E S E sequence (top row). The six slices of the six slice slab M E S E 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. 0 C 129 C H A P T E R 4. NON-180 R E F O C U S I N G 4.4. PULSES 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 . 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, x d global density) which is why the two phantoms on the right side were thresholded. C 2 a n 130 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 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.3.2 4.4. MULTI-SLICE, MULTI-ECHO Regional Analysis Regions were selected in white matter and gray matter areas of the earliest echo image. The signal intensity of each region was averaged, for each echo, to calculate a T decay curve, the 2 N N L S algorithm was applied to the decay curve. The short T and medium T parameters were Q 2 2 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 sequence, was always higher C 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 decay curve acquired with the slab selected, multi-slice M E S E the T 2 C 2 sequence decayed faster than decay curve acquired with the M E S E sequence. The M E S E sequence may have subtle 0 C 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 Structure frontal_wm_l frontaLwmJ frontal_wm_l frontal_wm_r frontaLwm_r frontal_wm_r genu_wmJ genu-wmJ genu_wmJ genu_wm_r genu_wm_r genu_wm_r post_w'mJ post_wm_l post_wm_l post_wm_r , post_wm_r post_wm_r mid_wm_l mid_wmJ mid_wmJ mid_wm_r mid_wm_r mid_wm_r perivent_gm_l perivent_gm_l perivent-gmJ perivent_gm_r perivent_gm_r perivent_gm_r post_gm_l post_gm_l post_gm_l post_gm_r post_gm_r post_gm_r Sequence MESE MESE MESE MESE MESE MESE MESE MESE MESE MESE„ MESE MESE MESE MESE MESE MESE MESE MESE MESE MESE MESE MESE MESE MESE MESE MESE MESE MESE MESE MESE MESE MESE MESE MESE MESE MESE 0 C C 0 C C 0 C C C C 0 C C 0 C C 0 C C 0 C C 0 C C 0 C C 0 C C 0 C C (180) (180) (174) (180) (180) (180) (180) (180) (180) (180) (180) (173) (180) (180) (180) (180) (180) (180) (180) (180) (180) (180) (180) (178) (180) (180) (172) (180) (180) (180) (180) (180) (159) (180) (180) (172) Short p 47.258 ( 4.6%) 78.021 ( 8.0%) 83.473 ( 8.5%) 36.290 ( 3.5%) 138.347 (13.5%) 138.347 (13.5%) 94.612 ( 9.2%) 150.983 (15.2%) 150.983 (15.2%) 99.003 ( 9.3%) 179.500 (17.8%) 179.886 (17.8%) 88.007 ( 8.5%) 244.997 (24.4%) 244.997 (24.4%) 70.068 ( 6.7%) 224.328 (22.3%) 224.328 (22.3%) 177.964 (17.2%) 215.597 (22.4%) 215.597 (22.4%) 168.491 (16.7%) 228.307 (23.1%) 229.005 (23.1%) 0.000 ( 0.0%) 0.000 ( 0.0%) 0.000 ( 0.0%) 17.381 ( 1.5%) 0.000 ( 0.0%) 0.000 ( 0.0%) 0.000 ( 0.0%) 197.091 (19.1%) 139.822 (13.1%) 0.000 ( 0.0%) 199.367 (19.3%) 209.809 (20.3%) 4.4. NNLS Short T 2 15.000 36.109 34.430 15.000 24.656 24.656 22.819 37.115 37.115 15.127 35.845 35.065 17.971 35.921 35.921 15.291 40.148 40.148 29.219 39.133 39.133 29.417 30.727 30.709 — — — 15.000 — — — 42.066 22.560 — 42.229 41.925 MULTI-SLICE, MULTI-ECHO Med Med T p 931.387 (91.5%) 863.296 (88.1%) 873.974 (88.8%) 977.200 (93.7%) 885.655 (86.5%) 885.655 (86.5%) 908.342 (88.0%) 837.326 (84.4%) 837.326 (84.4%) 947.429 (89.4%) 798.541 (79.3%) 799.187 (78.9%) 942.529 (91.5%)686.287 (68.4%). 686.287 (68.4%) 808.996 (77.5%) 589.941 (58.5%) 589.941 (58.5%) 607.018 (58.8%) 500.218 (52.0%) 500.218 (52.0%) 606.233 (60.0%) 567.696 (57.3%) 567.653 (57.3%) 1161.671 (100.0%) 1196.713 (100.0%) 1200.621 (100.0%) 1173.212 (98.5%) 1163.000 (100.0%) 1163.000 (100.0%) 1018.847 (100.0%) 745.337 (72.2%) 926.926 (86.9%) 873.843 (83.6%) 607.816 (59.0%) 599.110 (58.0%) 2 79.337 80.566 81.122 79.755 83.284 83.284 77.340 78.994 78.994 - 78.345 81.258 80.791 92.041 96.742 96.742 87.281 86.406 86.406 82.458 79.342 79.342 85.955 86.775 86.831 76.568 68.768 68.533 81.559 73.643 73.643 86.499 79.251 79.282 84.274 79.542 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 and M E S E 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 < 50 ms and 50 < T < 120 ms, respectively. 0 C 2 2 133 C H A P T E R 4. NON-180 R E F O C U S I N G PULSES 100 10 1 co D •g CO CD CC / \ 0 -10 ( 200 4.4. TE (ms) 300 400 ! 1 : V i ii i i 50 100 150 200 , 4 t MULTI-SLICE, MULTI-ECHO 500 1 1 /«>^^.w-., .. ,v 1 1 250 300 1000 J 1 350 400 450 O Measured _ . M E S E Fitted (a=180) . . . . MESE Fitted (a=176.0) c 800 £ CD C O |> 1000 600 800 m u a 400 co 200 600 400 1 20 8 "0 50 100 150 200 250 40 60 » » p « 8 8 t e t a a a 300 350 400 450 TE (ms) Figure 4.17: Region (genu_wm_l) shown on the top figure, T decay curves acquired by the M E S E and M E S E (middle figure). N N L S fits to the M E S E data, one assuming the refocusing pulse flip angle was 180° and the other where the optimal flip angle was calculated (bottom figure). 2 C 0 C 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.3.3 4.4. MULTI-SLICE, MULTI-ECHO Flip Angle Estimation The slice from the slab M E S E that corresponded to the slice from the single-slice M E S E was C 0 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 (mwm ) and 2) estimating the refocusing 180 pulse flip angle per voxel (mwm ). e Figure 4.18: Maps of the refocusing pulse flip angle calculated from the 6 slice, 42 echo pulse sequence M E S E . The two outside slices (top left and bottom right) were not expected to be reasonable as they had significant (as expected) aliasing artifacts. C 4.4.4 Conclusions Moving from single-slice to multiple slices (slab selection) per acquisition will greatly improve the utility of the myelin water mapping. To acquire multiple slices per acquisition slice-selective refocusing pulses are necessary. The flip angle of slice-selective refocusing pulses is inaccurate and must be accounted for during 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 . The T estimated from the decay curves of regions within the six 0 2 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 ), the refocusing pulse ( a ) and the tissue parameters, p, T , T i . The goal of this work is to e r 2 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 , T i , p and the pulse sequence parameters T E , a and a , the decay 2 e r 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] and was then rotated by an angle a based on the rotation matrix: T cos (a/2) sin (a/2) sin (a) 0 sin (a/2) cos (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) 2 2 R(cx) = 4.5.2 2 2 (4.19) Simulation Methods Three sets of simulations were done to determine the effect of knowing or not knowing the excitation pulse flip angle a . The first simulation (Section 4.5.3) had an excitation pulse that varied e and the decay curve parameters were calculated assuming the excitation pulse was a = 90° and e 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 a e assumed the excitation pulse a e but r = 90°. The third simulation (Section 4.5.5) had the excitation and refocusing pulses varied (but still a = a /2) e 4.5.3 = a /2 r and calculate both a and a . e T Simulation 1: Inaccurate Excitation and Accurate Refocusing, Assume both Accurate This simulation will assess the effect of having an inaccurate excitation pulse but a perfect refocusing pulse a = 180°. r Methods A set of true decay curves were calculated with the excitation pulse flip angle set' to a — e 90°,85°,...,50°. Other parameters were: a r = 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 S N R = 200 were added to the true decay curve. For each decay curve with noise, a r was calculated, as well as, p, T and x 2 2 (the misfit). Results and Discussion This simulation assessed the parameter estimation when o n l y 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 a e ^ 90° as when a for p. The average p decreased as a function of s i n ( a ) . The x e e 2 = 90°, except misfit remained constant at approximately 27, therefore the decrease in p corresponds very well to a decrease in a . e This would imply that if the only a or p changed, they could not be distinguished. e 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 True a e (°) 90 85 80 75 70 65 60 55 50 a r 176 176 176 176 176 176 176 176 175 4.5. NON 90° EXCITATION PULSE T ( ms) 79.2 (2.3) 79.2 (2.0) 79.1 (2.0) 79.0 (1.7) 79.0 (1.9) 79.1 (2.6) 79.0 (2.7) 79.0 (2.6) 78.9 (3.5) P 1011.2 (17.2) 1008.9 (18.8) 996.4 (16.8) 978.3 (18.2) 951.4 (16.7) 917.4 (17.5) 877.6 (18.3) 830.4 (17.7) 777.8 (17.4) (°) (4.3) (4.5) (4.5) (4.5) (4.5) (5.0) (4.8) (4.6) (5.2) 2 X 27.9 (7.5) 27.6 (7.5) 27.9 (7.5) 27.4 (7.0) 27.7 (7.7) 27.8 (7.7) 26.8 (6.9) 27.3 (7.2) 26.8 (7.0) 2 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 i m u l a t i o n 2: Inaccurate E x c i t a t i o n a n d Inaccurate R e f o c u s i n g , Assume Accurate Excitation This simulation will assess the effect of having a e = a /2 T but assuming a e = 90° for the curve fitting. This simulation closer to reality when acquiring multi-exponential data, as spatial nonuniformities in the transmit power will affect both the excitation and refocusing pulses. Methods The parameters for the simulation are: 32 echoes, T j = 1000 ms, T = 80 ms, T R = oo, echo 2 spacing = 1 0 ms and p — 1000. The excitation pulse was varied: a e a = 2a in every case). For each setting of a -a , r e e = 90°, 8 5 ° , . . . , 50° (and 500 realizations of quadrature noise (Gaussian r on two channels) with S N R = 200 were added to the true decay curve. For each decay curve with noise, a was calculated as well as p, T r and x 2 2 (the misfit). The excitation pulse flip angle was assumed to be a = 90°. e Results and Discussion The estimated parameters are shown in Table 4.14. The refocusing pulse flip angle a T was estimated reasonably accurately, but all other parameters had problems. As in the first simulation, p decreased almost as a function of sin (a ). e increased as a e As well, T decreased as a decreased. The x misfit 2 2 e 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 a /a e r a (°) (°) T (ms) P 176 (4.3) 1011.2 (17.2) 173 (6.0) 1004.4 (17.3) 167 (9.6) 983.5 (16.5) 151 (6.3) 974.5 (20.7) 138 (2.3) 955.0 (15.2) 127 (1.9) 925.7 (9.2) 117 (1.7) 902.2 (8.5) 107 (1.4) 871.6 (9.0) 852.8 (10.2) 97 (1.3) r 90/180 85/170 80/160 75/150 70/140 65/130 60/120 55/110 50/100 x 27.9 (7.5) 27.9 (7.6) 38.2 (9.2) 57.9 (11.7) 77.0 (16.0) 92.1 (17.4) 103.1 (18.7) 113.6 (18.7) 123.1 (20.4) 2 2 79.2 (2.3) 79.6 (2.4) 80.1 (2.4) 79.4 (2.0) 78.9 (2.7) 78.6 (2.0) 77.2 (1.3) 76.3 (1.7) 74.0 (0.9) Table 4.14: The excitation and refocusing pulses were varied with a = 90°, 8 5 ° , . . . , 50° and a = 2a . The curve fitting calculated a , T , p and x assuming a = 90°. The true T = 80 ms and p = 1000. e 2 r e r 2 e 2 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 estimation. This simulation will assess the effect of having a = a /2 2 e r and calculating both a and a along with the other parameters. e r Methods The parameters for the simulation were: 32 echoes, T i = 1000 ms, T = 80 ms, T R = oo, echo 2 spacing = 10 ms and p = 1000. The excitation pulse was varied: a e = 90°, 8 5 ° , . . . , 50° (and a = 2a in every case). For each setting of a , 500 realizations of quadrature noise (Gaussian on r e e two channels) with S N R = 200 were added to the. true decay curve. For each decay curve with noise, the a , a , p, T and % (the misfit) were calculated. The decay curve was fit assuming 2 r the ratio a /a e e r 2 = 1/2. Results and Discussion A l l parameters (shown in Table 4.15) were estimated accurately, including a and a . The stane r dard deviation of the estimated T and p increased by a factor of two as a = 90° to a = 50°. 2 e e 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 True a /a e (°) r 90/180 85/170 80/160 75/150 70/140 65/130 60/120 55/110 50/100 a (°) 88 (2.4) 86 (2.3) 80 (1.7) 76 (0.8) 70 (1.0) 65 (0.6) 60 (0.7) 55 (0.5) 50 (0.6) S a r 176 172 161 152 140 130 120 110 100 (°) (4.8) (4.5) (3.3) (1.6) (2.0) (1.1) (1.5) (1.0) (1.3) 4.5. NON 90° EXCITATION PULSE P 1013.1 (17.8) 1007.2 (18.0) 1004.9 (17.4) 998.1 (14.7) 1013.9 (24.8) 1008.3 (18.2) 1011.5 (23.5) 1005.4 (20.2) 1011.3 (28.3) T 2 79.2 79.3 79.5 79.7 78.8 79.3 79.1 79.5 79.3 (ms) (1.8) (2.1) (1.8) (1.7) (1.7) (2.7) (2.4) (3.0) (3.6) x 1 28.0 27.7 27.8 28.1 27.9 27.6 27.5 28.0 26.9 (7.6) (7.6) (7.5) (7.6) (7.4) (7.3) (7.6) (7.5) (6.8) Table 4.15: The excitation and refocusing pulses were varied with a = 90°, 8 5 ° , . . . , 50° and a = 2a . The curve fitting calculated a , a , T , p and x assuming a /a = 1/2. e 2 r 4.5.6 e e r 2 e r 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 , when the excitation pulse flip angle was assumed to be a = 90°, was as accurate 2 e and precise (Table 4.14) as the T estimated when the excitation pulse flip angle was estimated 2 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.5.7 4.6. NON-CONSTANT T R A I N Conclusions The excitation pulse flip angle, as well as the refocusing pulse flip angle, is affected by the inhomogeneity 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 a = 90° and a = 180°. e 4.6 r 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 S N R will lead to increased variability of the measured decay curve parameters. The short T is of primary interest and is determined from the signal in the decay curve with 2 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 S N R 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 s Q 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 • s L = ' i . o , . r., i.o i} 16 • s i = 0.5,. 2) , 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, s when testing , is a pulse train of 90° refocusing pulses. Both of these were used for reference ,(4) and s$. The third set of scalings (Figure 4 . 1 9 ) , is an amalgamation of the a first two, it retains 1 8 0 ° refocusing pulses for the first part of the train to retain the higher S N R and the last part of the train is 9 0 ° to decrease the power deposition. The fourth set of scalings (Figure S Q \ 4.20), was defined as Sa\i) for = exp(—i/10) i = 0 , 1 , . . . , 15. The denominator within the exponential was chosen so the decay in the scalings was not too fast. This last one resulted in a smooth decrease in the amplitude of the refocusing pulses throughout the train. 1 1 1 1 1 A""." A—A . A — - j x — / L (1 -4 U J m m L4 1 1 „i i 1 I : RF G fl G r fl fl u M M II • ... fl fl Ml 1 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 • ! 1 fl G L -4 I r 1 1 NON-CONSTANT TRAIN 1 M /A—„..A.. . A — - A \~; A RF 4.6. u L-4 f--4 „ 1 I I , •fl fl M i--4 l_-4 I I l i „ I I M Ml v ../ DACQ \ v - • 10 • 20 30 40 50 60 70 Time (ms) 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 and flip angle from a 2 decay curve that was "collected" with non-constant refocusing pulses through the train. 4.6.1.1 Methods I created a set of simulations with p = 1000, T 2 = 80 ms, a e = 90°. The four sets of refocus- ing pulse trains, defined above, were used. For each decay curve with noise, the decay curve parameters were calculated and compared to truth. For the decay curves with non-constant refocusing pulses, the relative scale between each pulse was assumed to be known. For example, the relative scale between pulses in the third refocusing pulse train was assumed to be 1.0,1.0,1.0,0.5,... , 0.5. 4.6.1.2 Results The parameters a, p, and T were estimated from the simulated data for each of the scalings and 2 are shown in Table 4.16. There was negligible difference between refocusing pulse train s|, ' and 3 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 SQ ^ compared to the constant refocusing pulse train of a = 1 8 0 ° . As expected, p and T had a (2) higher standard deviation in simulation s . r 2 a Simulation sk x a p T2 (ms) 176 (4.5) 1010.6 (19.6) 78.7 (2.6) 13.0 (5.4) r 2 90 (1.4) 1005.2 (27.3) 78.4 (5.5) 12.8 (4.9) ] 176 (4.5) 1010.6 (19.8) 78.8 (2.8) 12.6 (4.8) sL 176 (4.4) 1010.0 (20.1) 78.8 (2.9) 12.8 (5.1) 2) sa ( 4) 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, a 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. r 4.6.2 Scanning The four sets of scalings were implemented i n 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 refocusing 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 i n a file on the scanner computer. A t the beginning of the s c a n O routine, the amplitude of the i th 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 i n Figures 4.21 and 4.22. The fitted decay curve was very close to the measured decay curve for both s^) and 3 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: 18016 a: 9016 800 600 ^600 CO c CD | 400 co c a) ^ 200 0 "0 50 100 150 TE (ms) a: 180, 162.9, 200 50 100 150 TE (ms) a: 180, 180, 180, 90 40 800 800 ^600 ^600 c 5 400 400 cvj c g> CO c \\ c ^ 200 200 0 "0 13 c/> c CO 5 200 50 100 150 TE (ms) 200 0 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° (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. N 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: 18016 a: 90 16 600 350 500 300 | 250 -t—t cn c 400 CD I 200 _c To 150 c I * _c 300 "co o)200 100 CO 100 0 "0 50 50 100 TE (ms) a: 180, 162.9, 150 00 200 600 500 500 >^ "w 400 c 'co 400 c TE (ms) 150 200 13 4 -—* CD • 100 a: 180, 180, 180, 90 40 600 | 50 0 c 300 IB S)200 300 "CO D)200 cn CO 100 0 "0 100 50 100 TE (ms) •(t-C^-JHll 150 200 0 50 100 TE (ms) 150 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 T R A I N The T of each phantom was calculated from each pulse sequence and is shown in Table 4.17. 2 The exponentially decaying refocusing pulse train typically underestimated the T 2 by approxi- mately 10%. Phantom bottomJeft bottomJeft bottomJeft bottomJeft bottom_middle bottom_middle bottom_middle bottom_middle bottom_right bottom_right bottom_right bottom_right topJeft topJeft topJeft topJeft Sequence s (3) S ( ) 2 (ms) 24.6 20.9 24.6 21.2 83.2 70.0 82.8 75.0 23.0 17.3 22.5 20.3 299.2 109.9 374.9 280.4 4 s(D S( ) ( ) 2 3 S S( ) 4 (D S T S( ) 3 S( ) 4 W S S( ) 2 S( ) 3 S( ) 4 x 2.0 1.0 7.0 2.8 1.0 0.6 38.9 7.7 1.5 1.1 3.4 3.6 0.4 1.4 58.1 16.0 2 162.0 81.0 154.0 142.0 157.0 85.0 156.0 145.0 166.0 80.0 156.0 148.0 156.0 83.0 155.0 146.0 Table 4.17: The T and refocusing pulse flip angle (for the first echo) calculated from each train of refocusing pulses. The T , from the M E S E sequence were: 25 ms (bottom, left), 85 ms (bottom, middle), 24 ms (bottom, right), and 280 ms (top, left). 2 2 0 One thing that was a little odd was the curves for s^ ) did not follow as well as I would 4 have expected. I noticed from previous work that the calculated refocusing pulse flip angle given ao = 180° is approximately a i = 157° and given ao — 90° is approximately ft ] c a c ca 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 P U L S E S - 4.6. S90 x 0.50 0.53 0.55 0.57 0.60 0.62 0.65 3.434 2.551 1.851 1.387 0.987 0.735 0.603 22.47 22.39 22.38 22.48 22.50 22.54 22.74 0.68 0.596 22.81 0.70 0.72 0.75 0.647 0.780 0.969 23.00 23.21 23.53 2 T 2 NON-CONSTANT T R A I N (ms) Table 4.18: The resulting x misfit given the scaling of the R F pulses was 1.0,1.0,1.0,590 from s^ ). 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 x is in bold. 2 4 2 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.6.3 4.7. CONCLUSIONS Conclusion The amplitude of the R F pulses in the train can be scaled to have any amplitude and the T tends 2 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 correction due to B i inhomogeneities have relied on the acquisition of a B i map or 2 based on uniform phantoms. Other methods required two acquisitions to calculate the corrected mono-exponential T , one acquisition to calculate the multi-echo decay curve and the second 2 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 and a. Therefore, the T and B i are inherently 2 2 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 distributions of guinea pig brain at 1.5T and 4T and found both the short and medium T 2 2 peaks had a faster decay at the higher field strength. The decrease in T at higher field would require 2 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 field systems and will account for the echo artifacts. 149 2 measurements to be done on high 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 will likely be other filters which could prove to be better. The primary goals for a noise reduction filter, applied to myelin water quantification, should include: ability to smooth over homogeneous areas and retain the boundaries between the tissues, produce little or no bias in the mean myelin water fraction over a region, and reduce the variability of the myelin water fraction within a local region. 5.2 Linear Combination The linear combination method to calculate the myelin water fraction was shown to be robust and very fast, in addition to its extreme ease of implementation after the coefficients are determined. The short T filter was to select for T near 20 ms only, and suppressed all T above 70 ms. The 2 2 2 myelin water maps were comparable to myelin water maps calculated by the N N L S algorithm. The linear combination method to analyze T from a multi-echo M R I dataset is restricted primarily to 2 150 C H A P T E R 5. CONCLUSIONS 5.2. LINEAR COMBINATION the short T range. The restriction is due to the broad selection if longer T are of interest. The 2 2 broad selection would negate the use of linear combination for T longer than 50 ms. Therefore, 2 the curve fitting algorithms which provide information about the T distribution are still required 2 if the T of interest is greater than 50 ms. 2 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 nonorthogonal as a function of T . Therefore, it may not be possible to create a better filter than 2 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 < 70 ms is less sharp than the filter for 10 < T < 20 ms. 2 2 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 toward higher 2 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. 5.3 CONCLUSIONS 5.3. Bi/T 2 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). algorithm should be able to estimate the T 2 Therefore, a modified T 2 curve fitting 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 , to estimate the T distribution and a from a multi-exponential a 2 decay curve. Noise simulations showed that the accuracy and precision were similar to the standard method to estimate T . A decrease in the precision of the estimation was shown to be a 2 function of the decreased a. Several extensions to the method were discussed: slab selected multi-slice, train of nonconstant 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 important 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 second strength of the method is the curve fit algorithm is able to estimate the T from a train of 2 non-180 refocusing pulses. A pulse sequence, with a train of non-180 refocusing pulses deposits 0 0 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 of the short and medium components. 2 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. CONCLUSIONS 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 S N R is to do multiple averages. Decreasing the scan time, without a large decrease in the S N R , would increase the clinical feasibility of myelin water imaging. The standard assumption in quantitative T is the pulse sequence must have a long T R 2 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 r e f e r r e d T2 Q u a n t i f i c a t i o n 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 spinecho 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 S N R 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, multiecho 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, A p r 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, A p r i l 2000. [12] V . Vasilescu, V . Simplaceanu E . Katona and, and D . Demco. Water compartments i n 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 relaxation 2 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. Spinspin 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, J u n 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 distributions and water contents in normal human brain. Magn. Reson. 2 Med., 37:34-43, 1997. [20] C Beaulieu, F R Fenrich, and P S 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 multicomponent T 2 relaxation measurements with histopathologic correlation in an experimental model of M S . 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 quantitation for clinical applications. J. Magn. 2 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 i n T estimation using multislice multiple-echo 2 imaging. Magn. Reson. Med., 4:34-47, 1987. [36] I. M . Vavasour, K . P. Whittall, D. K . L i , and A . L . MacKay. Different magnetization transfer effects exhibited by the short and long T components i n human brain. Magn Reson Med, 2 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 . Y u . 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 . Whittall 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. Whittall 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 relaxation data. Magn. Reson. Med., 35(3):370-378, 1996. 2 157 Bibliography [45] K . P. Whittall, A . L . MacKay, and D . K . B . L i . Are mono-exponential fits to a few echoes sufficient to determine T relaxation for in vivo human brain? Magn. Reson. Med., 41 (6): 12552 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 . Whittall, J . Adler, D . K . B . L i , D . W . Paty, and D . Graeb. alization of myelin water in brain by magnetic resonance. In vivo visu- 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 vessel 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, 4 2 ( 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. P h D 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 S N R . 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 methods 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. T 2 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 maps of phantoms. 2 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 A p r 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. Magn. Reson. A, 109:117-120, July 1994. 160 J. Bibliography [80] C . Seong, D . Jones, W . E . Reddick, R . J . Ogg, and R . G . Steen. Establishing norms for agerelated 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. Magn Reson Med, Imaging of the active B l field in vivo. 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. ments of multi-component T2 relaxation behaviour in guinea pig brain. In vivo measure- Magn. Reson. Imag., 17(9):1319-1325, 1999. [86] S. Meiboom and D . G i 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., Whittall, K . P . , MacKay, A . L . , "Validation of Anisotropic Filtering for Myelin Water Maps". Proceedings of the International Society of Magnetic Resonance in Medicine. Glasgow, Scotland, A p r i l 21-27, (2001). Page 820. • Chapter 3: Myelin Water Fraction from Linear Combination - Paper: Jones, C.K., Xiang, Q.S., Whittall, K . P . , Mackay, A . L . , "Short T Selection 2 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 and 2 B\ from Decay Curves Collected with non-180 Refocusing Pulses", Magn Reson Med 0 (submitted), 2003. - Presentation: Jones, C.K., Xiang, Q.S., Whittall, K . P . , Mackay, A . L . , "Calculating T and B\ from Decay Curves Collected with non-180 Refocusing Pulses", Proceedings 0 2 of the International Society of Magnetic Resonance in Medicine. Toronto, O N , July 10-16, (2003). 162
- Library Home /
- Search Collections /
- Open Collections /
- Browse Collections /
- UBC Theses and Dissertations /
- T2 decay curve acquisition and analysis in MRI noise...
Open Collections
UBC Theses and Dissertations
Featured Collection
UBC Theses and Dissertations
T2 decay curve acquisition and analysis in MRI noise considerations, short T2, and B1 field encoding Jones, Craig K. 2003-12-31
pdf
Page Metadata
Item Metadata
Title | T2 decay curve acquisition and analysis in MRI noise considerations, short T2, and B1 field encoding |
Creator |
Jones, Craig K. |
Date | 2003 |
Date Issued | 2009-11-13T18:21:52Z |
Description | The myelin sheath is a lipid bilayer wrapped around axons in the human brain. The myelin allows faster conduction and uses less energy than along non-myelinated axons. Diseases, such as multiple sclerosis, break down the myelin sheath resulting in cognitive and physical disability. Water trapped between the myelin bilayers has a T2 ≈ 15 ms compared to intra/extra-cellular water with a T2 ≈ 80 ms and cerebrospinal fluid which has a T2 ≈ 2000 ms. A multi-echo MRI pulse sequence is used to acquire a T2 decay curve which is fit using a non-negative least-squares (NNLS) curve fitting algorithm to calculate the signal as a function of the T2 (T2 distribution). The myelin water fraction (MWF) is calculated as the signal with a T2 < 50 ms divided by the total signal in the T2 distribution. Quantification of the MWF in vivo, is slow, prone to errors due to artifacts in the T2 decay curve, and requires many data acquisition averages to attain the necessary high signal-to-noise ratio. Each of these areas were addressed in this thesis. First, I compared the accuracy and precision of MWF estimates from a set of simulated and in vivo multi-echo data with and without noise reduction filters applied. The MWF estimated from anisotropically filtered multi-echo data had the smallest variability, over homogeneous regions, compared to MWF estimates from other filtered data. Second, I created a new technique to optimize coefficients that when linearly combined with multi-echo data, result in fast, accurate estimates of the MWF. Simulations and in vivo estimates of the MWF were as accurate as those from the NNLS algorithm and 20,000 times faster. Finally, the standard multi-echo sequence was optimized to remove artifacts due to inaccurate refocusing pulses. I created a novel T2 curve fitting algorithm, based on NNLS, to accurately estimate the T 2 distribution and refocusing pulse flip angle from a T2 decay curve. Simulations and data acquired on phantoms showed the new technique was as accurate as the NNLS algorithm in quantifying the T2 distribution parameters. In vivo myelin water maps were as good as those estimated from the optimized multi-echo pulse sequence. |
Extent | 8349221 bytes |
Genre |
Thesis/Dissertation |
Type |
Text |
File Format | application/pdf |
Language | eng |
Collection |
Retrospective Theses and Dissertations, 1919-2007 |
Series | UBC Retrospective Theses Digitization Project |
Date Available | 2009-11-13 |
Provider | Vancouver : University of British Columbia Library |
Rights | For non-commercial purposes only, such as research, private study and education. Additional conditions apply, see Terms of Use https://open.library.ubc.ca/terms_of_use. |
DOI | 10.14288/1.0085721 |
URI | http://hdl.handle.net/2429/14941 |
Degree |
Doctor of Philosophy - PhD |
Program |
Physics |
Affiliation |
Science, Faculty of Physics and Astronomy, Department of |
Degree Grantor | University of British Columbia |
Graduation Date | 2003-05 |
Campus |
UBCV |
Scholarly Level | Graduate |
Aggregated Source Repository | DSpace |
Download
- Media
- [if-you-see-this-DO-NOT-CLICK]
- ubc_2003-859312.pdf [ 7.96MB ]
- [if-you-see-this-DO-NOT-CLICK]
- 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
- Original Record: 1.0085721 +original-record.json
- Full Text
- 1.0085721.txt
- Citation
- 1.0085721.ris
Full Text
Cite
Citation Scheme:
Usage Statistics
Country | Views | Downloads |
---|---|---|
United States | 49 | 3 |
China | 28 | 2 |
France | 24 | 2 |
Canada | 17 | 4 |
Finland | 15 | 2 |
Russia | 9 | 0 |
United Kingdom | 6 | 0 |
Germany | 5 | 7 |
Japan | 5 | 0 |
Spain | 2 | 0 |
India | 2 | 0 |
Sweden | 2 | 92 |
Italy | 2 | 0 |
City | Views | Downloads |
---|---|---|
Unknown | 41 | 64 |
Hangzhou | 18 | 0 |
Oulu | 15 | 2 |
Mountain View | 12 | 0 |
Champigneulles | 11 | 1 |
Edmonton | 9 | 1 |
Ashburn | 7 | 0 |
Penza | 6 | 0 |
Shenzhen | 5 | 1 |
Bethesda | 4 | 0 |
Vancouver | 4 | 2 |
Jackson | 4 | 0 |
Alexandria | 3 | 0 |
{[{ mDataHeader[type] }]} | {[{ month[type] }]} | {[{ tData[type] }]} |
Share
Embed
Customize your widget with the following options, then copy and paste the code below into the HTML
of your page to embed this item in your website.
<div id="ubcOpenCollectionsWidgetDisplay">
<script id="ubcOpenCollectionsWidget"
src="{[{embed.src}]}"
data-item="{[{embed.item}]}"
data-collection="{[{embed.collection}]}"
data-metadata="{[{embed.showMetadata}]}"
data-width="{[{embed.width}]}"
async >
</script>
</div>
Our image viewer uses the IIIF 2.0 standard.
To load this item in other compatible viewers, use this url:
http://iiif.library.ubc.ca/presentation/dsp.831.1-0085721/manifest