University of Southern California Dissertations and Theses

High-dimensional magnetic resonance imaging of microstructure

(USC Thesis Other)

High-dimensional magnetic resonance imaging of microstructure

PDF

Download a page range

Download transcript

Copy asset link

Request this asset

Transcript (if available)

Content HIGH-DIMENSIONAL MAGNETIC RESONANCE IMAGING OF MICROSTRUCTURE by DAEUN KIM A dissertation submitted to the USC graduate school in partial fulllment of the requirements for the degree of Doctor of Philosophy in Electrical and Computer Engineering University of Southern California Los Angeles, California August 2019 Doctoral Committee: Associate Professor Justin P. Haldar, Chair Professor Krishna S. Nayak Professor John C. Wood c 2019 Daeun Kim. All rights reserved. Abstract MRI is a powerful and widely used biomedical imaging technique, but it suers from limited spatial resolution. This prevents MRI from being used to directly probe microscopic tissue information. Fortunately, indirect estimation of microscopic tissue features is still possible by modeling each voxel as a mixture of dierent microstructural components with distinct signal characteristics, and encoding microstructural information by manipulating MRI physics. However, the conven- tional indirect approaches are often associated with a highly ill-posed inverse problem, which can severely limit the ability to accurately unmix the microstructural components. This dissertation proposes and evaluates a new high-dimensional MRI approach that combines multidimensional spa- tial encoding with multidimensional contrast encoding (e.g., simultaneous diusion encoding and relaxation encoding) and multidimensional regularized estimation to provide substantially improved microstructural unmixing capabilities. The new approach can be justied theoretically, and this dissertation demonstrates an estima- tion theoretic analysis that shows the multidimensional approach has substantial advantages over conventional lower-dimensional approaches. This estimation theoretic framework is also shown to enable optimized experiment design. The dissertation also evaluates the new approach in two dier- ent congurations (one using diusion-relaxation contrast encoding and the other using relaxation- relaxation contrast encoding) in a wide range of experiments, including numerical simulations and real experiments with phantoms, ex vivo animals, and in vivo humans. The results demonstrate that the proposed approach provides unique information that is not available with conventional methods, and suggests that the proposed approach has the potential to have a high impact if translated to practical applications in medicine, neuroscience and biology. Acknowledgments Foremost, I would like to express the deepest appreciation to my advisor, Professor Justin P. Haldar for his consistent guidance and invaluable advice throughout my graduate studies. Your enthusiasm for research has always been contagious, which inspired me to investigate something new, dig it into, and nally develop my philosophy. Thank you for all your support and encouragement on my research and career. I would like to thank my doctoral committee, Professors Krishna S. Nayak and John C. Wood for their time, eort and insightful questions. I would also like to thank the other two members of my qualifying exam committee, Professors Richard M. Leahy and Antonio Ortega for their contribution and thoughtful comments. I am also grateful to people at the Biomedical Imaging Group Dr. Anand A. Josh, Dr. Syed Ashrafulla, Dr. Yanguang Lin, Dr. Wentao Zhu, Dr. Sergul Aydore, Dr. Chitresh Bhushan, Dr. Divya Varadarajan, Dr. Minqi Chong, Dr. Jian Li, Tae Hyung Kim, Soyoung Choi, Rodrigo Lobos, Hossein Shahabi, Yunsong Liu, Haleh Akrami, Yijun Liu and Jiayang Wang for their valuable and stimulating discussions. I am also thankful for collaborators Dr. Jessica L. Wisnowski, Dr. Christopher Nguyen, Dr. Joong Hee Kim and Dr. Bo Zhao for their support that has made my research more productive. I would also like to express my gratitude to my former advisors and mentors. I would like to thank Dr. Jongho Lee and Dr. Jeongtae Kim for their support, encouragement and mentorship that led me to pursue my Ph.D. degree. I would also like to thank two pioneers, Dr. Zang-Hee Cho and Dr. Seiji Ogawa for their help and encouragement that made me dive into the MRI research eld. Finally, I would like to thank my family for all their love, support and encouragement. I am grateful to my parents and my parents-in-law for their sacrices so that I could have good achievements. I am also thankful for my brother, Hyo-Eun Kim for cheering me up. Special thanks to my beloved husband, Soowang Park for his all time support, devotion and trust that allowed me to overcome tough times and continue my Ph.D. degree. To my adorable baby Jenna Park, I am grateful that you make me always happy and smile. Thank you everyone who helped me get to this nal stage of my Ph.D., although I can't list all of their names. iii Contents Abstract ii Acknowledgments iii List of Figures vii List of Tables xiii Abbreviations xiv 1 Introduction 1 1.1 Problem Statement . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 1 1.2 Motivation . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 2 1.3 Main Results . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 4 1.4 Organization of the Dissertation . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 5 2 Background 7 2.1 Physics of Nuclear Magnetic Resonance . . . . . . . . . . . . . . . . . . . . . . . . . 7 2.1.1 Nuclear Spin System and Magnetization . . . . . . . . . . . . . . . . . . . . . 7 2.1.2 RF Excitation . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 10 2.1.3 Relaxation . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 13 2.1.4 The Bloch Equation . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 15 2.2 Basic Principles of MR Imaging . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 16 2.2.1 The Bloch Equation with Gradient Fields . . . . . . . . . . . . . . . . . . . . 17 2.2.2 Signal Detection . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 17 2.2.3 Fourier Interpretation . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 20 2.3 Basic Imaging Methods . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 20 2.3.1 Free Induction Decay (FID) . . . . . . . . . . . . . . . . . . . . . . . . . . . . 20 2.3.2 Spin-Echo Imaging . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 23 2.3.3 Inversion-Recovery Imaging . . . . . . . . . . . . . . . . . . . . . . . . . . . . 28 2.3.4 Diusion-weighted Imaging . . . . . . . . . . . . . . . . . . . . . . . . . . . . 32 2.3.5 High-dimensional Contrast Encoding . . . . . . . . . . . . . . . . . . . . . . . 34 3 Multidimensional Correlation Spectroscopic Imaging of Exponential Decays 37 iv Contents v 3.1 Overview . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 37 3.2 Compartment Modeling Background . . . . . . . . . . . . . . . . . . . . . . . . . . . 40 3.3 From Correlation Spectroscopy To Correlation Spectroscopic Imaging . . . . . . . . 43 3.3.1 Estimation Theoretic Analysis of Encoding in Higher Dimensions . . . . . . . 43 3.3.2 Correlation Spectroscopic Imaging . . . . . . . . . . . . . . . . . . . . . . . . 46 3.3.2.1 Spatial-Spectral Modeling . . . . . . . . . . . . . . . . . . . . . . . . 46 3.3.2.2 High-dimensional Correlation Spectroscopic Image Estimation . . . 48 3.4 Summary . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 50 4 Diusion-Relaxation Correlation Spectroscopic Imaging (DR-CSI) 52 4.1 Diusion-Relaxation Correlation Modeling . . . . . . . . . . . . . . . . . . . . . . . . 52 4.2 Experiments . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 54 4.2.1 Numerical Simulations . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 54 4.2.2 Phantom Experiments . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 56 4.2.3 Ex Vivo Mouse Spinal Cord Experiments . . . . . . . . . . . . . . . . . . . . 62 4.3 Discussion . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 69 4.4 Summary . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 70 5 T 1 Relaxation - T 2 Relaxation Correlation Spectroscopic Imaging (RR-CSI) 72 5.1 T 1 Relaxation - T 2 Relaxation Correlation Modeling . . . . . . . . . . . . . . . . . . 72 5.2 Experiments . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 73 5.2.1 Numerical Simulations . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 74 5.2.2 Pumpkin Experiments . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 80 5.2.3 Orange Experiments . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 83 5.2.4 In Vivo Human Brain Experiments . . . . . . . . . . . . . . . . . . . . . . . . 86 5.3 Discussion . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 94 5.4 Summary . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 96 6 Toward Improved Eciency 98 6.1 Optimized Experiment Design Approach . . . . . . . . . . . . . . . . . . . . . . . . . 98 6.1.1 CRLB for the Multidimensional CSI Model . . . . . . . . . . . . . . . . . . . 99 6.1.2 Optimal Experiment Design . . . . . . . . . . . . . . . . . . . . . . . . . . . . 99 6.1.3 Application Examples . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 100 6.1.3.1 DR-CSI Experiment . . . . . . . . . . . . . . . . . . . . . . . . . . . 100 6.1.3.2 RR-CSI Experiment . . . . . . . . . . . . . . . . . . . . . . . . . . . 105 6.2 RR-CSI with a MR Fingerprinting Acquisition . . . . . . . . . . . . . . . . . . . . . 108 6.2.1 In Vivo Human Brain Experiments . . . . . . . . . . . . . . . . . . . . . . . . 108 6.3 Summary . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 114 7 Conclusion 115 A ADMM Based Solution of Eq. (3.11) 118 Contents vi B ADMM Based Solution of Eq. (3.12) 121 Bibliography 123 List of Figures 2.1 Nuclear spins and magnetic moments. (a) A nuclear spin with nonzero magnetic moment. (b) Nuclear spins that are randomly oriented in the absence of an external eld. (c) Nuclear spins that are aligned in the presence of an external eld. . . . . . 8 2.2 Precession of the magnetization about the applied magnetic eld. . . . . . . . . . . . 11 2.3 Behavior of the net magnetization in (a) the laboratory frame (i.e., the static frame) and (b) the Larmor-rotating frame. In the rotating frame, the magnetization is tipped into the transverse plane by the application of the RF magnetic eld. . . . . 13 2.4 Relaxation curves in (a) the longitudinal direction and in (b) the transverse plane. . 14 2.5 Block diagram of the quadrature phase-sensitive detection. . . . . . . . . . . . . . . 19 2.6 A simulated FID signal from Eq. 2.54. The signal envelop decays with a time constant T 2 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 22 2.7 Graphical description of (a) 90 excitation and (b) dephasing of isochromats in a spin system after the excitation in the rotating frame. . . . . . . . . . . . . . . . . . 22 2.8 Graphical description of spin echo generation. . . . . . . . . . . . . . . . . . . . . . . 23 2.9 Generation of a spin echo signal. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 24 2.10 The saturation-recovery pulse sequence and the time evolution of the longitudinal magnetization. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 25 2.11 Spin-echo imaging sequence. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 26 2.12 The CPMG pulse sequence. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 27 2.13 T 2 -weighted brain images acquired using a CPMG-based imaging sequence with TR=5000 ms and 15 TEs ranging from 7.5 ms to 217.5 ms in 15 ms increments. Four representative images are displayed and each image is acquired with (a) TE = 22.5 ms, (b) TE = 67.5 ms, (c) TE = 112.5 ms and (d) TE = 217.5 ms. . . . . . . 28 2.14 The inversion-recovery pulse sequence and the time evolution of the longitudinal magnetization. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 29 2.15 Inversion-recovery (IR) spin-echo (SE) imaging sequence. . . . . . . . . . . . . . . . 31 2.16 T 1 -weighted brain images acquired using a IR-SE imaging sequence with TR=5000 ms and TE = 7.5 ms. Each image is acquired with (a) TI = 100 ms, (b) TI = 200 ms, (c) TI = 400 ms and (d) TI = 700 ms. . . . . . . . . . . . . . . . . . . . . . . . 32 2.17 Stesjkal-Tanner PFG diusion-weighted spin-echo sequence. . . . . . . . . . . . . . . 33 2.18 Inversion-Recovery Multi-echo Spin-Echo (IR-MSE) imaging sequence. . . . . . . . . 36 2.19 Diusion-weighted Spin-Echo imaging sequence. . . . . . . . . . . . . . . . . . . . . . 36 vii List of Figures viii 3.1 The framework of the proposed multidimensional correlation spectroscopic imaging approach. While a 4D spectroscopic imaging with diusion and T 2 relaxation is shown in the gure, other contrast mechanisms are possible even in higher dimension using the same framework. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 40 3.2 Toy illustration of the advantage of 2D multidimensional correlation spectroscopy over 1D spectroscopy. Ground truth values of three spectral peaks are shown in (a) a 2D T 1 -T 2 spectrum (left: a 3D plot and right: a 2D contour plot), (b) the corresponding 1D T 1 spectrum, and (c) the corresponding 1D T 2 spectrum. While the three peaks can all be successfully resolved in these ground truth spectra, real experiments will experience degraded spectral resolution because of nite sampling and noise. When resolution is degraded, the three peaks can still be easily resolved in (d) the 2D spectrum, though are no longer well-resolved in (e,f) either of the 1D spectra. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 44 4.1 DR-CSI simulation. (a-c) Numerical simulation setup and estimation results for three dierent multi-compartment datasets with dierent diusion-relaxation corre- lation characteristics. The rst column shows ground truth 2D diusion-relaxation correlation spectra (averaged across all voxels). The second and third columns show representative simulated diusion- and relaxation-encoded images (corresponding to b = 0 s/mm 2 and TE = 99 ms, and b = 5000 s/mm 2 and TE = 400 ms, respec- tively) that illustrate the geometry of compartmental overlap and the noise level. The fourth column shows 2D spectra estimated using DR-CSI (averaged across all voxels of the spectroscopic image), while the last three columns show spatial maps of the integrated spectral peaks from the estimated spectroscopic images. . . . . . . 55 4.2 Phantom . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 56 4.3 DR-CSI phantom results. (a-c) Example data and DR-CSI estimation results ob- tained with this phantom. The rst and second columns show representative diusion- and relaxation-encoded images (corresponding to b = 0 s/mm 2 and TE = 99 ms, and b = 5000 s/mm 2 and TE = 400 ms, respectively). The third column shows 2D spectra estimated using DR-CSI (averaged across all voxels of the spectroscopic image), while the last three columns show spatial maps of the integrated spectral peaks from the estimated spectroscopic images. . . . . . . . . . . . . . . . . . . . . 58 4.4 Spatially-varying DR-CSI spectra from the voxels corresponding to the red box drawn on the multi-compartment image in Fig. 4.3(c). Dierent peaks are color- coded by their corresponding compartments: green (S), red (R) and blue (I). . . . . 59 4.5 Estimation results from (a) conventional 1D diusion, (b) conventional 1D relax- ation, and (c) voxel-by-voxel DR-CSI (no spatial regularization). Each subgure shows (left) the estimated spectra (averaged across all voxels of the spectroscopic image) and (right) spatial maps of the integrated spectral peaks from the estimated spectroscopic image. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 61 4.6 The full set ofP = 28 diusion- and relaxation-encoded images from a single slice of a representative control mouse spinal cord. The diusion encoding orientation was parallel to the axonal tracts. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 62 List of Figures ix 4.7 Reference images and spatially-averaged DR-CSI spectra (corresponding to the par- allel diusion encoding orientation) from all of the (a) control spinal cords and (b) injured mouse spinal cords. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 63 4.8 Spatially-varying DR-CSI spectra from the boundary between WM and GM in the control cord (subject 1), corresponding to the red box drawn on the anatomical image in Fig. 4.7. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 64 4.9 Representative spatially-averaged DR-CSI spectra (corresponding to the perpendic- ular diusion encoding orientation) from (a) control and (b) injured mouse spinal cords. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 65 4.10 Spatial maps of the integrated spectral peaks from the DR-CSI reconstruction from data acquired with parallel diusion orientation. (a) The spectral regions that are integrated to generate the spatial maps (red: comp.1, blue: comp.2, green: comp.3 and yellow: comp.4). (b) The spatial maps from the spectral regions. . . . . . . . . . 66 4.11 Spatial maps of the integrated spectral peaks from the DR-CSI reconstruction from data acquired with parallel diusion orientation. (a) The spectral regions that are integrated to generate the spatial maps (red: comp.1, blue: comp.2, green: comp.3 and yellow: comp.4). (b) The spatial maps corresponding to the spectral regions. . . 67 4.12 Estimation results from conventional 1D diusion, conventional 1D relaxation, and voxel-by-voxel DR-CSI (no spatial regularization). Each subgure shows (a) the estimated spectra (averaged across all voxels of the spectroscopic image) with the spectral regions to generate the spatial maps (red: comp.1, blue: comp.2, yellow: comp.3 and green: comp.4) and (b) spatial maps from the integrated spectral regions. 68 5.1 Numerical simulation results for 2D correlation spectrum estimation. (a) Ground truth used for numerical simulation: (left) compartmental spectra f c (T 1 ;T 2 ) and (right) compartment spatial maps a c (x;y). (b) Representative simulated images. The highest SNR image (corresponding to TI = 0 ms and TE = 7.5 ms) and the lowest SNR image (corresponding toTI = 400 ms andTE = 217.5 ms) are displayed. Estimation results are shown for (c) our spatial-spectral estimation approach and (d) the conventional voxel-by-voxel estimation approach (no spatial constraint). Each subgure shows (left) the 2D spectrum obtained by spatially-averaging the 4D spec- troscopic image, (middle) spatial maps obtained by spectrally-integrating the spec- tral peak locations, and (right) a 2D spectrum obtained from a single representative voxel which contains two compartments. . . . . . . . . . . . . . . . . . . . . . . . . 75 5.2 Estimation results corresponding to simulated (a) 1D T 1 relaxometry and (b) 1D T 2 relaxometry acquisitions. Each gure shows (left) the 1D spectra obtained by spatially-averaging the 3D spectroscopic image, (middle) representative single-voxel spectra, and (right) spatial maps obtained by spectrally-integrating the spectral peak locations. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 77 5.3 Numerical simulation results as a function of . Each subgure shows (left) the 2D spectrum obtained by spatially-averaging the 4D spectroscopic image, (middle) spa- tial maps obtained by spectrally-integrating the spectral peak locations, and (right) a 2D spectrum obtained from a single representative voxel. . . . . . . . . . . . . . . 79 List of Figures x 5.4 Estimation results from the pumpkin experiment. (a) The 2D spectrum obtained by spatially-averaging the estimated 4D spectroscopic image. (b) Representative individual spectra plotted from three dierent spatial locations. In this plot, each of the spectral peaks has been numbered and color-coded (red: component 1, green: component 2 and blue: component 3). The individual spatial locations corresponding to each of the peaks are depicted using the same color-coding scheme shown in (c) a representative image acquired atTI = 0 ms andTE = 7.5 ms. (d) A high-resolution (0.5 mm 0.5 mm) reference image. (e) Spatial maps obtained by spectrally- integrating the three spectral peaks. Each of the map is color-coded based on the color-coding scheme described in (b), with the composite image shown on the right. 81 5.5 Pumpkin experiment results from (a) conventional 1D T1 relaxometry and (b) con- ventional 1D T2 relaxometry. Each gure shows (left) the estimated spectra averaged across all voxels, and (right) spatial maps of the integrated spectral peaks. . . . . . 82 5.6 The full set (P = 7 6 = 42) of T 1 -T 2 relaxation-encoded images from a single slice of the orange. The images show the magnitude of the observed signal. . . . . . . . . 83 5.7 (a): Reference image. (b): The estimated spatially-averaged 2D T 1 -T 2 relaxation spectrum. (c): The spatial maps of the integrated spectral peaks from the estimated spectroscopic image. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 84 5.8 Spatially-varyingT 1 -T 2 relaxation correlation spectra from the region-of-interest cor- responding to the green box drawn in the reference image from Figure 5.7(a). . . . . 85 5.9 A representative single-slice 4D dataset from an in vivo human brain (the axial slice from subject 1). . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 87 5.10 2DT 1 -T 2 correlation spectra estimated from an in vivo human brain (the axial slice from subject 1). (left) The 2D spectrum obtained by spatially-averaging the esti- mated 4D spectroscopic image. (middle) Representative individual spectra plotted from dierent voxels. Component 1 and component 6 are plotted from a white matter voxel, component 2 is plotted from a voxel in the putamen, component 3 is plotted from a voxel in the globus pallidus, component 4 is plotted from a voxel in gray matter, and component 5 is plotted from a voxel in the cerebral spinal uid. In this plot, each of the spectral peaks has been numbered and color-coded (red: component 1, green: component 2, cyan: component 3, blue: component 4, yellow: component 5, and magenta: component 6. This color coding scheme was also used to depict the individual voxel locations on (right) an anatomical reference image, although we do not mark the voxel for component 6 because it is the same as the voxel for component 1. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 88 5.11 Visualization of the estimated 4D spectroscopic image from an in vivo human brain (the axial slice from subject 1). Spatially-varying 2D spectra are shown from (left) the entire image (81 74 voxels), (middle) a subregion of the entire image corre- sponding to the green box (14 14 voxels), and (right) an even smaller subregion corresponding to the orange box (4 3 voxels). In each 2D spectrum, the horizontal axis corresponds to the T 2 dimension and the vertical axis corresponds to the T 1 dimension. Each spectrum is color-coded based on the six spectral peaks and the color-coding scheme described in Fig. 5.10. . . . . . . . . . . . . . . . . . . . . . . . . 89 5.12 Reference images (corresponding toTI = 0 andTE = 7.5 ms) and spatially-averaged 2D T 1 -T 2 spectra from dierent axial slices of dierent subjects. . . . . . . . . . . . . 90 List of Figures xi 5.13 Reference images (corresponding toTI = 0 andTE = 7.5 ms) and spatially-averaged T 1 -T 2 spectra obtained from our multidimensional correlation spectroscopic imaging approach from dierent coronal slices of dierent subjects. . . . . . . . . . . . . . . . 90 5.14 Spatial maps obtained from our multidimensional correlation spectroscopic imaging approach by spectrally-integrating the six spectral peaks from axial slices of dierent subjects. Each map is color-coded based on the scheme from Fig. 5.10, and the composite image is also shown on the right. . . . . . . . . . . . . . . . . . . . . . . . 91 5.15 Spatial maps obtained from our multidimensional correlation spectroscopic imaging approach by spectrally-integrating the six spectral peaks from coronal slices of dif- ferent subjects. Each map is color-coded based on the scheme from Fig. 5.10, and the composite image is also shown on the right. . . . . . . . . . . . . . . . . . . . . . 92 5.16 Results from conventional 1D T 2 relaxation. (a) The 1D spectrum obtained by spatially-averaging the 3D spectroscopic image. (b) Representative 1D spectrum from a voxel in white matter. (c) Spatial maps obtained by spectrally-integrating the three spectral peak locations. Each map is color-coded (magenta: comp.A, green: comp.B, and yellow: comp.C), and the composite image is also shown on the right. . 93 6.1 (a)The \fully sampled" diusion- and relaxation-encoded images from a control mouse spinal cord, containing 28 total images (every combination of 7 b-values and 4 TEs). (b) The sampling schemes we used with only 12 contrast encodings. . . . . . 101 6.2 Estimated 2D diusion-relaxation correlation spectra (averaged across all voxels) from representative (a) control and (b) injured spinal cords. The rst column shows representative diusion- and relaxation-encoded images. The second column shows the integrated 2D spectra from the ground truth. The right four columns show the integrated 2D spectra from the four dierent accelerated scans, corresponding to the sampling schemes shown in Fig. 6.1. . . . . . . . . . . . . . . . . . . . . . . . . . . . 102 6.3 (a) Spatially-varying DR-CSI spectra from (a) the \fully sampled" ground truth, (b) accelerated CRLB-based sampling, and (c) 4x3 Rectangular grid sampling. Spectra are shown from the transition region between injury and gray matter, corresponding to the red box drawn on the injured cord in Fig. 6.2. . . . . . . . . . . . . . . . . . . 103 6.4 Spatial maps of the integrated spectral peaks from the DR-CSI reconstruction from representative (a) control and (b) injured cords. The correspondences between spa- tial maps and spectral peaks are indicated using the color-coding scheme shown in the left column (red: comp. 1, blue: comp. 2, green: comp. 3 and yellow: comp. 4). 104 6.5 Illustration of the kind of 4D data acquisition we have previously used for RR-CSI, with 2D spatial encoding combined with 2D contrast encoding on a dense rectilinear grid. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 105 6.6 Subsampling patterns for (left) P = 36 contrast encodings and (right) P = 25 contrast encodings. Measured and unmeasured samples are marked with and symbols, respectively, and the dierence between P = 36 andP = 25 is additionally indicated with red color in the right image. . . . . . . . . . . . . . . . . . . . . . . . 106 6.7 RR-CSI subsampling results. (a): Averaged 2D RR-CSI spectra (summed over all voxels). (b): Spatial maps of the corresponding three peaks from (a) (red: component 1, green: component 2, blue: component 3). . . . . . . . . . . . . . . . . . . . . . . . 107 List of Figures xii 6.8 (a) Illustration of the ip angles used by the IR-FISP MRF sequence as a function of time. (b) Example IR-FISP MRF signal evolutions corresponding to parameter values that are typical of white matter (T1=800ms, T2=70ms) and gray matter (T1=1000ms,T2=90ms). (c) Example IR-MSE signal evolutions for reference, using the same parameter values as in (b). . . . . . . . . . . . . . . . . . . . . . . . . . . . 109 6.9 RR-CSI results from (a) the IR-MSE dataset and (b) the IR-FISP MRF dataset. In each gure, the top row shows (left) the average 2D RR-CSI spectrum (integrated across all voxels) and (middle) representative individual spectra plotted from (right) dierent spatial locations indicated on the anatomical reference image (each of the peaks and the corresponding voxel is color-coded). The bottom row shows spatial maps of the integrated spectral peaks from the estimated spectroscopic images. . . . 110 6.10 Illustration of a small region (3 x 3 voxels) from the 4D spectroscopic image ob- tained from the IR-FISP MRF dataset. The spatially-varying spectra are plotted from the transition region between putamen (comp.2) and globus pallidus (comp.3) corresponding to the orange box drawn on the reference image in Fig. 6.9. . . . . . . 111 6.11 CRLB results from the mixture of the six compartments. Each subgure shows (top) normalized CRLB values for T 1 ( p CRLB(T 1 )=T 1 ) and (bottom) normalized CRLB values for T 2 ( p CRLB(T 2 )=T 2 ). . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 112 6.12 CRLB results from three typical mixtures of three compartments. (a) The mixture of comp.1, comp.2 and comp.6. (b) The mixture of comp.2, comp.3 and comp.6. (c) The mixture of comp.3, comp.4 and comp.6. Each subgure shows (top) normalized CRLB values for T 1 ( p CRLB(T 1 )=T 1 ) and (bottom) normalized CRLB values for T 2 ( p CRLB(T 2 )=T 2 ). . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 113 List of Tables 4.1 Desired D and T 2 values and corresponding concentrations of PEG and gadobutrol for each compartment of the custom-built phantom. . . . . . . . . . . . . . . . . . . 56 xiii Abbreviations 1D 1-Dimensional 2D 2-Dimensional 4D 4-Dimensional ADC Apparent Diusion Coecient ADMM Alternating Directio Method of Multipliers CPMG Carr - Purr - Meiboom - Gill CRLB Cram er - Rao Lower Bound CSI Correlation Spectroscopic Imaging DR-COSY Diusion Relaxation - COrrelation SpectroscopY DR-CSI Diusion Relaxation - Correlation Spectroscopic Imaging DWI Diusion Weighted Imaging ILT Inverse Laplace Transform IR Inversion Recovery LT Laplace Transform MR Magnetic Resonance MRI Magnetic Resonance Imaging MSE Multi-echo Spin-Echo NMR Nuclear Magnetic Resonance RR-COSY T 1 Relaxation T 2 Relaxation - COrrelation SpectroscopY RR-CSI T 1 Relaxation T 2 Relaxation - COrrelation Spectroscopic Imaging SE Spin Echo SBS Sequential Backward Selection xiv Chapter 1 Introduction 1.1 Problem Statement Magnetic Resonance Imaging (MRI) is a powerful imaging technology that is safe and noninvasive, and can be controlled to provide images with a range of dierent contrasts. However, due to practical experimental limits, it is dicult to use MRI to directly acquire images with microscopic spatial resolution. MRI techniques that can investigate the microscopic structure of biological tissues have a range of applications in science and medicine, since many important tissue changes (due to development, pathology, aging, etc.) occur rst at a microscopic level. Luckily, it is still possible to indirectly estimate microscopic tissue features by modeling a measured MR signal as a mixture of multiple microscopic tissue compartments. Most conventional indirect estimation approaches are based on using a series of multi-contrast MR images that are encoded for one of MRI contrast mechanisms such as diusion or relaxation. However, such indirect estimation approaches often suer from ill-posedness of an inverse problem, preventing some microscopic tissue features from being distinguished. This dissertation addresses the ill-posedness of conventional indirect estimation approaches that use either diusion or relaxation contrast encoding. We deal with this problem by designing high-dimensional MRI experiments that utilize more than one contrast mechanism to infer micro- scopic tissue features. The primary goal of this dissertation is to develop a high-dimensional MRI approach that would provide substantially improved abilities of resolving microscopic sub-voxel compartments. The goal is achieved by considering and combining the following two factors: 1 Chap.1 Introduction 2 1. we consider the use of multidimensional contrast encoding that would enable separating more microscopic sub-voxel compartments that cannot be successfully resolved with traditional lower-dimensional approaches that use only either diusion or relaxation contrast encoding. 2. we consider the use of multidimensional regularized estimation that would enable high-quality spatial mapping of microscopic tissue compartments that are resolved from each image voxel. The main idea of the multidimensional regularization estimation is to leverage spatial infor- mation from neighboring voxels. Using spatial information would substantially improve the ill-posedness of the voxel-by-voxel estimation where sub-voxel compartments are estimated at each voxel independently. In this dissertation, synergistic eects on combining the two factors are demonstrated theoreti- cally and empirically. Furthermore, this dissertation investigates advanced approaches to improve experimental eciency of the new high-dimensional MRI approach guided by theoretical analysis. 1.2 Motivation MRI was rst introduced in the early 1970s [1], and has evolved since then to become one of the most important and popular biomedical imaging methods. A distinct advantage of MRI is that it is a safe and extremely exible tool that can be used to non-invasively explore the anatomical, functional, and biological properties of living tissues from a wide range of dierent perspectives. However, due to a combination of dierent practical limitations, most modern MRI methods are forced to use relatively large voxel sizes (e.g., on the order of millimeters for in-vivo human scans). This practical limit on spatial resolution prevents the use of MRI to directly investigate microscopic tissue features that may have clinical or biological signicance. For example, it would be much better for a clinical patient if pathology can be identied when it rst emerges at a microscopic scale, rather than waiting for progression to an advanced stage of disease where manifestations are clearly visible at macroscopic scales. The large imaging voxels used in MRI mean that each measurement contains a linear mix- ture of the signal contributions from the dierent microscopic compartments that coexist within the voxel. To investigate microstructure, many MRI data acquisition and parameter estimation methods have been developed to enable unmixing of the contributions from distinct microscopic Chap.1 Introduction 3 compartments. These unmixing approaches are often based on the fact that the physics of MRI data acquisition can be manipulated in a way that sensitizes the observed MRI signal to microstruc- tural tissue features that may vary between dierent microscopic compartments. For example, diusion MRI [2{4] makes the observed signal sensitive to the random microscopic diusion of water molecules, while MRI relaxometry [5{8] makes the observed signal sensitive to relaxation phenomena that are in uenced by an assortment of features from the local physical and chemical microenvironment. A variety of dierent multi-compartment unmixing methods have been proposed for relax- ation encoding alone [5{12] and diusion encoding alone [13{20]. These methods have been suc- cessful in many applications, although also have some well-known limitations. Specically, the associated unmixing problems often take the form of an \inverse Laplace transform" (ILT), i.e., an inverse problem in which the observed signal is modeled as a linear combination of exponential decay curves with unknown decay constants. The ILT is a notoriously ill-posed inverse problem [21], which can lead to unsuccessful unmixing of distinct microscopic compartments whenever the diusion or relaxation characteristics of dierent microscopic compartments are too similar to each other. Many groups have previously shown that the 2D (or higher dimensional) ILT problem is much easier than the 1D ILT [22{32], and our work on the high-dimensional MRI approach is inspired in part by earlier pioneering works in this area. Previous 2D ILT methods focused on estimating the ensemble-average microstructural features of a large volume of interest. In contrast, our primary motivation for developing the high-dimensional MRI approach is the desire to do more than identify and quantify dierent microscopic compartments, but also to produce spatial maps that reveal their distribution as a function of spatial position. In this sense, our approach is strongly inspired by the pioneering early imaging work of Professor Paul Lauterbur [1, 33], who revolutionized the MR eld by demonstrating that spatial mapping was feasible and could provide extremely valuable information that was not available from contemporary MR methods that did not attempt to produce spatial maps. Chap.1 Introduction 4 1.3 Main Results We have proposed a novel high-dimensional MRI approach to imaging microstructure. The basic idea of the method is that 1) it acquires multidimensional datasets with simultaneous and nonsep- arable encoding of multiple MRI contrast characteristics (e.g., diusion-relaxation or relaxation- relaxation) and 2) it estimates a distribution of the corresponding contrast parameters at each image voxel using a spatially-constrained estimation technique. The distribution of the contrast parameters is often called \a correlation spectrum" where multiple distinct peaks can be identi- ed. The peaks observed in each voxel are interpreted as corresponding to dierent microscopic sub-voxel compartments. Finally, spatial distributions of each peaks are provided with the form of an image (i.e., spatial maps of microscopic tissue compartments). We have evaluated the novel high-dimensional MRI method with theoretical and empirical aspects. Related publications include Refs. [34{39]. Our theoretical and empirical evaluations are summarized below: Theoretical evaluation: we have provided estimation-theoretic analysis where estimation ef- ciency is evaluated using Cram er-Rao lower bound (CRLB). The results have conrmed that our high-dimensional spatially-constrained estimation is signicantly more ecient than existing 1D/2D approaches. Empirical evaluation: the high-dimensional MRI approach can be implemented in multiple ways with a variety of choices of high-dimensional contrast encodings (e.g., diusion-T 2 , T 1 - T 2 , diusion-T 1 -T 2 , etc.). Among various options, we have implemented two dierent imaging methods. One is called diusion-relaxation correlation spectroscopic imaging (DR-CSI) that uses diusion-T 2 relaxation contrast encoding. The other is calledT 1 relaxation-T 2 relaxation correlation spectroscopic imaging (RR-CSI) that uses T 1 relaxation-T 2 relaxation contrast encoding. DR-CSI/RR-CSI methods were evaluated with a wide range of experiments in- cluding in-vivo human brain experiments that we believe the rst in-vivo demonstration with this kind of high-dimensional approaches. The results have shown that the high-dimensional method signicantly enhances resolving capabilities of microscopic sub-voxel compartments compared to lower-dimensional methods. The results have also conrmed that the high- dimensional method provides substantially improved spatial mapping of the sub-voxel com- partments compared to existing high-dimensional approaches. Chap.1 Introduction 5 We have also proposed an optimized experiment design approach that enables more ecient multidimensional data acquisitions. The method selects a small set of the most informative obser- vations guided by the estimation-theoretic analysis that we have provided. We have demonstrated the method with DR-CSI and RR-CSI experiments. The results have shown that our experiment design approach allows to substantially reduce the amount of data to acquire without a loss of information. Related publications include Refs. [40, 41]. In the context of improved eciency, we have also performed RR-CSI using a recent fast imaging technique called MR ngerprinting that enables ecient contrast encoding ofT 1 relaxation- T 2 relaxation [42]. Our theoretical and empirical evaluation on in-vivo human experiments have conrmed that a MR ngerprinting-type acquisition can be incorporated with our high-dimensional approach to make the experiments faster. A related publication includes Ref. [43]. 1.4 Organization of the Dissertation This dissertation is organized as below: Chapter 2 provides the basic knowledge of MRI that would be needed to understand following chapters. It includes the basic physics of Nuclear Magnetic Resonance (NMR) and the basic principles of MR imaging. It also introduces several basic MR imaging methods that would be helpful to understand our proposed imaging methods. Chapter 3 presents a novel high-dimensional MRI approach. It describes important features of the method including high-dimensional contrast encoding and multidimensional regularized esti- mation. It also provides estimation-theoretic analysis to justify advantages of the high-dimensional approach. Chapter 4 describes DR-CSI that makes use of diusion and T 2 relaxation contrast mecha- nisms. Evaluations of the DR-CSI method are shown with numerical simulations, phantom exper- iments and ex-vivo animal experiments. Chapter 5 describes RR-CSI that makes use ofT 1 relaxation andT 2 relaxation contrast mech- anisms. Evaluations of the RR-CSI method are shown with numerical simulations, fruit experiments and in-vivo human brain experiments. Chapter 6 illustrates advanced approaches to improve eciency of the high-dimensional ap- proach. The rst part presents an optimized experiment design method that enables sub-sampling Chap.1 Introduction 6 of high-dimensional contrast encoding. The optimized experiment design approach is then demon- strated with ex-vivo animal experiments for DR-CSI and fruit experiments for RR-CSI. The second part describes a method that uses a MR ngerprinting acquisition for RR-CSI and demonstrates the method with in-vivo human brain experiments. Finally, Chapter 7 provides conclusions of this dissertation. Chapter 2 Background This chapter provides the basic knowledge of MRI that would be needed to understand later chapters in this dissertation. The basic physics of Nuclear Magnetic Resonance (NMR) is described in Section 2.1. The basic principles of MR imaging are presented in Section 2.2. Finally, several basic MR imaging methods are illustrated in Section 2.3. 2.1 Physics of Nuclear Magnetic Resonance MRI principles are based on the nuclear magnetic resonance (NMR) phenomenon which arises from the interactions of atomic nuclei with magnetic elds. Although the NMR physics is described by quantum mechanics, MRI principles can be accurately described in a macroscopic way by classical mechanics [44]. This section will brie y describe the NMR physics with a high-level view. Extensive descriptions of the NMR physics can be found in Refs. [45{48]. 2.1.1 Nuclear Spin System and Magnetization From the quantum mechanical point of view, an atomic nucleus has an intrinsic form of angular momentum which is referred to as a nuclear spin or spin. An atomic nucleus that has a non- zero spin generates magnetic moment. The magnetic moment of a spin is described by a vector ! = x ! i + y ! j + z ! k and is proportional to the spin angular momentum ! J as follows: ! = ! J; (2.1) 7 Chap.2 Background 8 where = 2 is the gyromagnetic ratio which is a physical constant that is unique to an atomic nucleus. For example, in MRI, Hydrogen ( 1 H) is the most commonly studied atomic nucleus be- cause it is abundant in biological samples and the gyromagnetic ratio for ( 1 H) is = 42:58 MHz/T. Based on quantum mechanics, the magnitude of the magnetic moment is dened by: = ~ p I(I+1); (2.2) where I is the spin quantum number of a nucleus and ~ is the Planck's constant. The spin quantum number for 1 H is I = 1 2 . In the absence of an external magnetic eld, the direction of the magnetic moment is random due to its thermal random motion. When an atomic nucleus is exposed to an external magnetic eld ! B 0 = B 0 ! k which is applied along the z-direction, the z-component of the magnetic moment is dened by: z = m I ~; (2.3) wherem I is the magnetic quantum number such thatm I =I; I+1;:::; I. For 1 H, the magnetic quantum number can have two values m I = 1 2 ; 1 2 , which leads to two possible energy eigenstates (\parallel" and \antiparallel"). While the spins are not required to fully be in one of the two states, the spins can exist in a superposition of the two eigenstates. Figure 2.1: Nuclear spins and magnetic moments. (a) A nuclear spin with nonzero magnetic moment . (b) Nuclear spins that are randomly oriented in the absence of an external eld. (c) Nuclear spins that are aligned in the presence of an external eld. To understand MRI principles, now we account for a nuclear spin system that comprises an ensemble of nuclei of the same type per unit volume, rather than an individual spin. We focus Chap.2 Background 9 on describing the behavior of a nuclear spin system from the classical point of view. The behavior of a nuclear spin system is described by the net magnetization vector (i.e., the bulk magnetization vector) ! M = M x ! i +M y ! j +M z ! k and the net magnetization vector represents a summation of magnetic moments that exist within an unit volume: ! M = N X n=1 ! n ; (2.4) whereN is the number of nuclear spins per an unit volume and ! n indicates the magnetic moment of the n-th nuclear spin in the volume. In the absence of an external magnetic eld, the net magnetization is zero since magnetic moments are randomly oriented (i.e., maximum entropy state). If the strong external magnetic eld ! B 0 = B 0 ! k is applied along the z-direction, spins are aligned in two opposite orientations pointing-up and pointing-down, and spins in two opposite orientations belong to dierent energy states interacting with the external magnetic elds. The energy dierence between two states is nonzero which is referred to as the Zeeman splitting phenomenon in quantum mechanics and it is given by: E =E # (pointing-down)E " (pointing-up) = 1 2 ~B 0 ( 1 2 ) ~B 0 = ~B 0 : (2.5) In the aspect of the net magnetization, this phenomenon can be macroscopically viewed as polarization and the energy state of the spin system is in a thermally equilibrium state. In the equilibrium state, the direction of the net magnetization vector is the same direction of the external magnetic eld in the z-axis and the magnitude of the net magnetization vector is given by: M 0 z =j ! Mj = 2 ~ 2 B 0 N 4KT ; (2.6) where~ is the Planck's constant, K is the Boltzmann's constant andT is an absolute temperature of the spin system. It should be noted that Eqs. (2.5) and (2.6) describe only the spin system that has a spin quantum numberI = 1 2 such as 1 H. For a spin system withI, in general, the magnitude of the net magnetization can be given by (the detailed description of the equation can be found in Ref [49]): M 0 = N 2 ~ 2 I (I + 1)B 0 3KT : (2.7) Chap.2 Background 10 2.1.2 RF Excitation Precession and Larmor Frequency In the presence of the external magnetic eld ! B 0 , the magnetization vector ! M 0 in thermal equilibrium is aligned in the same direction of ! M 0 . However, if the direction of the magnetization ! M is not aligned to the direction of the applied eld ! B , precession of the magnetization occurs about the applied eld ! B as described in Fig. 2.2. From classical mechanics, the torque is applied to a spin ! in the presence of ! B and is equivalent to the change rate of angular momentum as described: d ! J dt = ! ! B; (2.8) where indicates an operator for the cross product. According to Eq. (2.1), the change rate of the spin's magnetic momentum is given by: d ! dt = ! ! B: (2.9) Therefore, the change rate of the net magnetization becomes: d ! M dt = ! M ! B: (2.10) The precession frequency is given by so-called the Larmor frequency ! that represents the resonance frequency of a spin system. The Larmor frequency is proportional to the magnitude of the applied magnetic eld B: ! = B rad/s: (2.11) Chap.2 Background 11 Figure 2.2: Precession of the magnetization about the applied magnetic eld. RF pulse To obtain a MR signal, a radio frequency (RF) magnetic eld ! B 1 (t) is applied to the equi- librium state of the spin system in the presence of the main magnetic eld ! B 0 . Distinct from the static ! B 0 eld, the RF magnetic eld ! B 1 (t) is applied during a short time (e.g., a few milliseconds on a human imaging system) with small magnitude, and it is an oscillating magnetic eld that is tuned to the resonance frequency of the spin system (i.e., ! B 1 (t) is rotating in the same way that the net magnetization is precessing). In the quantum mechanical perspectives, applying the RF magnetic eld induces energy exchange between the two energy states described in Eq. (2.5) if the radiation energy due to the RF magnetic eld is equal to the energy dierence between the two states as follows: ~! rf = E: (2.12) According to Eq. (2.5), the resonance frequency of the spin system is given by: ! rf =! 0 = B 0 : (2.13) Typically, the RF magnetic eld is applied to the xy plane (i.e., the transverse plane that is perpendicular to the direction of the main magnetic eld) by taking the following form: ! B 1 (t) = 2B e 1 (t)cos(! rf t + ) ! i ; (2.14) Chap.2 Background 12 whereB e 1 (t) is a pulse envelope function and is an initial phase angle. This led is mathematically composed of two circularly polarized elds: ! B 1 (t) =B e 1 (t) h cos(! rf t + ) ! i sin(! rf t + ) ! i i +B e 1 (t) h cos(! rf t + ) ! i + sin(! rf t + ) ! i i : (2.15) In the above equation, since the second term rotates counter-clockwise which is the opposite di- rection of the precessing spins, it is negligible by Eq. (2.13). Hence, the RF magnetic eld can be approximated by the following complex form: ! B 1 (t) =B e 1 (t)e i(! rf t+ ) : (2.16) According to Eq. (2.16), the RF magnetic eld is characterized with three parameters: the envelope function B e 1 , the oscillating frequency ! rf and the initial phase angle . The RF magnetic eld is generated by so-called an RF pulse via RF transmitter coils. More detailed descriptions about an RF pulse can be found Ref. [48]. Excitation and Flip Angle Once the RF magnetic eld is applied to the perpendicular direction of the main magnetic eld, the net magnetization precesses about the z-axis, simultaneously precessing about the x-axis as shown Fig. 2.3 (a). Now we can view this behavior in the Larmor-rotating frame which is a coordinate system whose transverse plane is rotating at the Larmor frequency ! 0 as follows: 8 > > > < > > > : ! i 0 :=cos(! 0 t) ! i sin(! 0 t) ! j ! j 0 :=sin(! 0 t) ! i +cos(! 0 t) ! j ! k 0 := ! k: (2.17) As shown in Fig. 2.3 (b), in the rotating frame, the net magnetization is tipped away from z-axis, precessing about x-axis at the angular frequency ! 1 = B 1 . We call this precession excitation. In the rotating frame, the angle between z 0 -axis and the net magnetization is referred to as a ip angle whose value is determined by an RF pulse as follows: = Z p 0 ! 1 (t)dt = Z p 0 B e 1 (t)dt; (2.18) Chap.2 Background 13 where p indicates the duration of an RF pulse. In the case that a rectangular pulse with the magnitude B 1 is applied, =! 1 t = B 1 p : (2.19) To obtain MR signals, a 90 pulse is typically used for the excitation. As a result, the transverse component of the net magnetization is measurable through RF coils and can be the source of MR signals. Figure 2.3: Behavior of the net magnetization in (a) the laboratory frame (i.e., the static frame) and (b) the Larmor-rotating frame. In the rotating frame, the magnetization is tipped into the transverse plane by the application of the RF magnetic eld. 2.1.3 Relaxation Once a spin system is excited from its thermal equilibrium state by the RF magnetic eld, the energy state of the spin system is perturbed and it will return to the thermal equilibrium state according to thermodynamics. This returning process involves a precession of the magnetization ! M about the static eld ! B 0 , which is referred to as relaxation. Relaxation is characterized by two magnetization components; the longitudinal magnetizationM z in thez-axis and the transverse magnetization M xy in the xy plane. Chap.2 Background 14 Figure 2.4: Relaxation curves in (a) the longitudinal direction and in (b) the transverse plane. As can be seen in Fig. 2.4 (a), the longitudinal component of the excited magnetization is exponentially recovered to the equilibrium condition with a time constant T 1 , which is referred to as T 1 relaxation. The behavior of the longitudinal magnetization with time is described below: dM z dt = M z M 0 z T 1 ; (2.20) where M 0 z is the longitudinal magnetization at thermal equilibrium. The solution of the above is followed by: M z =M 0 z + M z (0)M 0 z e t=T 1 : (2.21) In MRI of Hydrogen ( 1 H), T 1 relaxation is associated with an energy exchange between water protons and its surrounding lattice (e.g., lipids, proteins and macromolecules). In this context, T 1 is also called spin-lattice relaxation time. As can be seen in Fig. 2.4 (b), the transverse component of the excited magnetization expo- nentially decays to zero with a decay time constant T 2 , which is referred to as T 2 relaxation. The behavior of the transverse magnetization is described by: dM xy dt = M xy T 2 ; (2.22) and the solution is followed by: M xy =M xy (0)e t=T 2 : (2.23) SinceT 2 relaxation is associated with an energy exchange between spins,T 2 is also called spin-spin Chap.2 Background 15 relaxation time. It should be noted that the solutions in Eqs. (2.20) and (2.22) are derived in the rotating frame. The relaxation time constants T 1 and T 2 are intrinsic biophysical tissue properties that can determine MR signal strength. By manipulating RF pulses, we can generate image contrasts between dierent tissues that have unique T 1 and T 2 characteristics. Generation T 1 -weighted and T 2 -weighted images will be discussed in Section 2.3. 2.1.4 The Bloch Equation In summary, in the presence of the static eld ! B 0 , the net magnetization vector in thermal equi- librium has perturbed by applying the RF magnetic eld ! B 1 (t) during a short time period. As a result, the net magnetization tipped to the transverse plane will start relaxation. This phenomenon can be described by the Bloch equation that comprises both the precession term (Eq. (2.10)) and the relaxation terms (Eqs. (2.20) and (2.22)): d ! M dt = ! M ! B M x ! i +M y ! j T 2 M z M 0 z ! k T 1 ; (2.24) where M 0 z is the longitudinal magnetization component in thermal equilibrium; ! B is an applied magnetic eld; and T 1 and T 2 are relaxation time constants. To solve the Bloch equation, we now account for the static magnetic eld ! B 0 =B 0 ! k as the simplest case. First, ignoring relaxation eects, the Bloch equation can be written by the precessing term in a matrix form as follows: 2 6 6 6 4 dMx dt dMy dt dMz dt 3 7 7 7 5 = 2 6 6 6 4 0 B 0 0 B 0 0 0 0 0 0 3 7 7 7 5 2 6 6 6 4 M x M y M z 3 7 7 7 5 : (2.25) The solution of the above equation is given by: M(t) = 2 6 6 6 4 M x (t) M y (t) M z (t) 3 7 7 7 5 = 2 6 6 6 4 cos! 0 t sin! 0 t 0 sin! 0 t cos! 0 t 0 0 0 1 3 7 7 7 5 M(0); (2.26) Chap.2 Background 16 where M 0 is an initial magnetization vector and ! 0 = B 0 is the Larmor frequency. The solution describes precessing of the magnetization vector aboutz-axis at the Larmor frequency. The solution is simply written by the rotation matrix R z () that performs a rotation aboutz-axis with an angle : M(t) = R z (! 0 t)M(0): (2.27) If we consider relaxation eects, the Bloch equation can be written in a matrix form as below: dM dt = 2 6 6 6 4 1=T 2 B 0 0 B 0 1=T 2 0 0 0 1=T 1 3 7 7 7 5 M + 2 6 6 6 4 0 0 M 0 z =T 1 3 7 7 7 5 ; (2.28) where M 0 z is the longitudinal magnetization component in thermal equilibrium. The solution of the above equation becomes: M(t) = 2 6 6 6 4 e t=T 2 0 0 0 e t=T 2 0 0 0 e t=T 1 3 7 7 7 5 R z (! 0 t)M(0) + 2 6 6 6 4 0 0 M 0 z (1e t=T 1 ) 3 7 7 7 5 : (2.29) We can also represent the solution as the transverse component M xy =M x +iM y and the longitu- dinal component M z : 8 < : M xy (t) =M(0)e t=T 2 e i! 0 t M z (t) =M 0 z + (M(0)M 0 z )e t=T 1 : (2.30) 2.2 Basic Principles of MR Imaging In previous sections, we discussed about the basics of nuclear magnetic resonance phenomenon focusing on how MR signals are generated. In this section, we will brie y describe how MR images can be generated from MR signals. Chap.2 Background 17 2.2.1 The Bloch Equation with Gradient Fields To generate MR signals, two types of magnetic elds are essential: the main static eld ! B 0 which is used for magnetization and the oscillating RF magnetic eld ! B 1 (t) which is used for excitation. To generate MR images, an additional magnetic eld B( ! r;t) called Gradient eld is required for spatial localization of the MR signals. The gradient eld B( ! r;t) is a spatially-varying magnetic eld that is typically constructed by: B( ! r;t) =G x (t)x +G y (t)y +G z (t)z = ! G(t) ! r; (2.31) where ! G(t) = G x (t) ! i +G y (t) ! j +G z (t) ! k is a vector that represents linear gradient elds and ! r = x ! i +y ! j +z ! k is a spatial coordinate vector. The gradient eld B( ! r;t) is applied with the main static eld ! B 0 , hence the total magnetic eld becomes: ! B = (B 0 + B( ! r;t)) ! k (2.32) = (B 0 + ! G(t) ! r ) ! k: (2.33) As previously mentioned, the transverse component of the net magnetization is the source of MR signals. Using the total magnetic eld in Eq. 2.33, the transverse component of the Bloch equation solution is given by: M xy ( ! r;t) =M 0 ( ! r )e t=T 2 ( ! r ) e i! 0 t e (i R t 0 ! G() ! rd) ; (2.34) whereM 0 ( ! r ) is an initial condition (i.e., M( ! r; 0)). In the expression, we assume that T 2 ( ! r ) can be spatially-varying corresponding to the local tissue sample that will be measured. 2.2.2 Signal Detection We now describe how the transverse magnetization is measured. The precessing of the magnetiza- tion causes changes of the magnetic ux in a receiver coil that surrounds an image object, and the electromagnetic force (emf) is induced according to Faraday's law of induction [50]: = @ @t : (2.35) Chap.2 Background 18 According to the principle of reciprocity [51], the incremental EMF is given by: d = @ @t h ! B 1 ( ! r ) ! M( ! r;t) i dV; (2.36) where ! B 1 ( ! r ) is the RF magnetic eld at the spatial location ! r and ! M( ! r;t) is the net magne- tization vector. The MR signal that we measure is the total EMF signal s r (t) that is derived by integrating the incremental EMF over the volumes: s r (t) = Z V d = Z V @ @t h ! B 1 ( ! r ) ! M( ! r;t) i dV: (2.37) If we assume that the RF magnetic eld is applied to the transverse plane and is homogeneous over the unit volume, then the EMF signal is given by: s r (t) = Z V B 1xy @ @t M xy ( ! r;t)dV; (2.38) where B 1xy is the transverse component of the homogeneous RF magnetic eld and M xy ( ! r;t) is the transverse magnetization component that is derived from the Bloch equation in Eq. (2.34). Using Eq. (2.34) and ignoring T 2 decay eects, the EMF signal becomes: s r (t) = Z V B 1xy M 0 ( ! r ) h i(! 0 + ! G(t) ! r )e i! 0 t e (i R t 0 ! G() ! rd) i dV: (2.39) Since the main magnetic eld is much bigger than the gradient eld in most cases, ! 0 ! G(t) ! r , s r (t) =i! 0 B 1xy Z V M 0 ( ! r )e i! 0 t e (i R t 0 ! G() ! rd) dV: (2.40) Ignoring the constant term i! 0 B 1xy , the nal form of the signal equation is: s r (t) = Z Z Z M 0 (x;y;z)e i! 0 t e (i R t 0 ! G() ! rd) dxdydz: (2.41) Chap.2 Background 19 A complex-valued received signal s r (t) can be demodulated to a complex-valued baseband signal s(t) that is desired to be measured: s(t) =s r (t)e i! 0 t (2.42) = Z Z Z M 0 (x;y;z)e (i R t 0 ! G() ! rd) dxdydz; (2.43) which can be described by: s(t) =(t)e (t) ; (2.44) where(t) indicates the amplitude of the baseband signal and(t) indicates the phase modulation of the baseband signal. While the expression is analyzed by complex-valued signals, in practice, the real-value of s r (t) is measured and demodulated thorough a process known as quadrature phase- sensitive detection. As can be seen in Fig. 2.5, the input RF signal goes through two channels: in the rst channel, the input is multiplied by cos(! 0 t) and low-pass ltered and in the second channel, the input is multiplied by by sin(! 0 t) and low-pass ltered. The output is composed of I(t) = (t)cos(! 0 t) (in-phase signal) and Q(t) =(t)sin(! 0 t) (quadrature signal), followed by digitization. Hence, the baseband signal is obtained by: s(t) =I(t) +iQ(t) =(t)e (t) : (2.45) Figure 2.5: Block diagram of the quadrature phase-sensitive detection. Chap.2 Background 20 2.2.3 Fourier Interpretation In MR imaging, the measured baseband signal s(t) can be associated with the Fourier transform, which is signicant to form MR images from the measured MR signals. For simplicity, we assume that excitation is selectively performed over a 2D plane that is centered at z = z 0 and whose thickness is z. The 2D function that will be imaged is then dened as follows: m(x;y) := Z z 0 +z=2 z 0 z=2 M 0 (x;y;z)dz; (2.46) and thereby the signal measured from the 2D plane becomes: s(t) = Z Z m(x;y)e i2[kx(t)x+ky (t)y] dxdy; (2.47) where k x (t) and k y (t) are time integrals of the gradient waveforms dened by: k x (t) = 2 Z t 0 G x ()d (2.48) k y (t) = 2 Z t 0 G y ()d: (2.49) In the expression, the signal can be viewed as the 2D Fourier transform of m(x;y) with spatial- frequency variables k x (t) and k y (t), and spatial-frequency domain of the 2D Fourier transform is called k-space [52]. In MR imaging, the Fourier relationship is important because it allows to recoverm(x;y) from the signals(t) that is measured/sampled under certain conditions. More detail about k-space sampling requirements and image reconstruction can be found in Refs. [13, 46, 48]. 2.3 Basic Imaging Methods 2.3.1 Free Induction Decay (FID) A free induction decay (FID) signal is the transient response of a spin system after excitation. `Free' represents that the signal is produced by the free precession of the net magnetization vector about the static ! B 0 eld, `Induction' stands for that the signal is generated as a result of the Faraday's law of electromagnetic induction and `Decay' indicates that the signal exponentially decreases with time. Chap.2 Background 21 To mathematically describe the FID signal, we rst introduce a spin isochromat that represents a group of spins that are precessing at the same Larmor frequency experiencing the same local eld environment. Assuming a heterogeneous spin system where multiple isochromats exist, the signal in Eq. 2.43 can be expressed in the precessional frequency frame: s(t) = Z 1 1 (!)e t=T 2 (!) e i!t d!: (2.50) where (!) is a spin spectral density function that characterizes the frequency distribution of a heterogeneous spin system. Note that in the expression, we account for T 2 decay eects although we have omitted the T 2 decay term in Eq. (2.43). The spin spectral density function can be often dened with a Lorentzian distribution for a heterogeneous spin system with multiple isochromats: (!) =M 0 z ( B 0 ) 2 ( B 0 ) 2 + (!! 0 ) 2 ; (2.51) where is the gyromagnetic ratio and ! 0 is the Larmor frequency dened by the Larmor equation ! 0 = B 0 . Using Eq. (2.51), the FID signal is given by: s(t) = Z 1 1 M 0 z ( B 0 ) 2 ( B 0 ) 2 + (!! 0 ) 2 e t=T 2 (!) e i!t d! (2.52) =M 0 z B 0 e B 0 t e t=T 2 e i! 0 t (2.53) =M 0 z B 0 e t=T 2 e i! 0 t : (2.54) In the expression, we introduce a new relaxation time constant T 2 which is often used in the MR literature: 1 T 2 = 1 T 2 + B 0 : (2.55) As shown in Fig. 2.6, the envelop of a FID signal decays with a time constant T 2 . It should be noted that this type of the FID signal is valid under the assumption that the distribution of a spin system follows Lorentzian distribution. It should be also noted that, while the FID signal equation in Eq.(2.54) is derived from a 90 excitation, the FID signal equation in the case of applying a ip Chap.2 Background 22 angle becomes: s(t) =M 0 z B 0 sin e t=T 2 e i! 0 t : (2.56) Figure2.6: A simulated FID signal from Eq. 2.54. The signal envelop decays with a time constant T 2 . Figure 2.7 graphically illustrates behaviors of the vectors of multiple isochromats. When a spin system is placed in a inhomogeneous magnetic eldB 0 + B, multiple isochromats can exist. As can be seen, isochromats are precessing at dierent angular frequencies, which causes phase dispersions of the vectors of multiple isochromats. As a result, the transverse component of the net magnetization vector decays with a time constant T 2 and this phenomenon is referred to as dephasing or T 2 decay. Figure 2.7: Graphical description of (a) 90 excitation and (b) dephasing of isochromats in a spin system after the excitation in the rotating frame. Chap.2 Background 23 2.3.2 Spin-Echo Imaging Spin Echo Echo is another form of MR signals that can be induced from a FID signal. An echo signal was discovered by Erwin L. Hahn in 1950 [53] and it can be generated by two RF pulses. In MR literature, 90 and 180 are commonly used to generate an echo signal and the echo signal generated from the two pulses is often called spin echo. Fig. 2.8 illustrates behaviors of isochromats as a result of the application of 90 and 180 pulses. After the application of a 90 excitation pulse along the x 0 -axis, isochromats are dephasing causing that the net magnetization disappears. At certain time, a 180 pulse is applied along the y 0 -axis and then all the isocromats are rotated by 180 (i.e., ipped) about the y 0 -axis. As the ipped isochromats continue to precess, the phases of the isochromats become coherent, which is referred to as rephasing. As a result of rephasing, the transverse component of the net magnetization is recreated. This recreated signal is a spin echo signal. Figure 2.8: Graphical description of spin echo generation. Chap.2 Background 24 Fig. 2.9 shows a sequence diagram to generate a spin echo signal. As can be seen, the FID signal is generated right after the application of a 90 pulse, and then a spin echo signal is generated after the application of a 180 pulse. The time duration between a 90 pulse and the center of a spin echo signal is called echo time denoted by TE. A 180 pulse is applied at t = TE=2. Since a 180 pulse is used to regenerate the signal, it is often called a refocusing pulse or rephasing pulse. According to the NMR physics, the amplitude at the center of a spin echo signal is given by: M(t) =M 0 z e TE=T 2 : (2.57) This expression shows that the amplitude depends onT 2 and TE, which is an important feature to generateT 2 -weighted image contrast. It should be also noted that an echo signal can be generated by a combination of arbitrary RF pulses 1 and 2 , generating the echo signal whose amplitude is dened by: M(t) =M 0 z sin 1 sin 2 2 e TE=T 2 : (2.58) Figure 2.9: Generation of a spin echo signal. Chap.2 Background 25 Spin-Echo Imaging We now describe how MR images are generated using spin echoes. As brie y discussed in Sec. 2.2.3, a 2D MR image m(x;y) can be reconstructed by the inverse 2D Fourier transform of k-space data. In spin-echo imaging, k-space data is lled with multiple spin echoes that are sequentially acquired with a series of 90 excitation pulses. For example, if we want to acquire 2D k-space data with a size of N x N y in spin-echo imaging, N y echoes should be acquired to ll k-space and each echo should be discretized with N x samples. To obtain multiple spin echoes, 90 excitation pulses are equally spaced with time interval TR which is called repetition time. Fig. 2.10 shows a sequence diagram of the multiple 90 excitation pulses and the longitudinal magnetization corresponding to the pulse sequence. The sequence is often called a saturation recovery sequence because spins are in a saturated state in this sequence. Under the assumption that TR T 2 , the transverse magnetization disappears at the end of each TR, meaning that the recovery of the longitudinal magnetization at the time of the next pulse takes place under an initial condition with no magnetization from the previous pulse. The longitudinal magnetization in the saturated state before the n-th excitation pulse is derivate using the longitudinal magnetization equation in Eq. (2.30) and M z ((n 1)TR + ) = 0: M z (nTR ) =M 0 z (1e TR=T 1 ): (2.59) Figure 2.10: The saturation-recovery pulse sequence and the time evolution of the longitudinal magnetization. Spin-echo imaging is commonly performed with the saturation recovery sequence. Fig. 2.11 represents a saturation-recovery spin-echo imaging sequence. Within each TR, a refocusing pulse is Chap.2 Background 26 applied and then a spin-echo signal is acquired. As previously mentioned in Sec. 2.2, gradient elds G x (t);G y (t) and G z (t) are applied along with RF pulses to spatial localization of the spin-echo signal 1 . Assuming TE TR, the amplitude of the spin echo from the sequence is given by: M(t) =M 0 z (1e TR=T 1 )e TE=T 2 : (2.60) which determines the image intensity acquired by the saturation-recovery spin-echo imaging se- quence. While this image intensity contains aT 1 -weighting,T 2 -weighting and spin-density-weighting simultaneously, we can also emphasize one of contrast mechanisms by controlling the sequence pa- rameters TR and TE. For example, if we use a long TR, the term (1e TR=T 1 ) can be ignored and thereby a T 2 -weighted image with the intensity M 0 z e TE=T 2 can be obtained. Similarly, a T 1 -weighted image can be acquired with a short TE and an appropriate TR and a spin-density- weighting image can be acquired with a short TE and a long TR. Figure 2.11: Spin-echo imaging sequence. 1 In Chapter 2, we didn't deal with k-space sampling and image reconstruction in detail. More detailed descriptions about those concepts can be found Refs. [45, 46, 48] Chap.2 Background 27 The CPMG echo train When a series of 180 pulses follow a 90 pulse, a train of spin echoes can be generated as shown in Fig. 2.12. The sequence was developed in 1950s [54] by Carr and Purcell, and it is called the CPMG (Carr-Purcell-Meiboom-Gill) sequence. An advantage of this method is that multiple spin-echo signals at dierent echo times are acquired within a TR, which reduce the acquisition time. This method has been widely used in fast imaging when multi-echo spin-echo images need to be acquired specically for measuring T 2 . As an example, Fig. 2.13 shows representative T 2 - weighted brain images acquired using a CPMG-based imaging sequence with TR=5000 ms and 15 TEs ranging from 7.5 ms to 217.5 ms in 15 ms increments. Figure 2.12: The CPMG pulse sequence. Chap.2 Background 28 Figure 2.13: T 2 -weighted brain images acquired using a CPMG-based imaging sequence with TR=5000 ms and 15 TEs ranging from 7.5 ms to 217.5 ms in 15 ms increments. Four representative images are displayed and each image is acquired with (a) TE = 22.5 ms, (b) TE = 67.5 ms, (c) TE = 112.5 ms and (d) TE = 217.5 ms. 2.3.3 Inversion-Recovery Imaging Inversion-Recovery (IR) imaging is a gold standard method to measure T 1 [55] and the method was developed in 1940s by Drain [56] and Hahn [57]. As shown in Fig. 2.14, the basic IR pulse sequence consists of a preparatory inversion pulse (180 pulse) and an excitation pulse (90 pulse). In the sequence, the time interval between the inversion pulse and the excitation pulse is called Inversion Time denoted by TI. Chap.2 Background 29 Figure 2.14: The inversion-recovery pulse sequence and the time evolution of the longitudinal magnetization. As in the saturation-recovery sequence, the longitudinal magnetization in the inversion- recovery sequence reaches a steady state before the second 90 pulse. The time evolution of the longitudinal magnetization is followed by: 1. After the rst 180 at t=0, M z (0 + ) =M 0 z : (2.61) 2. Recall the longitudinal magnetization derived from the Bloch equation in Eq. (2.30), M z (t) =M 0 z + (M z (0)M 0 z )e t=T 1 ;; (2.62) whereM 0 z is the longitudinal magnetization in equilibrium state and M z (0) is the initial lon- gitudinal magnetization . Then, using Eqs. (2.61) and (2.62), the longitudinal magnetization at t = TI becomes: M z (TI ) =M 0 z + (M z (0 + )M 0 z )e TI=T 1 (2.63) =M 0 z (1 2e TI=T 1 ): (2.64) Chap.2 Background 30 3. After the rst 90 pulse, the longitudinal magnetization is ipped into he transverse plane and the longitudinal magnetization recovers to: M z ((TRTI) ) =M 0 z (1e (TRTI)=T 1 ): (2.65) 4. After the second 180 at t = TR-TI, M z ((TRTI) + ) =M 0 z (1e (TRTI)=T 1 ): (2.66) 5. Eventually, the longitudinal magnetization recovers to the following value before the second 90 , M z (TI ) =M 0 z + (M z ((TRTI) + )M 0 z )e TI=T 1 (2.67) =M 0 z (1e TI=T 1 )M 0 z (1e (TRTI)=T 1 )e TI=T 1 (2.68) =M 0 z (1 2e TI=T 1 +e TR=T 1 ): (2.69) Hence, when the inversion-recovery sequence is combined with a spin-echo imaging as shown in Fig. 2.15, the amplitude of the spin echo is given by: M(t) =M 0 z (1 2e TI=T 1 +e TR=T 1 )e TE=T 2 : (2.70) This signal expression indicates that this sequence can provide T 1 -, T 2 - and spin-density-weighted image contrast as the saturation-recovery spin-echo does. For example, with a short TE and a long TR, the signal expression can be represented by: M(t) =M 0 z (1 2e TI=T 1 +e TR=T 1 ); (2.71) where the T 2 term can be omitted, resulting in a T 1 -weighted signal. Fig. 2.16 shows T 1 -weighted images using the inversion-recovery spin-echo imaging sequence with the parameters TR = 5000 ms, TE = 7.5 ms and 4 TIs = 100, 200, 400 and 700 ms. Chap.2 Background 31 It has been known that the inversion-recovery-based imaging produces greater contrast for T 1 , so the inversion-recovery spin-echo imaging is more commonly used in T 1 -weighted imaging. In contrast, the saturation-recovery spin-echo imaging is more commonly used inT 2 -weighted imaging. Figure 2.15: Inversion-recovery (IR) spin-echo (SE) imaging sequence. Chap.2 Background 32 Figure 2.16: T 1 -weighted brain images acquired using a IR-SE imaging sequence with TR=5000 ms and TE = 7.5 ms. Each image is acquired with (a) TI = 100 ms, (b) TI = 200 ms, (c) TI = 400 ms and (d) TI = 700 ms. 2.3.4 Diusion-weighted Imaging Water diusion is a random or Brownian motion of water molecules. In the presence of a magnetic eld, water diusion causes dephasing of the transverse magnetization, which leads to decay of the MR signal. The aect of diusion to the MR signal can be understood by simply considering the generation of a spin echo described in Fig. 2.9. In the absence of diusion or a random motion of water molecule, the dephased spins before the application of a refocusing (180 ) pulse will be exactly rephased after the application of the refocusing pulse. However, in the presence of diusion, Chap.2 Background 33 there has been a random drift of the spins between dephasing and rephasing. The refocusing of the spins is then partially done, resulting in the loss of the signal amplitude at t = TE. The amount of diusion attenuation can be measured by adding pulsed eld gradients (PFG). The PFG-based method was developed in 1960s by Stesjkal and Tanner [58] and Fig. 2.17 shows the Stesjkal-Tanner PFG diusion weighted spin echo sequence where a pair of the same two diusion gradients is applied before and after the 180 refocusing pulse. The amplitude g and duration of the gradients, and the time interval between the gradients determine the amplitude of the spin echo (i.e., the amount of the diusion attenuation) and these parameters can be represented by a diusion sequence parameter b-value: b = 2 2 1 3 g 2 : (2.72) While we illustrate just a pair of the diusion gradients in the sequence for simplicity, multiple pairs of the diusion gradients can be applied along the x-, y- and z-axes depending on the direction along which we want to measure the diusion attenuation 2 . Figure 2.17: Stesjkal-Tanner PFG diusion-weighted spin-echo sequence. 2 More detailed descriptions of designing diusion gradients can be found in Ref. [48]. Chap.2 Background 34 Along with the sequence parameterb-value, the amplitude of the spin echo is also determined by the diusion properties of tissue samples. The diusion property can be represented by a diusion coecient D which is dened by the Einstein equation 3 as follows: <R 2 >= 2Dt; (2.73) where R is the net vector of the displacement by a random motion of a molecule in time t and < R 2 > is the 1D root-mean-squared displacement. The unit of D is distance squared per time (e.g., mm 2 /sec). The amplitude of the spin echo produced by the diusion-weighted spin-echo sequence is give by: M(TE;b) =Ce bD e TE=T 2 ; (2.74) where C is a constant term including an initial magnetization. The diusion coecient can be calculated by two diusion-weighted signals obtained with two dierent b-values as below: ln M(TE;b 1 ) M(TE;b 2 ) =ln " Ce b 1 D e TE=T 2 Ce b 2 D e TE=T 2 # = (b 2 b 1 ): (2.75) The calculated diusion coecient is often called an Apparent Diusion Coecient (ADC) because it can be in uenced by tissue perfusion, partial volume eects and other experimental errors [60]. In addition toT 1 andT 2 , a diusion coecient is also an important MRI contrast parameter to probe microstructure, because the value of the diusion coecient can be aected by surrounding microenvironment (i.e., macromolecules, organelles, cell membrane and other cellular and subcel- lular structures) [60]. 2.3.5 High-dimensional Contrast Encoding In previous sections, we have described spin-echo imaging, inversion-recovery imaging and diusion- weighted imaging methods that are commonly used to obtainT 1 -weighted,T 2 -weighted and diusion- weighted MR signals, respectively. However, while either of the methods can be used for encoding 3 More detailed descriptions of the Einstein equation can be found in Ref. [59]. Chap.2 Background 35 one contrast mechanism (i.e., 1D contrast encoding), we can extend the methods for simultane- ously encoding more than one contrast mechanisms (i.e., high-dimensional contrast encoding). As this dissertation mainly deals with MR images that are acquired with high-dimensional contrast encoding, we introduce two imaging sequences that would be helpful to understand later chapters. Figure 2.18 shows the inversion-recovery multi-echo spin-Echo (IR-MSE) imaging sequence that can be used to acquire nonseparable and simultaneous encoding of T 1 and T 2 contrast infor- mation. Assuming TR T 1 , the image intensity of this sequence is given by: M(TI;TE) =M 0 (1 2e TI=T 1 )e TE=T 2 : (2.76) In the expression, the image intensity is controlled by the sequence parameters TI for T 1 contrast encoding andTE forT 2 contrast encoding. Fig. 5.9 shows a human brain dataset that was acquired with the IR-MSE imaging sequence with a combination of 7 TIs (0, 100, 200, 400, 700, 1000 and 2000 ms) and 15 TEs (from 7.5 ms to 217.5 ms in 15 ms increments) for a total of 105 images. Figure 2.19 shows the diusion-weighted spin-echo imaging sequence that can be used to acquire nonseparable and simultaneous encoding of diusion and T 2 contrast information. The image intensity of this sequence is given by: M(b;TE) =M 0 e bD e TE=T 2 : (2.77) In the expression, the image intensity is controlled by the sequence parameters b-value for diusion contrast encoding and TE for T 2 contrast encoding. Fig. 4.6 shows a rat spinal cord dataset that was acquired with the diusion-weighted spin-echo imaging sequence with a combination of 7 b- values 0, 500, 1000, 2000, 3000, 4000 and 5000 s/mm 2 ) and 4 TEs (40, 80, 120 and 160 ms) for a total of 28 images. The dataset is acquired by applying diusion gradients along the z-axis. Chap.2 Background 36 Figure 2.18: Inversion-Recovery Multi-echo Spin-Echo (IR-MSE) imaging sequence. Figure 2.19: Diusion-weighted Spin-Echo imaging sequence. Chapter 3 Multidimensional Correlation Spectroscopic Imaging of Exponential Decays This chapter presents a novel multidimensional MRI approach that we developed to enable improved spatial mapping of microscopic tissue compartments from multidimensional correlation spectra estimation [34{38] 1 . An overview of the proposed method is provided in Section 3.1. Compartment modeling background is provided in Section 3.2. In Section 3.3, theoretical analysis is presented and then modeling and estimation strategies of the new method is illustrated. Finally, the summary of this chapter is provided in Section 3.4 3.1 Overview It is well known that many important biological changes to living tissues (due to development, aging, injury, disease, scientic intervention, etc.) initially occur at microscopic spatial scales. However, due to the limited sensitivity of MRI, it is very challenging to generate high-resolution MR images in a reasonable amount of time [61], meaning that conventional MR data is acquired 1 Some of the text in this chapter have been previously published in [35] and are copyright of the ISMRM. Some of the text in this chapter have been previously published in [36] and are copyright of the SPIE. 37 Chap.3 Multidimensional Correlation Spectroscopic Imaging of Exponential Decays 38 using relatively large imaging voxels that are much too big to directly interrogate microscopic tissue features. Despite the apparent resolution barrier, the scientic community has still invested decades of eort to try and infer microscopic-scale tissue information from MRI data. Since the direct approach is generally impractical, existing MRI-based methods have focused on indirect approaches that leverage the fact that certain MRI contrast mechanisms are sensitive to microscopic structure. For example, diusion MRI [2{4] leverages the fact that the MR signal can be sensitized to the random microscopic diusion of water molecules through tissue, while MRI relaxometry [5{7] leverages the fact that the relaxation characteristics of the MR signal are sensitive to the local physical and chemical microenvironment. While there has been considerable progress in using diusion and relaxation information to infer microstructure, existing methods still suer from ambiguities, and there remains a pressing need for new MRI methods that can estimate microstructural tissue compartments with higher sensitivity and specicity. Multicomponent modeling of relaxation or diusion decay curves [5, 6, 8{10, 12{20, 62{67] is one of the most common indirect approaches for resolving sub-voxel microstructure. The basic assumption of these approaches is that a large macroscopic voxel can be modeled as containing multiple dierent \compartments" corresponding to the water pools from distinct tissue microen- vironments, where each compartment is likely to exhibit distinct diusion or relaxation decay characteristics. As a result, neglecting any inter-compartmental exchange, the measurements from a single voxel can be modeled as a partial volume mixture of the distinct relaxation or diusion signatures (which generally take the form of exponential decays) that would be observed from each of the sub-voxel compartments. While some methods have used a small preselected number of compartments based on prior assumptions about tissue characteristics [13, 14, 65, 66], a more exible approach (which can accommodate an arbitrary and a priori unknown number of compartments) is to model the signal from a single voxel as a continuous distribution (or \spectrum") of dierent exponential decays [6, 8{10, 19, 62{65, 67]. In practical applications, these spectra generally exhibit distinct peaks that are usually ascribed to distinct compartments, and it is common to use spectral peak integrals to measure the contributions from each compartment. In imaging experiments, it is common to form a \spectroscopic image" that consists of a distinct decay spectrum for every spatial location, Chap.3 Multidimensional Correlation Spectroscopic Imaging of Exponential Decays 39 and the spatial maps of the spectral peak integrals can provide important additional insights into the spatial organization of the tissue compartments [8{10]. With 1D contrast encoding (i.e., designed to provide spectral information about a single decay parameter such as the T 1 , T 2 , or T 2 relaxation parameters or the apparent diusion coecient D), the inverse problem associated with estimating the exponential decay spectrum for a given voxel is sometimes described as a kind of 1D limited-data Inverse Laplace Transform (ILT)[21]. Unfortunately, it has been recognized for centuries that this inverse problem is fundamentally ill- posed and dicult to solve [21, 68]. Practically, this ill-posedness means that it is extremely dicult to separate tissue compartments that have similar exponential decay characteristics. To avoid some of the ill-posedeness associated with 1D spectrum estimation, multidimen- sional contrast encoding has been investigated in many applications [22{25, 27, 31, 32, 64, 69{74]. These approaches acquire a high-dimensional dataset that nonseparably encodes two or more decay parameters (e.g., jointD-T 2 contrast encoding [22, 23] or jointT 1 -T 2 contrast encoding [69, 70]), and then solve an inverse problem to estimate a multidimensional correlation spectrum that describes the joint distribution of these multiple decay parameters within each voxel. This work has success- fully demonstrated that multidimensional encoding and multidimensional spectrum estimation has benets over 1D approaches. The early work on multidimensional correlation spectroscopy of exponential decays has usu- ally reported correlation spectra representing large spatial volumes. These were obtained by either exciting and collecting signal from a large spatial volume at once with no additional spatial encod- ing, or by acquiring imaging data and then averaging the measured signal over regions of interest. Notably, there were no attempts to perform spatial mapping of the integrated spectral peaks un- til very recently. In principle, it would have been straightforward to acquire imaging data and generate spatial maps by performing voxel-by-voxel estimation of the multidimensional correlation spectra. However in practice, while the multidimensional inverse problem is not as ill-posed as the 1D problem, it is still somewhat ill-posed. As a consequence, conventional spectrum estimation approaches are associated with very onerous data quantity and quality requirements that would be hard to satisfy in a relatively short-duration imaging experiment with relatively good spatial resolution. To address the ill-posedeness of conventional spectrum estimation approaches, we present Chap.3 Multidimensional Correlation Spectroscopic Imaging of Exponential Decays 40 a novel high-dimensional MRI approach. The proposed method combines the multidimensional contrast encoding with spatial encoding and multidimensional spatially-regularized spectrum es- timation techniques. As can be seen in Fig. 3.1, the method will be performed by the following three steps: 1) multidimensional data acquisition with spatial encoding and joint contrast encod- ing, 2) Estimation of multidimensional correlation spectroscopic image, and 3) Spatial mapping of representative spectral peaks, in which each map presents the resolved sub-voxel structures. Figure 3.1: The framework of the proposed multidimensional correlation spectroscopic imaging approach. While a 4D spectroscopic imaging with diusion andT 2 relaxation is shown in the gure, other contrast mechanisms are possible even in higher dimension using the same framework. 3.2 Compartment Modeling Background Conventional multi-compartment modeling is based on the assumption that the observed signal is a linear mixture of signals originating from distinct compartments, where each compartment has a distinct set of intrinsic parameters that in uence image contrast. Multi-compartment modeling can be described using either a discrete compartment model or a continuum model, and our description will focus on a continuum model for the sake of generality. Under this model, it is assumed that the observed MR signal m( ) for a given volume of interest in the presence of experimental contrast encoding parameters can be expressed as m( ) = Z f()k( ;)d; (3.1) Chap.3 Multidimensional Correlation Spectroscopic Imaging of Exponential Decays 41 where f() is a density function that describes the amount of signal contribution associated with the parameters , andk( ;) is the ideal signal we expect to observe when using contrast encoding parameters in the presence of the parameters . Under the continuum model, f() is often called a spectrum, and the peaks of this spectrum are often interpreted as arising from distinct microscopic compartments. The ability to resolve and quantify distinct spectral peaks thus provides an indirect mechanism for understanding the microstructural features of a macroscopic volume. For example, the T 2 -spectrum and the T 2 -spectrum (described below) are commonly used to identify and quantify the myelin water compartment and the intra/extracellular water compartments in white matter of the brain. The ability to separate and quantify the myelin compartment from surrounding compartments has a variety of important applications, ranging from understanding the progression of normal brain white matter maturation and development across the lifespan, to improved understanding of demyelinating white matter diseases. The multi-compartment model from Eq. (3.1) has been applied in many dierent 1D settings. For example: Transverse relaxation time constant, T 2 . When =fT 2 g, the signal observed from a spin-echo sequence at echo time =fTEg is sometimes modeled as a multi-exponential decay[5, 8, 10]: m(TE) = Z f(T 2 )e TE=T 2 dT 2 : (3.2) A similar model has been used when using gradient-echo sequences instead of the spin-echo sequence. In this case, the only change to the signal model is that the transverse relaxation parameter T 2 is replaced with the slightly dierent parameter T 2 . Longitudinal relaxation time constant, T 1 . When =fT 1 g, the signal observed from a saturation recovery spin-echo sequence with repetition time =fTRg is sometimes modeled as a multi-exponential signal recovery [11]: m(TR) = Z f(T 1 )(1e TR=T 1 )dT 1 : (3.3) Chap.3 Multidimensional Correlation Spectroscopic Imaging of Exponential Decays 42 A T 1 -spectrum can be also estimated from the signal measured from a inversion recovery spin-echo sequence with inversion time =fTIg using the following model [9]: m(TI) = Z f(T 1 )(1 2e TI=T 1 )dT 1 : (3.4) Apparent diusion coecient, D. When =fDg, the signal observed from a Stejskal- Tanner pulsed gradient spin-echo sequence with diusion encoding parameter =fbg is sometimes also modeled as a multi-exponential decay[13, 75]: m(b) = Z f(D)e bD dD; (3.5) Estimating the 1D spectrum from either Eqs. (3.2),(3.3), (3.4), and (3.5) can be described as a form of 1D limited-data ILT, which is classically ill-posed as noted previously. Because of this ill-posedness, it is common to use additional constraints when solving these inverse problems to help stabilize the solution. It is especially common to assume that the spectra f(T 1 ) and f(T 2 ) should be nonnegative, which leads to a nonnegative least squares (NNLS) [76] formulation of the inverse problem [63]. Multidimensional extensions are easy to construct by combining multiple 1D models together. For example in the 2D case, existing multidimensional methods have used the same contrast model from Eq. (3.1) under various choices of (e.g., =fD;T 1 g,fD;T 2 g,fT 1 ;T 2 g, etc.) and corre- sponding choices of (e.g., =fb;TRg,fb;TEg,fTR;TEg, etc.).[22{32] In the multidimensional case, the spectrum f() is often called a correlation spectrum to distinguish it from a conven- tional 1D spectrum. For example, a method designed for the case with =fD;T 2 g has been called diusion-relaxation correlation spectroscopy (DR-COSY) [22, 23] with f(D;T 2 ) called the diusion-relaxation correlation spectrum. It is important to note that multidimensional correlation spectra almost always contain more information content than 1D spectra. This is easy to see based on the fact that most 1D spectra can be interpreted as some form of tomographic projection of a corresponding higher-dimensional correlation spectrum (e.g., we should have that f(D) = R f(D;T 2 )dT 2 in the ideal case). These multidimensional correlation spectra are also often much easier to estimate than the 1D spectra (assuming a corresponding shift in the data acquisition strategy), and the improved compartmental Chap.3 Multidimensional Correlation Spectroscopic Imaging of Exponential Decays 43 resolving power of multidimensional methods over 1D methods has been demonstrated consistently in many applications.[22{25, 27, 31, 32, 64, 69{74] 3.3 From Correlation Spectroscopy To Correlation Spectroscopic Imaging 3.3.1 Estimation Theoretic Analysis of Encoding in Higher Dimensions Multidimensional correlation spectroscopy methods are based on the synergistic higher-dimensional combination of these kinds of lower dimensional models. For example, the signal for 2D T 1 -T 2 correlation spectroscopy using an inversion-recovery spin-echo pulse sequence [72] can be modeled as m(TE;TI) = Z Z f(T 1 ;T 2 )(1 2e TI=T 1 )e TE=T 2 dT 1 dT 2 ; (3.6) where it now becomes necessary to acquire 2D contrast-encoded data m(TE;TI) in order to esti- mate the 2D T 1 -T 2 correlation spectrum f(T 1 ;T 2 ). This 2D approach is easily generalized to even higher dimensions, e.g., by using 3D contrast encoding with a 3D T 1 -T 2 -D correlation spectrum or by using 4D contrast encoding with a 4DT 1 -T 2 -D-T 2 correlation spectrum. In all of these settings, classical multidimenstional spectrum estimation approaches still frequently rely on nonnegativity assumptions and NNLS tting, just like for the 1D case. An important advantage of multidimensional contrast encoding relative to 1D contrast en- coding is reduced ill-posedness of the inverse problem. This leads to an improved capability to successfully resolve multiple compartments, even if some of the compartments possess similar de- cay parameters. For illustration, consider the toy scenario depicted in Fig. 3.2 in which a voxel consists of three distinct compartments, where compartment 1 has T 2 = 70 ms and T 1 = 750 ms, compartment 2 has T 2 = 100 ms and T 1 = 700 ms, and compartment 3 has T 2 = 110 ms and T 1 = 1000 ms. These three spectral peaks are well-separated in the 2D T 1 -T 2 space, and therefore may be relatively straightforward to resolve. On the other hand, they are not well-separated when only considering the 1D T 2 dimension (i.e., compartments 2 and 3 may be hard to resolve because they have similar T 2 values) or only considering the 1D T 1 dimension (i.e., compartments 1 and 2 may be hard to resolve because they have similar T 1 values). Chap.3 Multidimensional Correlation Spectroscopic Imaging of Exponential Decays 44 Figure 3.2: Toy illustration of the advantage of 2D multidimensional correlation spectroscopy over 1D spectroscopy. Ground truth values of three spectral peaks are shown in (a) a 2D T 1 -T 2 spectrum (left: a 3D plot and right: a 2D contour plot), (b) the corresponding 1D T 1 spectrum, and (c) the corresponding 1DT 2 spectrum. While the three peaks can all be successfully resolved in these ground truth spectra, real experiments will experience degraded spectral resolution because of nite sampling and noise. When resolution is degraded, the three peaks can still be easily resolved in (d) the 2D spectrum, though are no longer well-resolved in (e,f) either of the 1D spectra. While previous literature has conrmed the advantage of multidimensional correlation spec- troscopy in this setting empirically [24, 27, 31, 32, 35, 64, 69{72, 77], we nd it instructive to examine the dierence between 1D and multidimensional approaches from an estimation theoretic perspective. Our theoretical characterization is based on the Cram er-Rao Lower bound (CRLB), which is a theoretical lower bound on the variance of an unbiased estimator for an unknown param- eter of interest [78], and is frequently used to compare and optimize dierent experiment protocols in a variety of quantitative MR applications [79{85]. Since the CRLB depends on the specic values of the model parameters, we perform an illustrative analysis for the case of estimating the same toy three-compartment model described above. In this case, the ideal noiseless data for a given set of encoding parameters (TE;TI) is given by m(TE;TI) = 3 X s=1 f s 1 2e TI=T s 1 e TE=T s 2 ; (3.7) Chap.3 Multidimensional Correlation Spectroscopic Imaging of Exponential Decays 45 where the spin density parameters f s were all set equal to 1, and the T s 1 and T s 2 parameters were set to the T 1 and T 2 parameters given above for the three compartment model. Assuming white Gaussian noise, the CRLB is obtained by rst computing the Fisher information matrix for the unknown parameters, and then computing the inverse of the Fisher information matrix [78]. When computing the CRLB for 2D T 1 -T 2 correlation spectroscopy, we assumed that the number of compartments was known a priori, such that there were 9 unknown parameters to be estimated (i.e., thef s ,T s 1 , andT s 2 parameters for each compartment). For 1DT 1 orT 2 spectroscopy, we instead computed CRLBs assuming that there were 6 unknown parameters to be estimated (e.g., only thef s andT s 1 need to be estimated for each compartment in theT 1 -relaxometry case). The 2D T 1 -T 2 correlation spectroscopy acquisition was assumed to use an inversion-recovery preparation for T 1 contrast encoding and a Carr-Purcell-Meiboom-Gill (CPMG) multi-spin echo sequence for T 2 contrast encoding, with data sampled at every combination of 7 inversion times (TI = 0, 100, 200, 400, 700, 1000 and 2000) and 15 echo times (TEs ranging from 7.5 ms to 217.5 ms in 15 ms increments) for a total of 7 15 = 105 contrast encodings. We compared against a conventional inversion recovery spin-echo sequence for T 1 relaxometry, using the same 7 inversion times used for the 2D case. In this case, the expected scan time would be the same as for 2D case when the both sequences use the same TR. We also compared against conventional T 2 relaxometry using a standard CPMG-based multi-spin echo sequence with 32 echo times (TEs ranging from 10 ms to 320 ms). We assumed that this data was averaged 7 times, so that the experiment duration will match the duration of both the 2D correlation spectroscopy and T 1 relaxometry experiments. For all three experiments, the noise standard deviation was assumed to be the same. Comparing 2DT 1 -T 2 correlation spectroscopy against 1DT 2 relaxometry, the CRLB calcula- tion indicates that the lower bound on the standard deviation (i.e., the square root of the CRLB) achieved for T 2 with 2D correlation spectroscopy was 3.7710 2 times smaller for the rst com- partment, 1:08 10 3 times smaller for the second compartment, and 2:29 10 3 times smaller for the third compartment. Comparing 2D T 1 -T 2 correlation spectroscopy against 1D T 1 relaxometry, the CRLB calculation indicates that the lower bound on the standard deviation achieved for T 1 with 2D correlation spectroscopy was 9.1110 4 times smaller for the rst compartment, 2:21 10 4 times smaller for the second compartment, and 2:1010 3 times smaller for the third compartment. Chap.3 Multidimensional Correlation Spectroscopic Imaging of Exponential Decays 46 As can been seen, the CRLB values for multidimensional correlation spectroscopy are orders-of- magnitude lower than either of the conventional 1D methods, and imply that the 1D experiments would require from millions (for T 2 ) to billions (for T 1 ) times more data averaging to achieve the same CRLBs as the 2D experiment. This calculation helps to quantify the substantial estimation theoretic improvements oered by multidimensional encoding relative to 1D encoding. 3.3.2 Correlation Spectroscopic Imaging The previous section demonstrated considerable advantages for multidimensional encoding over 1D encoding, and this advantage can be sucient for experiments that are designed to generate spec- tra from large spatial volumes. However, conventional multidimensional correlation spectroscopic imaging experiments would still have relatively onerous data quantity and quality requirements if spectral reconstruction is performed voxel-by-voxel with high spatial resolution. In this section, we will describe new modeling and estimation schemes of our proposed approach. 3.3.2.1 Spatial-Spectral Modeling If we assume that a multidimensional dataset is acquired with high-dimensional contrast encoding and spatial encoding (i.e., a set of multi-contrast images), then the spatial-spectral signal model is give by: m(r; ) = Z f(r;)k( ;)d; (3.8) where m(r; ) represents the image acquired with a set of contrast encoding parameters as a function of the spatial location r, and f(r;) represents the high-dimensional spectroscopic image that is comprised of a full multidimensional correlation spectrum of at every spatial location. In the case where 2D imaging experiments are performed with diusion, T 1 andT 2 contrast encoding, the measured data m(x;y;TI;TE;b) and the spectroscopic image f(x;y;T 1 ;T 2 ;D) would be 5D while the spectroscopic image would be 6D for a 3D imaging experiment. It is worth noting that the ideal contrast modelk( ;) is not required to be separable and exponential, and it can be any function that is consistent with the physics of MRI data acquisition. One potential benet of spatial-spectral modeling is that, if the spectral peak locations (but not necessarily the spin density associated with each spectral peak) are assumed to be similar for neighboring voxels, then using the data from multiple voxels to estimate the location of the shared Chap.3 Multidimensional Correlation Spectroscopic Imaging of Exponential Decays 47 spectral peak can substantially reduce the ill-posedness of the estimation problem. In particular, analysis of a related problem has shown that joint spatial-spectral modeling (i.e., with the spectra for all voxels estimated simultaneously) has the potential to reduce the the CRLB by an order of magnitude or more relative to uncoupled voxel-by-voxel estimation [86]. The benets of spatial- spectral modeling have also been observed empirically in PET parameter estimation [86] and 1D relaxometry [87, 88]. To illustrate this improvement from an estimation theoretic perspective, it is instructive to again revisit the three compartment toy example from Section 3.3.1. However, instead of considering a single voxel in isolation, we now consider a scenario involving the simultaneous estimation of 2D correlation spectra from three dierent voxels. For this particular toy example, we will assume that each of these three voxels has the same three compartments as in Section 3.3.1 with the sameT 1 and T 2 parameter values. However, we will assume that the volume fractions for each compartment are distinct (i.e., f 1 =f 2 =f 3 = 1 for voxel 1, f 1 = 0:8, f 2 = 0:6 and f 3 = 1:8 for voxel 2, and f 1 = 2, f 2 = 0:5, andf 3 = 0:5 for voxel 3). Using the CRLB and the same 2D experimental paradigm from Section 3.3.1, we compare the voxel-by-voxel estimation strategy (where the relaxation parameters are not assumed to be the same for dierent voxels) against a strategy that is aware that the voxels share the same relaxation parameters. Our CRLB analysis shows that the spatial-spectral approach has CRLBs that range from 48 to 220 lower (depending on the specic parameter) relative to the voxel-by-voxel case. The voxel-by-voxel approach would thus require hundreds of additional averages to match the good estimation theoretic characteristics of the spatial-spectral approach, again indicating a substantial reduction in the ill-posedness of the spectral estimation problem. While this toy example is exaggerated (since it depends on the unrealistic and very strong prior information that the compartments in all voxels have identical relaxation characteristics), it is still indicative of the potential benets of spatial-spectral modeling. Importantly, this additional reduction in ill-posedness implies that data quality and quantity requirements can be relaxed, which can enable high quality correlation spectroscopic imaging experiments from much shorter experiments with fewer averages and fewer contrast encodings. Chap.3 Multidimensional Correlation Spectroscopic Imaging of Exponential Decays 48 3.3.2.2 High-dimensional Correlation Spectroscopic Image Estimation Based on the theory presented in the previous section, there are clear potential advantages to estimating a high-dimensional spectroscopic image from high-dimensional data using an estimation procedure that incorporates some form of spatial constraints. However, we generally don't want to make the same very strong assumptions that were used in previous toy examples. In the following, we describe the regularization based approach that we have developed. While the spatial-spectral model in Eq. (3.8) is continuous, we will use a dictionary-based discrete model for practical numerical implementation, as in conventional 1D and 2D diusion and relaxation spectroscopy methods [6, 8{10, 24, 27, 31, 32, 62{64, 69{72] are replaced by standard Riemann sum approximations of the form: m(r n ; p ) = Q X q=1 w q f(r n ; q )k( p ; q ); (3.9) forn = 1;:::;N andp = 1;:::;P . In this expression, it is assumed that that there are N dierent voxel locations r n and that full images are acquired withP dierent values of the contrast encoding parameters p . The continuum distribution is discretized over a large number Q of preselected q values, andw q is the normalization term (i.e., the quadrature weight) associated with the Riemann sum. The Riemann sum approximation can be made arbitrarily accurate by increasing the value of Q. Equation (3.9) is equivalently represented in matrix-vector form as m n = Kf n (3.10) for n = 1;:::;N. We have introduced m n 2 R P to represent the vectorization of all the data from the nth voxel, i.e., m(r n ; p ) for p = 1;:::;P ; the matrix K2 R PQ has entries [K] pq = w q k( p ; q ); and f n 2R Q is the vectorization of the multidimensional correlation spectrum from the nth voxel, i.e., f(r n ; q ) for q = 1;:::;Q. Given this discrete model, we perform estimation using the same nonnegativity constraints from classical relaxation spectroscopy methods, while also using a spatial smoothness constraint on Chap.3 Multidimensional Correlation Spectroscopic Imaging of Exponential Decays 49 the reconstructed 2D spectra [34{38]: n ^ f 1 ; ^ f 2 ;:::; ^ f N o = arg min ff i 2R Q g N i=1 " N X i=1 t i km i Kf i k 2 2 + X l2i kf i f l k 2 2 # (3.11) subject to [f i ] q 0 for8q = 1;:::;Q and8i = 1;:::;N. The rst term in this expression is a standard data consistency term for each voxel. The second term is a spatial regularization term that encourages the 2D correlation spectrum from one voxel to be similar to the 2D correlation spectra from neighboring voxels, where i is the index set for the voxels that are directly adjacent to the ith voxel. The parameter is a user-selected regularization parameter that controls the strength of the spatial regularization. In the data consistency term, the variables t i correspond to a spatial mask for the object, and are equal to 0 if the ith voxel is outside the object and are otherwise equal to one. This spatial mask is used to avoid tting spectra to noise-only voxels of the image, and prevents spectra within the object from being contaminated by noise when spatial smoothness constraints are imposed. Smoothness-based spatial regularization, as used in Eq. (3.11), is a classical constraint that is used in a wide range of imaging inverse problems [89]. This constraint is based on the principle that the spectra and spatial maps are likely to be spatially smooth, and can be viewed as a \soft" way of imposing the spatial-spectral constraints described in Sec. 3.3.2.1. In particular, the constraint encourages spectral similarity between adjacent voxels without forcing exact correspondence, which accommodates situations in which the spectral peak locations or lineshape characteristics vary gradually from voxel to voxel. In addition, this approach is not expected to fail in problematic ways if the spectrum from one voxel is very dierent from its neighbors (e.g., as will happen frequently along compartmental boundaries in the examples we show in the next section). In such cases, the use of spatial regularization is expected to behave gracefully, e.g., by blurring a feature that may originally have been sharper [89, 90]. Importantly, the regularization parameter can be varied to achieve a good balance between the ill-posedness of the estimation problem and the loss of spatial resolution, and there exist theoretical tools that can be used to quantify the trade-o between the two [91, 92]. If we let F2R QN and M2R PN represent the matrices whose ith columns respectively correspond to f i and m i , the optimization problem from Eq. (3.11) can also be more compactly be Chap.3 Multidimensional Correlation Spectroscopic Imaging of Exponential Decays 50 represented as ^ F = arg min F2R QN kMT KFTk 2 F +kFC H k 2 F ; (3.12) subject to [F] pq 0 for8p = 1;:::;P and8q = 1;:::;Q. In this expression,kk F denotes the standard Frobenius norm, T2R NN is a diagonal matrix whose ith diagonal entry contains the value of t i , C is the matrix operator that computes spatial nite dierences, and () H denotes the Hermitian operator (conjugate transpose). This representation is convenient for numerical optimization. The optimization problems in Eqs. (3.11) and (3.12) are convex, and there are many convex optimization methods that will nd the global solution from arbitrary initializations. Our work has generally used an alternating directions method of multipliers (ADMM) algorithm [93] to solve the nonnegativity-constrained optimization problem. We describe the steps of this algorithm in the Appendices A and B. 3.4 Summary In this chapter, we described a novel multidimensional correlation spectroscopic imaging method of exponential decays. The method combines multidimensional contrast encoding with spatial en- coding and multidimensional regularized estimation technique to enable substantially improved abilities of resolving microscopic sub-voxel compartments. This chapter theoretically demonstrated benets of the novel multidimensional method compared to existing methods, and described model- ing and estimation strategies of the proposed multidimensional method in a general framework. In the following chapters, the novel multidimensional method will be demonstrated with two specic congurations. Related Publications [Journal] Daeun Kim, Eamon K. Doyle, Jessica L. Wisnowski, Joong H. Kim, and Justin P. Haldar, \Diusion-Relaxation Correlation Spectroscopic Imaging: A Multidimen- sional Approach for Probing Microstructure," Magnetic Resonance in Medicine 78(6); p2236-2249, 2017. doi: http://dx.doi.org/10.1002/mrm.26629 Chap.3 Multidimensional Correlation Spectroscopic Imaging of Exponential Decays 51 [Preprint] Daeun Kim, Jessica L. Wisnowski, Christopher T. Nguyen, Jusitn P. Haldar, \Relaxation-Relaxation Correlation Spectroscopic Imaging (RR-CSI): Leverag- ing the Blessings of Dimensionality to Map In Vivo Microstructure," 2019. doi: https://arxiv.org/abs/1806.05752 [Conference procedding] Daeun Kim, Jessica L. Wisnowski, Justin P. Haldar, \MR Corre- lation Spectroscopic Imaging of Multidimensional Exponential Decays: Probing Microstructure with Diusion and Relaxation," SPIE Wavelets and Sparsity XVI, San Diego, 2017, p.10394D. doi:http://dx.doi.org/10.1117/12.2272833 Chapter 4 Diusion-Relaxation Correlation Spectroscopic Imaging (DR-CSI) There are multiple ways to implement a multidimensional correlation spectroscopic imaging method based on a variety of choices for multidimensional contrast encoding. As one specic example, this chapter presents diusion-relaxation correlation spectroscopic imaging (DR-CSI) where data is acquired with nonseparable encoding of diusion and T 2 relaxation 1 . We present spatial-spectral modeling for diusion-relaxation in Section 4.1. We evaluate the DR-CSI method with a wide range of experiments including numerical simulations, phantom experiments and ex-vivo animal experiments in Section 4.2. The discussion of this chapter is provided in Section 4.3 and the summary of this chapter is provided in Section 4.4 4.1 Diusion-Relaxation Correlation Modeling There are several dierent diusion decay models to choose from (e.g, isotropic versus anisotropic diusion, Gaussian versus non-Gaussian diusion within each compartment, etc.) that will have more or less sensitivity to various parameters of the acquisition pulse sequence (i.e., the diusion encoding b-value, the diusion time , the diusion encoding orientation, etc.). For the sake of simplicity and concreteness and without loss of generality, we will describe DR-CSI assuming that we are interested in estimating 2D diusion-relaxation correlation spectra, where the diusion decay 1 Some of the text and gures in this chapter have been previously published in [35] and are copyright of the ISMRM. 52 Chap.4 DR-CSI 53 is described using a simple 1D exponential decay model characterized by the apparent diusion coecient D, and the relaxation decay is similarly described by a 1D exponential decay model characterized by the transverse relaxation constant T 2 . Assuming 2D imaging (without loss of generality) with spatial coordinates (x;y), then the ideal signal model used by DR-CSI is given by m(x;y;b;TE) = ZZ f(x;y;D;T 2 )e bD e TE=T 2 ; (4.1) wherem(x;y;b;TE) is the 2D image acquired with diusion contrast encoding parameter b andT 2 contrast encoding parameter TE (echo time) andf(x;y;D;T 2 ) is the 4D spectroscopic image that comprises a full 2D diusion-relaxation correlation spectrum at every voxel. The 4D spectroscopic image is then estimated using the multidimensional regularized estimation method described in Section 3.3.2.2. For the estimation, the DR-CSI model is reformed by a dictionary-based discrete model m(x i ;y i ;b p ;TE p ) = Q X q=1 w q f(x i ;y i ;D q ;T q 2 )e bpDq e TEp=T q 2 ; (4.2) for i = 1;:::;N and p = 1;:::;P , which can be equivalently written in matrix-vector form as m i = Kf i (4.3) for i = 1;:::;N. In these expressions, we have assumed that there are N dierent voxel locations (x i ;y i ),i = 1;:::;N; the data is collected withP dierent combinations of diusion and relaxation encoding (b p ;TE p ), p = 1;:::;P ; we have sampled the continuum distribution at a large number Q of preselected (D q ;T q 2 ) values,q = 1;:::;Q; andw q is the density normalization term associated with the Riemann sum. The Riemann sum approximation of the continuum integral is well known to have arbitrarily good accuracy as we allow Q to be larger. In our matrix-vector expression, the vector m i 2R P is the set of all observed data samples from theith voxel, the matrix K2R PQ has entries [K] pq =w q e bpDq e TEp=T q 2 , and the vector f i 2R Q is the vectorized 2D diusion-relaxation correlation spectrum f(x i ;y i ;D q ;T q 2 ) from the ith voxel. From the matrix-vector formulation, the 4D spectroscopic image matrix F2R QN whose columns correspond to f i is estimated by solving the optimization problem in Eq. (3.12). Chap.4 DR-CSI 54 4.2 Experiments 4.2.1 Numerical Simulations To demonstrate the DR-CSI method, we constructed three dierent numerical simulation datasets, with each simulation corresponding to estimation of a 3-compartment model. The simulation setup is visualized in the left three columns of Fig. 4.1. In each case, the compartments were each given a dierent ideal single-peak D-T 2 spectrum f c (D;T 2 ) with a logarithmic 2D Gaussian lineshape, and a distinct binary spatial mask for each compartment a c (x;y). The ideal three compartment model was generated using simple summation: f(x;y;D;T 2 ) = 3 X c=1 a c (x;y)f c (D;T 2 ): (4.4) Noiseless DR-CSI measurements were simulated from this ideal model according to Eq. (4.1), and data was sampled at every combination of 7 b-values (0, 200, 500, 1000, 1500, 2500 and 5000 s/mm 2 ) and 7 TEs (99, 120, 150, 190, 240, 300 and 400 ms), for a total ofP = 49 images. Gaussian noise was subsequently added to these images. Finally, the image phase was discarded, resulting in magnitude images following a Rician distribution. The highest-SNR image (i.e., the least amount of diusion and relaxation decay, with b = 0 s/mm 2 and TE = 99 ms) had an SNR of 20, with SNR dened as the ratio between the average signal intensity and the noise standard deviation. For each multi-compartment simulation, 2D spectra were estimated at every voxel using the following parameters: the dictionary K was constructed using every combination of 80 D values (logarithmically distributed from 0.001 to 5m 2 /ms) and 80T 2 values (logarithmically distributed from 20 to 1000 ms) for a total of Q = 80 80 = 6400 dictionary elements. DR-CSI estimation was performed using = 0.1, = 0.1, and zero initialization. Figure 4.1 shows the results of the numerical simulations. DR-CSI estimated a 2D spectrum for every voxel in the image, and the average spectra (integrated across all voxels) are plotted in the fourth column of Fig. 4.1 for each simulation. In each case, the average spectra have three distinct peaks that are are largely consistent with the peaks from the ground truth spectra shown in the rst column of Fig. 4.1. Compared to the ground peaks, the estimated peaks have broader lineshapes. This lineshape broadening is expected to be a consequence of the nite spectral resolution associated Chap.4 DR-CSI 55 with nite data sampling and noise. The last 3 columns of Fig. 4.1 illustrate the spatial variations of each spectral peak appearing in the spatially-averaged spectrum. Specically, each image shows a spatial map of an integrated spectral peak. Despite the relatively high noise level in some of the simulated images, the three dierent original compartments are successfully separated in all three simulations. Ground truths a 20 10 2 10 3 T2 (msec) 10 -2 10 -1 10 0 D (μm 2 /msec) comp.1 comp.2 comp.3 Es!mated 2D spectra 20 10 2 10 3 T2 (msec) 10 -2 10 -1 10 0 D (μm 2 /msec) Spa!al maps comp.1 comp.2 comp.3 2.5 0 1.5 0 2 0 b c 20 10 2 10 3 T2 (msec) 10 -2 10 -1 10 0 D (μm 2 /msec) comp.1 comp.2 comp.3 4 0 4 0 20 10 2 10 3 T2 (msec) 10 -2 10 -1 10 0 D (μm 2 /msec) comp.1 comp.2 comp.3 20 10 2 10 3 T2 (msec) 10 -2 10 -1 10 0 D (μm 2 /msec) 20 10 2 10 3 T2 (msec) 10 -2 10 -1 10 0 D (μm 2 /msec) 10 0 1.5 0 5 0 1.5 0 2 0 5 0 4 0 Mul!-compartment images b = 0 , TE = 99 b = 5000, TE = 400 0.2 0 0.2 0 0.2 0 Figure4.1: DR-CSI simulation. (a-c) Numerical simulation setup and estimation results for three dierent multi-compartment datasets with dierent diusion-relaxation correlation characteristics. The rst column shows ground truth 2D diusion-relaxation correlation spectra (averaged across all voxels). The second and third columns show representative simulated diusion- and relaxation- encoded images (corresponding to b = 0 s/mm 2 and TE = 99 ms, and b = 5000 s/mm 2 and TE = 400 ms, respectively) that illustrate the geometry of compartmental overlap and the noise level. The fourth column shows 2D spectra estimated using DR-CSI (averaged across all voxels of the spectroscopic image), while the last three columns show spatial maps of the integrated spectral peaks from the estimated spectroscopic images. Chap.4 DR-CSI 56 4.2.2 Phantom Experiments To validate the DR-CSI method, we conducted phantom experiments. For comparison against the numerical simulation results, we custom-built a phantom possessing similar geometric structures as shown in Fig. 4.2. The letter-shaped structures were 3D printed, and following recommendations from the Wald group (https://phantoms.martinos.org/), we used liquid rubber for waterproong. Each compartment of the phantom was lled with a dierent mixture of Polyethylene Glycol (PEG) (Up and up powderlax; Target Corporation, Minneapolis, MN, USA; 3350 g/mol) and gadobutrol (Gadovist; Bayer Healthcare, Leverkusen, Germany; 1 mmol/ml) to respectively adjust the D and T 2 values [94]. The specications of each solution are given in Table 4.1. Figure 4.2: Phantom Table 4.1: DesiredD andT 2 values and corresponding concentrations of PEG and gadobutrol for each compartment of the custom-built phantom. S B C G M I U/bottle R Desired D (m 2 /ms) 0.1 0.3 0.8 2.1 PEG (mM) 210 145 73 0.18 Desired T 2 (ms) 80 160 350 80 350 80 160 350 gadobutrol (mM) 0.18 0.003 0 1.14 0.15 2.4 0.99 0.46 Chap.4 DR-CSI 57 DR-CSI datasets were acquired using a 3 T human MRI system (Achieva; Philips Healthcare, Best, The Netherlands) using a diusion-weighted spin-echo imaging sequence and an 8-channel head coil with the following parameters: TR = 11000 ms, voxel size = 3 mm 3 mm, matrix size = 64 40, slice thickness = 5 mm, and 33 slices. Contrast encoding used the same set of (b;TE) parameters that were used in the numerical simulations (P = 49). The data was acquired with thin slices so that there was no overlap between compartments within a single slice. Three dierent 3-compartment datasets were generated by overlaying and summing single-compartment data from dierent slices. Representative 3-compartment images are shown in the rst and second columns of Fig. 4.3(b-d). DR-CSI estimation was performed using the same parameters as the numerical simulations. We also compared DR-CSI against voxel-by-voxel DR-CSI (no spatial regularization), con- ventional 1D relaxation, and conventional 1D diusion methods. Voxel-by-voxel DR-CSI is similar to DR-COSY, in that it estimates a 2D spectrum for each voxel without the benets of spatial information from neighboring voxels. Voxel-by-voxel DR-CSI was implemented using the same problem formulation as Eq. (3.11) but without spatial regularization ( = 0). For the conventional 1D methods, relaxation spectra were estimated from the 7 dierent TEs acquired with no diusion weighting, and diusion spectra were estimated from the 7 dierent b-values acquired at the short- est echo time. The 1D relaxation and 1D diusion spectra were both estimated using Eq. (3.11), with the dictionary matrix K modied to include only relaxation or diusion decays respectively. For enhanced performance [87, 95], the 1D spectra were both estimated using spatial regularization. Fig. 4.3(b-d) shows DR-CSI results from the phantom experiment, with the spatially-averaged 2D spectra shown in the third column and the spatial maps of integrated spectral peaks shown in the last three columns. Accurate separation of the three superposed compartments is also achieved in all three of these cases, and the estimated spectral peak characteristics are largely consistent with the characteristics of the PEG-gadobutrol solutions listed in Table 4.1. Chap.4 DR-CSI 58 4 0 4 0 c b a 0 4 0 0 0 2 20 10 2 10 3 T2 (msec) 10 -2 10 -1 10 0 D (μm 2 /msec) comp.1 comp.2 comp.3 7 0 1.5 0 5 0 20 10 2 10 3 T2 (msec) 10 -2 10 -1 10 0 D (μm 2 /msec) comp.1 comp.2 comp.3 1.5 0 2 0 5 0 20 10 2 10 3 T2 (msec) 10 -2 10 -1 10 0 D (μm 2 /msec) comp.3 comp.2 comp.1 Es mated 2D spectra Spa al maps comp.1 comp.2 comp.3 2.5 0 1.5 0 2 0 4 Mul -compartment images b = 0 , TE = 99 b = 5000, TE = 400 0.1 0 0.1 0 0.1 0 Figure4.3: DR-CSI phantom results. (a-c) Example data and DR-CSI estimation results obtained with this phantom. The rst and second columns show representative diusion- and relaxation- encoded images (corresponding to b = 0 s/mm 2 and TE = 99 ms, and b = 5000 s/mm 2 and TE = 400 ms, respectively). The third column shows 2D spectra estimated using DR-CSI (averaged across all voxels of the spectroscopic image), while the last three columns show spatial maps of the integrated spectral peaks from the estimated spectroscopic images. Chap.4 DR-CSI 59 Figure 4.4: Spatially-varying DR-CSI spectra from the voxels corresponding to the red box drawn on the multi-compartment image in Fig. 4.3(c). Dierent peaks are color-coded by their corresponding compartments: green (S), red (R) and blue (I). Figure 4.4 shows individual DR-CSI spectra from the dataset of Fig. 4.3(c) at voxels where dierent compartmental geometries are superposed. This result shows an ability of DR-CSI to identify multiple components within each voxel. Notably, the spectral characteristics of real data are not as regular as they were in the previous numerical simulation results shown in Fig. 4.1, and demonstrate irregular lineshapes. This is likely due to a number of practical factors, including truly non-Gaussian spectral lineshapes, spatial variations in the concentrations of the PEG-gadobutrol mixtures, chemical shift artifacts from the PEG spectrum [94], and external factors like B0 inhomogeneity. Due to the way that the phantom was constructed, eld inhomogeneity was observed to be qualitatively much more substantial in this phantom compared to what is observed in typical biological tissues, which may contribute additional spatially-varying diusion weighting that is not modeled by our estimation scheme. Nevertheless, Chap.4 DR-CSI 60 we are still able to robustly separate distinct spectral peaks which are successfully mapped back to reveal the original compartmental geometry of the phantom. Figure 4.5 shows the results obtained with conventional 1D diusion, conventional 1D re- laxation, and voxel-by-voxel DR-COSY corresponding to the same multi-compartment data used in Fig. 4.3(c). The spatially-integrated spectra for the conventional 1D methods only show two distinct peaks, which is expected based on the diusion-relaxation characteristics of the compart- ments. Specically, as seen in Table 4.1 and Fig. 4.3(c), the `S' and `I' compartments have very similar relaxation parameters, while the `R' and `I' compartments have very similar diusion coef- cients. The spatial maps of the integrated peaks demonstrate that, as expected, it is dicult to separate compartments when the decay parameters are too similar to one another. The voxel-by- voxel DR-CSI results are more successful than the 1D approaches, but the spatially-averaged 2D spectrum contains more peaks than the original number of compartments, and the spatial maps of the integrated spectral peaks are not as accurate at reconstructing the true geometries of the orig- inal compartments. These results strongly demonstrate that both the multidimensional contrast encoding and multidimensional spectrum estimation components of DR-CSI are important, and that DR-CSI can oer substantially enhanced ability to resolve multiple overlapping compartments compared to the conventional approaches. Chap.4 DR-CSI 61 10 0 2 0 1.5 0 2 0 Integrated 2D spectra Conven!onal 1D Diffusion Integrated 1D spectrum Spa!al maps comp.1 comp.2 0 a 10 -1 D (μm 2 /msec) 10 0 0.02 0.06 0.04 0.08 10 -2 0 comp.1 comp.2 Conven!onal 1D Relaxa!on Integrated 1D spectrum Spa!al maps comp.1 comp.2 0 b DR-CSI without Spa!al Regulariza!on Spa!al maps comp.1 comp.2 c 0.04 10 2 T2 (msec) 10 3 0.01 0.03 0.02 comp.1 comp.2 0 20 20 10 2 10 3 T2 (msec) 10 -2 10 -1 10 0 D (μm 2 /msec) comp.1 comp.2 comp.4 comp.3 comp.5 6 0 5 0 3 0 1.5 0 10 0 comp.3 comp.4 comp.5 Figure4.5: Estimation results from (a) conventional 1D diusion, (b) conventional 1D relaxation, and (c) voxel-by-voxel DR-CSI (no spatial regularization). Each subgure shows (left) the estimated spectra (averaged across all voxels of the spectroscopic image) and (right) spatial maps of the integrated spectral peaks from the estimated spectroscopic image. Chap.4 DR-CSI 62 4.2.3 Ex Vivo Mouse Spinal Cord Experiments As an biological application example, six ex vivo mouse spinal cords (three sham controls and three with traumatic spinal cord injury as described in [96]) were scanned with a DR-CSI protocol. Diusion in the mouse spinal cord is highly oriented [97], which enables the use of 1D diusion encoding despite the anisotropic nature of the diusion process. We obtained two DR-CSI datasets for each cord, one with diusion encoding parallel and one with diusion encoding perpendicular to the axonal tracts. Datasets were acquired using a 4.7 T animal system (Agilent Inc., Palo Alto, CA) using a diusion-weighted spin-echo imaging sequence with the following parameters: TR = 2000 ms, voxel size = 0.78 m 0.78 m, matrix size = 96 96, slice thickness = 1 mm, and 5 slices. Data was sampled at every combination of 7 b-values (0, 500, 1000, 2000, 3000, 4000 and 5000 s/mm 2 ) and 4 TEs (40, 80, 120 and 160 ms) for a total of P = 28 images. Fig. 4.6 shows a representative full set of 28 DR-CSI contrast encodings from a single slice of a control cord. b = 0 unit: s/mm 2 TE = 40 TE = 80 TE = 120 TE = 160 unit: ms b = 500 b = 1000 b = 2000 b = 3000 b = 4000 b = 5000 Figure 4.6: The full set of P = 28 diusion- and relaxation-encoded images from a single slice of a representative control mouse spinal cord. The diusion encoding orientation was parallel to the axonal tracts. Chap.4 DR-CSI 63 For each dataset, 2D spectra were estimated using the following parameters: the dictionary K was constructed using every combination of 70 D values (logarithmically distributed from 0.01 to 5 m 2 /ms), and 70 T 2 values (logarithmically distributed form 3 to 300 ms) for a total of Q = 70 70 = 4900 dictionary elements. Similar to the previous section, DR-CSI based results were compared against voxel-by-voxel DR-CSI, conventional 1D relaxation, and conventional 1D diusion. Figure 4.7 shows DR-CSI results generated from all six of the ex vivo mouse spinal cord datasets acquired with a parallel diusion encoding orientation. The 2D spectra from the control cords consistently have two distinct well-resolved peaks, as well as a third weaker peak in between. The third peak may not be very visible in the spatially-averaged spectra, although its existence is more clear in many of the 2D spectra obtained from individual voxels. In contrast, the 2D spectra for the injured cords contain an additional peak that was not present in the control spectra. 10 1 10 2 T2 (msec) 10 -1 10 0 D (μm 2 /msec) 10 1 10 2 T2 (msec) 10 -1 10 0 D (μm 2 /msec) a Control sub. 1 Control sub. 2 Control sub. 3 b Injured sub. 1 Injured sub. 2 Injured sub. 3 10 1 10 2 T2 (msec) 10 -1 10 0 D (μm 2 /msec) 10 1 10 2 T2 (msec) 10 -1 10 0 D (μm 2 /msec) 10 1 10 2 T2 (msec) 10 -1 10 0 D (μm 2 /msec) 10 1 10 2 T2 (msec) 10 -1 10 0 D (μm 2 /msec) Figure 4.7: Reference images and spatially-averaged DR-CSI spectra (corresponding to the par- allel diusion encoding orientation) from all of the (a) control spinal cords and (b) injured mouse spinal cords. To highlight the spatially-varying nature of the estimated DR-CSI spectroscopic images, Fig. 4.8 shows spatially-varying spectra from a region from the white matter (WM) - gray matter (GM) boundary of the rst control cord, indicated with a red box in Fig. 4.7(a). As can be seen, the Chap.4 DR-CSI 64 spatial distribution of the spectra clearly depicts the the transition between WM and GM, with one distinct peak in the WM region, a dierent distinct peak in the GM region, and a combination of the two peaks in the partial volume region at the boundary. Notably, it is not possible to gain this kind of insight from 2D spectroscopy methods like DR-COSY, which do not use spatial encoding. 10 1 10 2 T2 (ms) 10 -1 D (μm 2 /ms) 10 0 10 1 10 2 T2 (ms) 10 -1 D (μm 2 /ms) 10 0 10 1 10 2 T2 (ms) 10 -1 D (μm 2 /ms) 10 0 10 1 10 2 T2 (ms) 10 -1 D (μm 2 /ms) 10 0 10 1 10 2 T2 (ms) 10 -1 D (μm 2 /ms) 10 0 10 1 10 2 T2 (ms) 10 -1 D (μm 2 /ms) 10 0 10 1 10 2 T2 (ms) 10 -1 D (μm 2 /ms) 10 0 10 1 10 2 T2 (ms) 10 -1 D (μm 2 /ms) 10 0 10 1 10 2 T2 (ms) 10 -1 D (μm 2 /ms) 10 0 10 1 10 2 T2 (ms) 10 -1 D (μm 2 /ms) 10 0 10 1 10 2 T2 (ms) 10 -1 D (μm 2 /ms) 10 0 10 1 10 2 T2 (ms) 10 -1 D (μm 2 /ms) 10 0 10 1 10 2 T2 (ms) 10 -1 D (μm 2 /ms) 10 0 10 1 10 2 T2 (ms) 10 -1 D (μm 2 /ms) 10 0 10 1 10 2 T2 (ms) 10 -1 D (μm 2 /ms) 10 0 10 1 10 2 T2 (ms) 10 -1 D (μm 2 /ms) 10 0 GM WM Figure 4.8: Spatially-varying DR-CSI spectra from the boundary between WM and GM in the control cord (subject 1), corresponding to the red box drawn on the anatomical image in Fig. 4.7. Chap.4 DR-CSI 65 Representative spatially-averaged DR-CSI spectra from the perpendicular diusion encoding orientation are shown in Fig. 4.9. The spectral peaks in this case are more closely spaced than they were for the parallel orientation, and, for example, it is more dicult to visually separate multiple peaks from the spatially-averaged 2D spectrum from the control cord. Nevertheless, there is a clear additional peak in the spectrum from the injured cord that was not present in the control cord spectrum. 10 1 10 2 T2 (msec) 10 -1 10 0 D (μm 2 /msec) 10 1 10 2 T2 (msec) 10 -1 10 0 D (μm 2 /msec) a Control sub. 1 Control sub. 2 Control sub. 3 b Injured sub. 1 Injured sub. 2 Injured sub. 3 10 1 10 2 T2 (msec) 10 -1 10 0 D (μm 2 /msec) 10 1 10 2 T2 (msec) 10 -1 10 0 D (μm 2 /msec) 10 1 10 2 T2 (msec) 10 -1 10 0 D (μm 2 /msec) 10 1 10 2 T2 (msec) 10 -1 10 0 D (μm 2 /msec) Figure 4.9: Representative spatially-averaged DR-CSI spectra (corresponding to the perpendicu- lar diusion encoding orientation) from (a) control and (b) injured mouse spinal cords. Representative spatial maps generated by integrating spectral regions from the DR-CSI spec- troscopic images with the parallel diusion encoding orientation are shown in Fig. 4.10 (The results from the perpendicular diusion encoding orientation are shown in Fig. 4.11). There is no ground truth for this case, and we cannot denitively associate each of these estimated components with distinct components of the tissue microstructure without additional investigations that are beyond the scope of this paper (e.g., histology). However, as can be seen, the spatial maps for the con- trol cords seem consistent with the known anatomy of the spinal cord. Specically, in all cases, component 1 appears to correspond to white matter, component 2 appears to correspond to gray matter, and component 3 also appears to correspond to gray matter, but with a larger signal from Chap.4 DR-CSI 66 the dorsal gray matter than from the ventral gray matter (except in cases of injury). Notably, component 4 indicates a compartment that is substantial in the injured cords but is not present in the control cords, and is likely to re ect a microstructural change resulting from the injury. In addition, as can be seen from the composite images, the spatial maps for dierent components have considerable spatial overlap, which suggests that DR-CSI is successfully disentangling partial vol- ume contributions from multiple tissue compartments within the same voxel, as would be expected based on the previous numerical simulation and phantom experiment results. 0 1 0 0 1 0 0.6 0.3 10 1 10 2 T2 (msec) 10 -1 10 0 D (μm 2 /msec) 1 2 3 4 10 1 10 2 T2 (msec) 10 -1 10 0 D (μm 2 /msec) 1 2 3 4 comp.1 comp.2 comp.3 comp.4 composite control injured Spa al maps Spectral regions a b control injured Figure 4.10: Spatial maps of the integrated spectral peaks from the DR-CSI reconstruction from data acquired with parallel diusion orientation. (a) The spectral regions that are integrated to generate the spatial maps (red: comp.1, blue: comp.2, green: comp.3 and yellow: comp.4). (b) The spatial maps from the spectral regions. Chap.4 DR-CSI 67 comp.1 comp.2 comp.3 comp.4 composite control injured 0 1 0 0 1 0 1 0.3 Spa al maps Spectral regions 2 3 10 1 10 2 T2 (msec) 10 -1 10 0 D (μm 2 /msec) 1 4 2 3 10 1 10 2 T2 (msec) 10 -1 10 0 D (μm 2 /msec) 1 4 a b control injured Figure 4.11: Spatial maps of the integrated spectral peaks from the DR-CSI reconstruction from data acquired with parallel diusion orientation. (a) The spectral regions that are integrated to generate the spatial maps (red: comp.1, blue: comp.2, green: comp.3 and yellow: comp.4). (b) The spatial maps corresponding to the spectral regions. Chap.4 DR-CSI 68 For comparison, Fig. 4.12 shows the compartment estimation results from conventional meth- ods. There is considerable ambiguity in the 1D diusion and 1D relaxation spectra, which resolve substantially fewer peaks than DR-CSI. Voxel-by-voxel DR-CSI (no spatial regularization) yields spectra with considerably less structure than the proposed approach, and yields spatial maps that are dicult to interpret in a meaningful way. Spaal maps comp.1 comp.2 a 1 0.6 10 -1 D (μm 2 /msec) 10 0 0.01 0.02 0.03 10 -2 0 comp.3 comp.2 comp.1 comp.3 comp.4 Integrated 1D spectrum b 10 1 T2 (msec) 10 2 0.02 0.04 0.06 0 comp.3 comp.2 comp.1 0 1 0 1 0 1 Integrated 2D spectrum 10 1 10 2 T2 (msec) 10 -1 10 0 D (μm 2 /msec) 1 2 4 3 composite 0 1 0 1 0 0 0 0 0 0 1 1 1 0 1 0 1 0 1 0 1 0 1 0 1 diffusion relaxaon voxel-by -voxel DR-CSI diffusion relaxaon voxel-by-voxel DR-CSI Figure 4.12: Estimation results from conventional 1D diusion, conventional 1D relaxation, and voxel-by-voxel DR-CSI (no spatial regularization). Each subgure shows (a) the estimated spectra (averaged across all voxels of the spectroscopic image) with the spectral regions to generate the spatial maps (red: comp.1, blue: comp.2, yellow: comp.3 and green: comp.4) and (b) spatial maps from the integrated spectral regions. Chap.4 DR-CSI 69 4.3 Discussion Our extensive simulation, phantom, and ex vivo mouse spinal cord results strongly conrm the hypothesis that DR-CSI can oer substantial advantages in resolving overlapping tissue compart- ments compared to traditional methods. However, a notable limitation of DR-CSI is that it requires substantially more data than conventional methods, which can lead to long data acquisition times. The phantom data shown in this paper took about seven hours to acquire, while each mouse spinal cord dataset took about one and half hours to acquire. However, it should be noted that our present pulse sequence implementations were quite basic, and did not use fast imaging methods like echo-planar imaging, parallel imaging, or compressed sensing [98], or fast contrast encoding methods like the CPMG pulse sequence [99]. For example, using a simple EPI readout would be expected to make the acquisition time 40-96 times faster for the specic experiments described in this work. In addition, we have not yet made any attempts at optimizing the distribution of (b;TE) contrast parameters that we sample data at, and we anticipate that there are substantial opportunities for improvement, especially given the success of compressed-sensing type methods in non-imaging multi-dimensional spectroscopy experiments [31, 32] and emerging fast contrast en- coding methods like multi-compartment MR ngerprinting [100, 101]. We anticipate that optimal experiment design for DR-CSI and related methods will be a promising area for future research. Another limitation of DR-CSI is that, while we can observe spectral peaks that correspond to dierent microstructural compartments, it can be dicult to associate a given spectral peak with a specic feature of the tissue microstructure without additional prior information. For example, in the results shown in Fig. 4.10 and Supporting Fig. 4.11, we had to rely on prior knowledge of mouse spinal cord anatomy in order to hypothesize about the association between dierent spectral peaks with dierent apparent tissue types. This limitation is common to almost every existing MR- based compartment estimation method, and a popular approach to resolving this correspondence problem is to also perform histological analysis of the tissue. We believe that combining DR-CSI with histological analysis is a promising direction for future research. While we presented DR-CSI using a simple 2D model involving a 1D diusion coecient and a 1D T 2 relaxation coecient with contrast manipulated through the b-value and the TE, it is important to emphasize that this model was only assumed for the sake of providing a simple proof-of-principle demonstration of the power of the technique. DR-CSI can potentially be used in Chap.4 DR-CSI 70 a substantially broader range of settings. Alternative types of multidimensional spectral estimation are also possible within this framework, including T 1 -T 2 , D-D, D-T 1, and D-T 1 -T 2 spectroscopic imaging [24, 31, 32, 64]. It is also straightforward to include more complicated data acquisition and biophysical compartment models that account for non-exponential signal variations, diusion anisotropy, diusion time dependence, water exchange, imperfect ip angles, B0 eld inhomogene- ity, etc. We believe that future explorations of these opportunities will expand the capabilities of DR-CSI even further. 4.4 Summary This chapter presented and evaluated DR-CSI, a novel diusion-relaxation correlation spectroscopic imaging method that combines ideas from 2D diusion-relaxation spectroscopy with imaging gra- dients and an advanced 4D joint spectral estimation scheme. We demonstrated that DR-CSI has powerful capabilities for resolving spatially-overlapping tissue compartments using numerical sim- ulations, phantom experiments, and data from an animal model. We expect that the DR-CSI technique, along with its future evolutions, may substantially expand the role of MRI in probing important features of tissue microstructure that have previously been inaccessible to traditional MR methods. Related Publications [Journal] Daeun Kim, Eamon K. Doyle, Jessica L. Wisnowski, Joong H. Kim, and Justin P. Haldar, \Diusion-Relaxation Correlation Spectroscopic Imaging: A Multidimen- sional Approach for Probing Microstructure," Magnetic Resonance in Medicine 78(6); p2236-2249, 2017. doi: http://dx.doi.org/10.1002/mrm.26629 [Conference abstract] Daeun Kim, Eamon K. Doyle, Jessica L. Wisnowski, Justin P. Haldar, \Phantom Validation of Diusion-Relaxation Correlation Spectroscopic Imaging (DR-CSI)," International Society for Magnetic Resonance in Medicine, 25th Annual Meet- ing, Honolulu, 2017, p609. Chap.4 DR-CSI 71 [Conference abstract] Daeun Kim, Joong H. Kim, Justin P. Haldar, \Diusion-Relaxation Correlation Spectroscopic Imaging: Diusion-Relaxation Correlation Spectro- scopic Imaging (DR-CSI): An Enhanced Approach to Imaging Microstructure," International Society for Magnetic Resonance in Medicine, 24th Annual Meeting, Singapore, 2016, p660. Chapter 5 T 1 Relaxation - T 2 Relaxation Correlation Spectroscopic Imaging (RR-CSI) While we have demonstrated the novel multidimensional approach with DR-CSI in the previ- ous chapter, we can also demonstrate the method with other contrast mechanisms. This chapter presents T 1 relaxation -T 2 relaxation correlation spectroscopic imaging (RR-CSI) 1 . We describe spatial-spectral modeling for T 1 relaxation-T 2 relaxation in Section 5.1. We evaluate the RR-CSI method with numerical simulations and real experiments including in-vivo human brain experiments that we believe the rst demonstration of multidimensional spectroscopic imaging approaches in Section 5.2. The discussion of this chapter is provided in Section 5.3 and the summary of this chapter is provided in Section 5.4. 5.1 T 1 Relaxation - T 2 Relaxation Correlation Modeling RR-CSI can be achieved by choosing =fT 1 ;T 2 g and selecting an MR acquisition based on an inversion recovery spin echo (IR-SE) sequence with contrast encoding parameters =fTI;TEg in Eq. (3.1). Assuming 2D imaging experiments, this leads to the following ideal RR-CSI signal 1 Some of the text and gures are previously published in [36, 37] and are copyright of the SPIE and the IEEE. 72 Chap.5 RR-CSI 73 model without loss of generality: m(x;y;TE;TI) = Z Z f(x;y;T 1 ;T 2 )(1 2e TI=T 1 )e TE=T 2 dT 1 dT 2 ; (5.1) where m(x;y;TE;TI) represents the 2D image acquired with contrast encoding parameters TI (inversion time) andTE (echo time), andf(x;y;T 1 ;T 2 ) represents the 4D spectroscopic image that is comprised of a full 2D T 1 -T 2 relaxation correlation spectrum at every spatial location. Similar to DR-CSI, we will use a dictionary-based discrete model given by m(x i ;y i ;TE p ;TI p ) = Q X q=1 w q f(x i ;y i ;T q 1 ;T q 2 )(1 2e TIp=T q 1 )e TEp=T q 2 ; (5.2) for8i = 1;:::;N and8p = 1;:::;P . In this expression, N is the number of voxels in the image; it is assumed that we acquire P dierent contrast encodings (TE p ;TI p ); we have modeled the relaxation distribution using a dictionary with Q elements, where the qth element corresponds to the relaxation parameters (T q 1 ;T q 2 ); and w q is the density normalization term (i.e., the numerical quadrature weights) required for accurate approximation of the continuous integral using a nite discrete sum. This discrete model can be equivalently represented in matrix form as m i = Kf i ; (5.3) for i = 1;:::;N, where the vector m i 2 R P contains all the measured contrast-encoded data corresponding to the ith voxel and has pth entry [m i ] p = m(r i ;TE p ;TI p ); the dictionary matrix K 2 R PQ has entries [K] pq = w q (1 2e TIp=T q 1 )e TEp=T q 2 ; and f i 2 R Q is the 2D spectrum corresponding to the ith voxel of the high-dimensional spectroscopic image and has qth entry [f i ] q =f(r i ;T q 1 ;T q 2 ). Given this discrete model, a 4D spectroscopic image is subsequently estimated by solving the optimization problem in Eq. (3.12) same as in DR-CSI. 5.2 Experiments As illustrative examples of multidimensional correlation spectroscopic imaging, we will demonstrate 2D T 1 -T 2 correlation spectroscopic imaging using numerical simulations, real experiments with a Chap.5 RR-CSI 74 pumpkin, and several real experiments with in vivo human brains. While advanced experimen- tal protocols would likely enable improved experimental eciency, we have focused on a simple proof-of-principle implementation in which the 2D T 1 -T 2 data is acquired with a basic inversion- recovery multi-echo spin-echo (IR-MSE) pulse sequence and 2D spatial encoding that are described in Fig. 2.18. 5.2.1 Numerical Simulations Numerical simulations are valuable for understanding the characteristics of multidimensional corre- lation spectroscopic image estimation, because (unlike for real experiments), we have a ground truth we can compare the estimation results against. In this simulation, a gold standard spectroscopic image was constructed using a three-compartment model according to f(x;y;T 1 ;T 2 ) = 3 X c=1 a c (x;y)f c (T 1 ;T 2 ); (5.4) where a c (x;y) is the spatial distribution and f c (T 1 ;T 2 ) is the spectrum for the cth compartment, as shown in Figure 5.1(a). The spectra were generated using a 2D Gaussian spectral lineshape with dierent spectral peak locations (compartment 1: (T 1 ;T 2 )=(70ms, 750ms); compartment 2: (T 1 ;T 2 )=(100ms, 700ms); compartment 3: (T 1 ;T 2 )=(110ms, 1000ms)). Note that the three compartments are dicult to distinguish using 1DT 1 relaxation orT 2 relaxation approaches because the three compartments have similar T 1 orT 2 values to one another. Noiseless data was generated using every combination of 7 inversion times (TI = 0, 100, 200, 400, 700, 1000 and 2000 ms) and 15 echo times (TEs ranging from 7.5 ms to 217.5 ms in 15 ms increments) for a total of P = 105 contrast encodings. Subsequently, Gaussian noise was added to the noiseless data, and magnitudes were taken leading to Rician-distributed data. In the dataset, image SNRs range from 3.83 to 200 (SNRs were computed separately for each contrast-encoded image as the ratio between the average per-pixel signal intensity within the image support and the noise standard deviation). Figure 5.1(b) shows representative images from the highest SNR image (SNR = 200) with TI = 0 andTE = 7.5 ms and the lowest SNR image (SNR = 3.83) with TI = 400 ms and TE = 217.5 ms. Chap.5 RR-CSI 75 SNR = 200 (TI = 0 ms, TE = 7.5 ms) Reference images SNR = 3.83 (TI = 400 ms, TE = 217.5 ms) (b) comp.3 2D spectrum comp.1 comp.2 Spa!al maps (a) Ground truth 30 400 T 2 (msec) 10 2 T 1 (msec) 10 3 comp.2 comp.1 comp.3 (d) Averaged 2D spectrum comp.1 comp.2 comp.3 Spa!al maps Conven!onal voxel-by-voxel es!ma!on approach (no spa!al constraint) 40 400 T 2 (msec) 10 2 T 1 (msec) 10 3 comp.2 comp.1 comp.3 Single-voxel spectrum 40 400 T 2 (msec) 10 2 T 1 (msec) 10 3 comp.2 comp.1 (c) Averaged 2D spectrum comp.1 comp.2 comp.3 Spa!al maps Our spa!al-spectral es!ma!on approach 40 400 T 2 (msec) 10 2 T 1 (msec) 10 3 comp.2 comp.1 comp.3 Single-voxel spectrum 40 400 T 2 (msec) 10 2 T 1 (msec) 10 3 comp.2 comp.1 Figure 5.1: Numerical simulation results for 2D correlation spectrum estimation. (a) Ground truth used for numerical simulation: (left) compartmental spectra f c (T 1 ;T 2 ) and (right) compart- ment spatial maps a c (x;y). (b) Representative simulated images. The highest SNR image (corre- sponding toTI = 0 ms andTE = 7.5 ms) and the lowest SNR image (corresponding toTI = 400 ms andTE = 217.5 ms) are displayed. Estimation results are shown for (c) our spatial-spectral estima- tion approach and (d) the conventional voxel-by-voxel estimation approach (no spatial constraint). Each subgure shows (left) the 2D spectrum obtained by spatially-averaging the 4D spectroscopic image, (middle) spatial maps obtained by spectrally-integrating the spectral peak locations, and (right) a 2D spectrum obtained from a single representative voxel which contains two compartments. Chap.5 RR-CSI 76 For spectroscopic image estimation, a dictionary matrix K was formed with every combination of 100T 1 values (ranging from 100 ms to 3000 ms spaced logarithmically) and 100T 2 values (ranging from 2 ms to 300 ms spaced logarithmically) for a total of Q = 10,000 dictionary elements. This type of dictionary construction is typical of previous relaxation spectroscopy methods. (Note that we also tried other dictionary constructions with dierent ranges for T 1 andT 2 , and withQ values ranging from 10,000 to 40,000. However, we only report results from this single dictionary for simplicity, because the results did not change in consequential ways when other dictionaries were used.) Optimization was performed using = 0:01, = 1, and zero initialization. The optimization was performed using in-house MATLAB software. To demonstrate the importance of spatial constraints, we also estimated 2D spectra voxel- by-voxel using conventional 2D correlation spectroscopy techniques (i.e., without the spatial con- straint). In addition, for comparison against 1D relaxometry, we simulated 1D T 1 relaxometry with the same 7 inversion times as in the high-dimensional case and 1D T 2 relaxometry with 32 echo times (TEs ranging from 10 ms to 320 ms) as in conventional approaches, with appropriate data averaging so that experiment durations were always matched. Each 1D spectroscopic image was estimated using the same basic optimization formulation from Eq. (3.12) (including spatial regularization to improve the estimation results), but modied so that we only had a 1D T 1 or T 2 spectrum at each voxel. Figures 5.1 (c) and (d) show 2D correlation spectroscopic imaging results from our spatial- spectral approach as well as conventional voxel-by-voxel 2D spectrum estimation. As can been seen in Fig. 5.1(c), two strong peaks and one weak peak are discernible in the spatially-averaged 2D T 1 -T 2 spectrum from our approach, and spatial maps obtained by spectral integration of these peaks are well matched to the ground truth spatial maps. However, the reconstructed spectral peak widths were broader than the ground-truth peaks, as should be expected based on the use of nite sampling and the resolution limits imposed by the ill-posedness of multi-exponential signal estimation [102]. In contrast, as can been seen in Fig. 5.1(d), we observe that the 2D correlation spectrum is not as well depicted when using voxel-by-voxel 2D spectrum estimation, with potentially a fourth peak emerging in between the peaks ascribed to components 1 and 2. As can be seen, the spatial maps for each component also exhibit cross-contamination, where the spatial details have bled from Chap.5 RR-CSI 77 one component to another, suggesting a lack of adequate spectral resolution. Another important observation is that our spatial-spectral estimation approach successfully resolves the ne spatial details of compartment 3, while voxel-by-voxel estimation was substan- tially less successful. This occurred despite the fact that compartment 3 is not very spatially smooth, while spatial smoothness constraints are the only dierence between the our approach and the conventional voxel-by-voxel method. These results empirically demonstrate the benets of the spatial-spectral approach to this inverse problem, and underscore the fact that the true spectroscopic image does not actually need to be very spatially smooth for these constraints to be useful. For comparison, results from the 1D relaxometry simulations are shown in Fig. 5.2. As expected based on our previous analyses, both of the 1D relaxometry approaches fail to resolve three distinct spectral peaks, and were substantially less successful than the 2D approaches at recovering the spatial maps of the original three components. T 1 (msec) 10 2 10 3 Averaged 1D T 1 spectrum comp.1 comp.2 Spa!al maps comp.1 comp.2 Spa!al maps 10 1 T 2 (msec) 10 2 Averaged 1D T 2 spectrum 10 1 T 2 (msec) 10 2 Representa!ve 1D T 2 spectra comp.1 comp.2 10 2 T 1 (msec) 10 3 Representa!ve 1D T 1 spectra comp.1 comp.2 (a) T 1 relaxometry (b) T 2 relaxometry Figure 5.2: Estimation results corresponding to simulated (a) 1D T 1 relaxometry and (b) 1D T 2 relaxometry acquisitions. Each gure shows (left) the 1D spectra obtained by spatially-averaging the 3D spectroscopic image, (middle) representative single-voxel spectra, and (right) spatial maps obtained by spectrally-integrating the spectral peak locations. Chap.5 RR-CSI 78 In this work, we used a xed value of the regularization parameter (i.e., = 0:01). This value was chosen heuristically based on the visual quality of the resulting spatial maps. However, it is important to note that dierent choices of the regularization parameter can impact the RR- CSI reconstruction result in dierent ways. For example, Fig. 5.3 shows the results of dierent choices of in the case of the higher-SNR numerical simulation. As can be seen, larger values of improve the apparent SNR of the resulting spatial maps, but are also associated with a loss in both spatial resolution and spectral resolution. While automatic \data-driven" regularization parameter selection methods have been advocated in the relaxometry literature [72, 103] and are also potentially applicable to RR-CSI, we believe that there may be strong benets to characterizing the resolution characteristics of RR-CSI using tools from the theoretical literature [91, 92], such that the trade-o between noise and resolution can be understood more directly and tuned based on the needs of the particular application (e.g., as we have done in our previous work on regularized diusion MRI [104]). This promises to be an interesting topic for further research. Chap.5 RR-CSI 79 (d) RR-CSI: λ = 1 comp.1 comp.2 comp.3 Spa"al maps Averaged 2D spectrum 40 400 T 2 (msec) 10 2 T 1 (msec) 10 3 comp.2 comp.1 comp.3 Single-voxel spectrum 40 400 T 2 (msec) 10 2 T 1 (msec) 10 3 comp.2 comp.1 (c) RR-CSI: λ = 0.1 comp.1 comp.2 comp.3 Spa"al maps Averaged 2D spectrum 40 400 T 2 (msec) 10 2 T 1 (msec) 10 3 comp.2 comp.1 comp.3 Single-voxel spectrum 40 400 T 2 (msec) 10 2 T 1 (msec) 10 3 comp.2 comp.1 (e) RR-CSI: λ = 10 comp.1 comp.2 comp.3 Spa"al maps Averaged 2D spectrum 40 400 T 2 (msec) 10 2 T 1 (msec) 10 3 comp.2 comp.1 comp.3 Single-voxel spectrum 40 400 T 2 (msec) 10 2 T 1 (msec) 10 3 comp.2 comp.1 (a) RR-CSI: λ = 0.0001 Single-voxel spectrum 40 400 T 2 (msec) 10 2 T 1 (msec) 10 3 comp.2 comp.1 Averaged 2D spectrum 40 400 T 2 (msec) 10 2 T 1 (msec) 10 3 comp.2 comp.1 comp.3 comp.1 comp.2 comp.3 Spa"al maps (b) Averaged 2D spectrum comp.1 comp.2 comp.3 Spa"al maps RR-CSI: λ = 0.001 40 400 T 2 (msec) 10 2 T 1 (msec) 10 3 comp.2 comp.1 comp.3 Single-voxel spectrum 40 400 T 2 (msec) 10 2 T 1 (msec) 10 3 comp.2 comp.1 Figure 5.3: Numerical simulation results as a function of . Each subgure shows (left) the 2D spectrum obtained by spatially-averaging the 4D spectroscopic image, (middle) spatial maps obtained by spectrally-integrating the spectral peak locations, and (right) a 2D spectrum obtained from a single representative voxel. Chap.5 RR-CSI 80 5.2.2 Pumpkin Experiments Real MRI data of a small pumpkin was acquired at every combination of 7 inversion times (TI = 0, 100, 200, 400, 700, 1000 and 2000 ms) and 15 echo times (TEs ranging from 7.5ms to 217.5 ms in 15 ms increments) for a total of P = 105 contrast encodings. We used an IR-MSE sequence on a 3T MRI scanner (Achieva; Philips Healthcare, Best, The Netherlands). For each TI, this sequence uses a train of spin-echoes to acquire data from all 15 dierent TEs in a single shot after the initial inversion recovery preparation. Acquisition used the following imaging parameters: 2 mm 2 mm in-plane resolution, 4mm slice thickness, TR = 5000 ms, a 32-channel receiver array coil, and SENSE parallel imaging with an acceleration factor of 2. A high-resolution image of the same slice was additionally acquired with 0.5 mm 0.5 mm in-plane resolution for reference. For spectroscopic image estimation, we used the same dictionary matrix K and the same optimization parameters that were used in the numerical simulation. It should be noted that TI = 0 was not practical to implement, but it would theoretically produce the same sequence of magnitude images as a standard multi-echo spin-echo sequence with- out an inversion pulse. As a result, we acquired data corresponding to TI = 0 without using an inversion pulse and manually inverted the signal polarity. This procedure assumes perfect inversion of the longitudinal magnetization. However, the rest of data was acquired with a real inversion pulse, for which inversion eciency may not be perfect due to various factors. In addition, the scanner also used a dierent (and unknown) scaling factor when saving images with an inversion pulse than it did when saving images without. As a result, the data corresponding to TI = 0 had dierent scaling compared to the rest of the data. To correct for this unknown scale factor, we t a single-compartment (monoexponential) inversion recovery curve to the data acquired with TI > 0, and used this model to synthesize what the signal should have looked like atTI = 0. We compared this synthesized data against the measured data atTI = 0 to generate a scale factor for each voxel. A single global scale correction was then obtained by averaging these voxelwise scale factors, and applied to the measured data at TI = 0. To compare our multidimensional correlation spectroscopic imaging approach against 1D methods, we also performed both 1D T 1 relaxation and 1D T 2 relaxation spectrum estimation. In particular, 1D T 1 relaxation spectra were estimated for every voxel from the seven TIs at TE = 7.5 ms, and 1D T 2 relaxation spectra were estimated for every voxel from the fteen TEs. Both Chap.5 RR-CSI 81 3D spectroscopic images (i.e., 1D spectra at every voxel) were estimated using the same estimation parameters described in the numerical simulation, including spatial regularization to improve the estimation results. (b) Representave Individual 2D spectra (c) Individual voxel locaons (a) Averaged 2D spectrum T 2 (msec) 10 1 10 2 10 2 T 1 (msec) 10 3 comp.1 comp.2 comp.3 T 1 (msec) 10 3 10 1 10 2 T 2 (msec) 10 2 comp.1 comp.2 comp.3 (d) Reference (High resoluon) comp.1 comp.2 comp.3 composite (e) Spaal maps 0 1 0 0.8 0 0.4 Figure 5.4: Estimation results from the pumpkin experiment. (a) The 2D spectrum obtained by spatially-averaging the estimated 4D spectroscopic image. (b) Representative individual spectra plotted from three dierent spatial locations. In this plot, each of the spectral peaks has been numbered and color-coded (red: component 1, green: component 2 and blue: component 3). The individual spatial locations corresponding to each of the peaks are depicted using the same color- coding scheme shown in (c) a representative image acquired at TI = 0 ms and TE = 7.5 ms. (d) A high-resolution (0.5 mm 0.5 mm) reference image. (e) Spatial maps obtained by spectrally- integrating the three spectral peaks. Each of the map is color-coded based on the color-coding scheme described in (b), with the composite image shown on the right. Figure 5.4 shows the results from our multidimensional correlation spectroscopic imaging approach. As shown in Fig. 5.4(a), we visually (subjectively) identied one strong peak and two weak peaks are resolved in the spatially-averaged 2D spectrum, and these three peaks are even more clearly distinguished when looking at the spectra from representative individual voxels as shown in Fig. 5.4(b). By spectrally-integrating these three peaks, spatial maps were generated as shown in Chap.5 RR-CSI 82 Fig. 5.4(e). The three compartments that are observed in each of the spatial maps are consistent with high-resolution features from the reference image in Fig. 5.4(d), and while we do not claim to be pumpkin experts, they appear to be consistent with our understanding of the anatomical structure of the pumpkin. For comparison, Fig. 5.5 shows results from conventional 1D methods. While the three spectral peaks were well-resolved in 2D, only two peaks are observed for both of the 1D methods. Spatial maps of these peaks further illustrate that the conventional 1D methods do not resolve compartments successfully as our multidimensional correlation spectroscopic imaging approach. comp.B comp.A T 1 (msec) 10 3 10 2 (a) T 1 relaxa!on 0 0 comp.B comp.A Spa!al maps Averaged 1D spectrum comp.D comp.C T 2 (msec) 10 2 10 1 (b) T 2 relaxa!on comp.D comp.C Spa!al maps Averaged 1D spectrum composite composite 1 0 0 0.3 1 0 0 0.8 Figure 5.5: Pumpkin experiment results from (a) conventional 1D T1 relaxometry and (b) con- ventional 1D T2 relaxometry. Each gure shows (left) the estimated spectra averaged across all voxels, and (right) spatial maps of the integrated spectral peaks. Chap.5 RR-CSI 83 5.2.3 Orange Experiments We also scanned an orange fruit using an inversion recovery turbo spin-echo imaging sequence on a 3T human MRI system (Achieva; Philips Healthcare, Best, The Netherlands). This acquisition used the following sequence parameters: TR = 33000 ms, voxel size = 1.8 mm 1.8 mm, slice thickness = 4 mm, and 33 slices. Data was sampled at every combination of 7 inversion times (TI = 0, 100, 200, 500, 1000, 1500 and 2000 ms) and 6 echo times (TE = 23, 40, 70, 100, 150 and 200 ms) for a total of P = 42 contrast encodings. Figure 5.6 shows all 42 encodings from a single slice of the orange. TI = 0 TI = 100 TI = 200 TI = 500 TI = 1000 TI = 1500 TI = 2000 TE = 23 TE = 40 TE = 70 TE = 100 TE = 150 TE = 200 unit: ms unit: ms Data set Figure 5.6: The full set (P = 7 6 = 42) of T 1 -T 2 relaxation-encoded images from a single slice of the orange. The images show the magnitude of the observed signal. Chap.5 RR-CSI 84 For this case, estimation was performed using a dictionary K that was constructed using every combination of 100 T 1 values (logarithmically distributed from 10 to 10,000 ms) and 100 T 2 values (logarithmically distributed from 2 to 10,000 ms) for a total of Q = 100 100 = 10; 000 dictionary elements. Optimization was performed using = 0:1 and zero initialization. Figure 5.7 shows results from the real data. The averaged 2D spectrum has two distinct peaks, and these peaks generate the spatial maps associated with two dierent compartments of the orange. There is no ground truth in this case, and we are not experts on the microstructural characteristics of oranges. However, the spatial maps that were identied seem to be consistent with our basic understanding of orange anatomy. T2 (ms) 10 1 10 2 T1 (ms) 10 3 10 2 (a) Orange image (b) Es!mated T1-T2 relaxa!on correla!on spectrum comp.1 comp.2 composite (c) Spa!al maps comp.1 comp.2 Figure 5.7: (a): Reference image. (b): The estimated spatially-averaged 2D T 1 -T 2 relaxation spectrum. (c): The spatial maps of the integrated spectral peaks from the estimated spectroscopic image. It is important to emphasize that the compartments that were estimated have spatial overlap, and our correlation spectroscopic imaging approach appears to be successfully resolving the asso- ciated partial volume artifacts. To illustrate this point as well as to emphasize the spatial-spectral nature of the spectroscopic image, Fig. 5.8 shows plots of the full 4D spectroscopic image from a subregion of the image (indicated with a green box in Fig. 5.7(a)). As can be seen, there are regions in which a single prominent spectral peak is present, but these are separated by a transition region in which two spectral peaks are present (likely representing partial voluming between the two compartments). Chap.5 RR-CSI 85 10 1 10 2 T2 (ms) T1 (ms) 10 3 10 2 10 1 10 2 10 3 10 2 T2 (ms) T1 (ms) 10 1 10 2 10 3 10 2 10 1 10 2 10 3 10 2 10 1 10 2 10 3 10 2 10 1 10 2 10 3 10 2 10 1 10 2 10 3 10 2 10 1 10 2 10 3 10 2 10 1 10 2 10 3 10 2 10 1 10 2 10 3 10 2 10 1 10 2 10 3 10 2 10 1 10 2 10 3 10 2 10 1 10 2 10 3 10 2 10 1 10 2 10 3 10 2 10 1 10 2 10 3 10 2 10 1 10 2 10 3 10 2 Spaally varying T1-T2 relaxaon correlaon spectra T2 (ms) T1 (ms) T2 (ms) T1 (ms) T2 (ms) T1 (ms) T2 (ms) T1 (ms) T2 (ms) T1 (ms) T2 (ms) T1 (ms) T2 (ms) T1 (ms) T2 (ms) T1 (ms) T2 (ms) T1 (ms) T2 (ms) T1 (ms) T2 (ms) T1 (ms) T2 (ms) T1 (ms) T2 (ms) T1 (ms) T2 (ms) T1 (ms) Figure 5.8: Spatially-varyingT 1 -T 2 relaxation correlation spectra from the region-of-interest cor- responding to the green box drawn in the reference image from Figure 5.7(a). Chap.5 RR-CSI 86 5.2.4 In Vivo Human Brain Experiments We also acquired in vivo human brain data using the same imaging protocol and sequence param- eters from the pumpkin experiment in Section 5.2.2. Axial slices with 4 mm thickness and coronal slices with 2 mm thickness were acquired from four healthy subjects (3 females and 1 male, and age: 30 55 years). We acquired one axial and one coronal slice from subject 1, two axial slices and two coronal slices from subject 2, and two axial slices from both subjects 3 and 4. Contrast encoding used the same 105 (TE,TI) combinations as for the pumpkin experiment. Each single-slice dataset was acquired within 20 minutes. Just like for the pumpkin experiment, scale correction was performed for the data atTI = 0, and then spectroscopic image reconstruction was performed using the same parameters described previously. To enable comparison against a 1D relaxometry method, we also used a multi-echo spin- echo sequence to acquire a T 2 relaxometry dataset from one subject with 32 TEs ranging from 10 ms to 320 ms in 10 ms increments, and otherwise using the same imaging parameters described previously. This set of sequence parameters is typical forT 2 -based myelin water imaging [8, 10], and we chose to compare against this case because multicomponent T 2 relaxometry is more common in the literature than multicomponent T 1 relaxometry. Estimation of the 3D spectroscopic image was performed using the same basic optimization formulation from Eq. (3.12) (including spatial regularization), but the spectroscopic image and the dictionary matrix were modied for 1D T 2 relaxation. For illustration, a representative single-slice 4D dataset is shown in Fig. 5.9 with correspond- ing 2D T 1 -T 2 correlation spectra shown in Fig. 5.10 and a visualization of the spectroscopic image shown in Fig. 5.11. Chap.5 RR-CSI 87 Figure 5.9: A representative single-slice 4D dataset from an in vivo human brain (the axial slice from subject 1). Chap.5 RR-CSI 88 Representave Individual 2D spectra subject 1 Individual voxel locaons Averaged 2D spectrum T 1 (msec) 10 3 10 1 10 2 T 2 (msec) 10 2 T 2 (msec) 10 1 10 2 10 2 T 1 (msec) 10 3 comp.1 comp.3 comp.2 comp.4 comp.5 comp.6 comp.1 comp.3 comp.2 comp.4 comp.5 comp.6 Figure 5.10: 2DT 1 -T 2 correlation spectra estimated from an in vivo human brain (the axial slice from subject 1). (left) The 2D spectrum obtained by spatially-averaging the estimated 4D spectro- scopic image. (middle) Representative individual spectra plotted from dierent voxels. Component 1 and component 6 are plotted from a white matter voxel, component 2 is plotted from a voxel in the putamen, component 3 is plotted from a voxel in the globus pallidus, component 4 is plotted from a voxel in gray matter, and component 5 is plotted from a voxel in the cerebral spinal uid. In this plot, each of the spectral peaks has been numbered and color-coded (red: component 1, green: component 2, cyan: component 3, blue: component 4, yellow: component 5, and magenta: component 6. This color coding scheme was also used to depict the individual voxel locations on (right) an anatomical reference image, although we do not mark the voxel for component 6 because it is the same as the voxel for component 1. As can be seen in Fig. 5.10, we visually (subjectively) identied six spectral peaks from the estimated 4D spectroscopic image. (While it may be dicult to distinguish the number of peaks that are present in the spatially-averaged spectrum, the individual peaks are actually much easier to identify when looking at the spectra from individual voxels, as can be seen in both Figs. 5.10 and 5.11.) We believe that the capability to resolve six peaks is encouraging, since conventional 1D relaxometry methods generally only resolve two or three dierent compartments in the brain. The ability of our multidimensional correlation spectroscopic imaging approach to resolve substantially more spectral peaks is consistent with our expectations about the superiority of high-dimensional encoding and spatial-spectral estimation relative to lower-dimensional approaches. Chap.5 RR-CSI 89 Figure 5.11: Visualization of the estimated 4D spectroscopic image from an in vivo human brain (the axial slice from subject 1). Spatially-varying 2D spectra are shown from (left) the entire image (81 74 voxels), (middle) a subregion of the entire image corresponding to the green box (14 14 voxels), and (right) an even smaller subregion corresponding to the orange box (4 3 voxels). In each 2D spectrum, the horizontal axis corresponds to the T 2 dimension and the vertical axis corresponds to theT 1 dimension. Each spectrum is color-coded based on the six spectral peaks and the color-coding scheme described in Fig. 5.10. As can be seen in Fig. 5.11, we also frequently observe multiple peaks coexisting within a single voxel, and the peaks each have their own distinct spatial distributions. If we ascribe the dierent spectral peaks to dierent microstructural tissue compartments, then we can interpret this 4D spectroscopic image as demonstrating the ability to resolve partial voluming and to spatially map the spatial variations of each compartment. 2D spectra obtained by spatially averaging the 4D T 1 -T 2 spectroscopic image are shown for dierent slices and subjects in Figs. 5.12 (axial slices) and 5.13 (coronal slices). The number and spectral locations of the observed spectral peaks are largely the same as observed in the spectra shown in Fig. 5.10, demonstrating that our multidimensional correlation spectroscopic imaging approach appears to yield robust and consistent results across a range of dierent subjects and slice orientations. Chap.5 RR-CSI 90 Figure5.12: Reference images (corresponding toTI = 0 andTE = 7.5 ms) and spatially-averaged 2D T 1 -T 2 spectra from dierent axial slices of dierent subjects. T 2 (msec) 10 1 10 2 10 2 Averaged 2D spectrum T 1 (msec) 10 3 Anatomy subject 1 T 2 (msec) 10 1 10 2 10 2 Averaged 2D spectrum T 1 (msec) 10 3 Anatomy subject 2 T 2 (msec) 10 1 10 2 10 2 Averaged 2D spectrum T 1 (msec) 10 3 Anatomy subject 2 (a) (b) (c) comp.1 comp.3 comp.2 comp.4 comp.5 comp.6 comp.1 comp.3 comp.2 comp.4 comp.5 comp.6 comp.1 comp.2 comp.4 comp.5 comp.6 Figure5.13: Reference images (corresponding toTI = 0 andTE = 7.5 ms) and spatially-averaged T 1 -T 2 spectra obtained from our multidimensional correlation spectroscopic imaging approach from dierent coronal slices of dierent subjects. Spatial maps obtained by spectrally integrating the six previously-identied spectral peaks are shown in Figs. 5.14 (axial slices) and 5.15 (coronal slices). We observe that these maps are also largely consistent with one another. Chap.5 RR-CSI 91 Lower slice Upper slice Figure 5.14: Spatial maps obtained from our multidimensional correlation spectroscopic imaging approach by spectrally-integrating the six spectral peaks from axial slices of dierent subjects. Each map is color-coded based on the scheme from Fig. 5.10, and the composite image is also shown on the right. Chap.5 RR-CSI 92 Figure 5.15: Spatial maps obtained from our multidimensional correlation spectroscopic imaging approach by spectrally-integrating the six spectral peaks from coronal slices of dierent subjects. Each map is color-coded based on the scheme from Fig. 5.10, and the composite image is also shown on the right. Importantly, the spatial maps also appear to qualitatively match well with known brain anatomy: component 1 seems to correspond to a white matter (WM) compartment; component 2 seems to correspond to GM structures with relatively high myelin content, including subcortical GM, putamen, thalamus and brainstem nuclei (as seen on the wall of the fourth ventricle in subject 1 in Fig. 5.15) as well as cortical GM; component 3 seems to correspond to brain structures with high iron content including the globus pallidus, subthalamic nucleus and substantia nigra; component 4 is similar to component 2, but seems to represent the GM content absent any myelin-content and notably does not include the subcortical GM; component 5 seems to correspond to cerebrospinal uid (CSF); and component 6 resembles the myelin water compartment that has been observed in previous myelin water imaging experiments [8{10]. It should be noted that the component 3 is not observed in every slice, which we believe is reasonable because the third compartment seems localized to gray matter structures like the globus pallidus, and these structures are not present in all of the slices we acquired data from. It should also be noted that the relaxation parameter values we estimated for component 5 (which seems to correspond to CSF) and component 6 (which resembles myelin water) do not match the parameter values for these tissue types reported in previous literature [5, 11, 64, 105]. Chap.5 RR-CSI 93 This is somewhat expected for a variety of reasons (including simplistic modeling assumptions as will be discussed in the next section), but is especially expected because the range of contrast en- coding parameters we used may be insucient to accurately estimate very quickly-relaxing tissues like myelin water or very slowly-relaxing tissues like CSF. Nevertheless, while the specic relax- ation parameter values we've estimated are unlikely to be accurate, we are encouraged by the fact that it appears that our multidimensional correlation spectroscopic imaging approach may still be successfully resolving the spatial maps of these components, and that the resulting maps are still consistent from slice to slice and subject to subject. Figure 5.16: Results from conventional 1D T 2 relaxation. (a) The 1D spectrum obtained by spatially-averaging the 3D spectroscopic image. (b) Representative 1D spectrum from a voxel in white matter. (c) Spatial maps obtained by spectrally-integrating the three spectral peak loca- tions. Each map is color-coded (magenta: comp.A, green: comp.B, and yellow: comp.C), and the composite image is also shown on the right. For comparison, Fig. 5.16 shows the results obtained from the 1DT 2 relaxometry experiment. As can been seen (and consistent with previous literature [8, 10]), only three spectral peaks are resolved, which is substantially fewer than the number of peaks resolved by our multidimensional correlation spectroscopic imaging approach. While the spatial maps corresponding to these three peaks all appear to be anatomically reasonable, we believe that the interpretation of these maps is less straightforward than the interpretation of the spatial maps from our approach. In particular, Chap.5 RR-CSI 94 we believe that the 1D relaxometry results are still substantially confounded by partial volume contributions, which appear to be more successfully resolved by 2D T 1 -T 2 correlation spectroscopic imaging. 5.3 Discussion This chapter described our approach to multidimensional correlation spectroscopic imaging and demonstrated empirically that this approach can have better compartmental resolving power than lower-dimensional approaches. However, using higher-dimensional contrast encoding is associated with practical increases in data acquisition time. While our human brain experiments were reason- ably fast, a 20 minute acquisition may still be too long for routine practical use of this technique. However, there are still plenty of opportunities to make the scan faster. For example, while our acquisition used a relatively large number (i.e., 105) of dierent contrast encodings, it should be noted that this number of samples was not optimized, and was chosen based on the maximum number of contrasts we could t within a 20 minute acquisition time. We have recently presented preliminary work that uses CRLB theory to optimize experimental protocols for both 2D T 1 -T 2 [41] and 2D D-T 2 [40] correlation spectroscopic imaging that will be described in Chapter 6. This preliminary work has demonstrated that we can obtain similar-quality results with substantially less than 105 encodings, which may be leveraged to enable substantial improvements in data acqui- sition time. In addition, there have been other recent constrained reconstruction approaches that have been proposed to reduce the contrast-encoding requirements of high-dimensional relaxation spectroscopy [31, 32, 77], as well as constrained reconstruction approaches that have been proposed to reduce the k-space sampling and averaging requirements of multi-contrast imaging [106{109]. Simultaneous-multislice imaging [110] could also be used to increase the volume coverage of an acquisition without increasing the acquisition time. In addition, while we relied on an IR-MSE sequence for the results shown in this paper, our approach can also be used with other pulse se- quences that may have higher eciency such as MR ngerprinting [42], inversion recovery balanced steady-state free-precession [111, 112], or a recent advanced high-dimensional contrast encoding technique [113]. Multicomponent MR ngerprinting methods have also recently been described in the literature [114, 115] that have similar objectives to 2DT 1 -T 2 correlation spectroscopic imaging, although it does not appear that these approaches were able to resolve as many tissue compartments Chap.5 RR-CSI 95 as our approach does. Recently, we have presented preliminary results of 2DT 1 -T 2 correlation spec- troscopic imaging with an MR ngerprinting acquisition with very promising results [43]. Any of these kinds of approaches could potentially be synergistically combined to make multidimensional correlation spectroscopic imaging experiments even faster. Although our experimental results appear to demonstrate successful decomposition of sub- voxel compartments from in vivo human subjects, some of the estimated relaxation parameters we estimated do not match closely with previous literature values. As explained previously, some of these discrepancies should be expected because our range of TEs and TIs may not be sucient to accurately estimate very long or very short relaxation parameters. However, it is also known that quantication of relaxation parameters can be signicantly aected by a variety of factors including experimental conditions, signal modeling, and optimization parameters, leading to mismatches even between dierent standard approaches [55]. Recognizing these issues, we have focused on qualitative evaluation in real data scenarios rather than quantitative validation because of the lack of a gold standard. Nevertheless, we were able to show promising results in terms of reproducibility and consistency throughout the experiment of multiple subjects. Although we were able to get consistently-similar and qualitatively-reasonable correlation spectroscopic imaging results, the ill-posedness of the inverse problem still means that we do not expect that our reconstructions will necessarily be the unique spectroscopic images that are con- sistent with the measured data. As a result, it would likely be fruitful in the future to do further exploration of the uncertainty of our solutions to the inverse problem, e.g., as might be achieved using optimization tools[116] or Monte Carlo methods [117, 118].) Similarly, it is also worth noting that our application examples with the IR-MSE sequence used a relatively simple acquisition and a substantially simplied physics model that does not account for the eects of water exchange, magnetization transfer, B1 inhomogeneity, B0 inhomo- geneity, slice prole eects, etc. Thus, while our current implementation appears to successfully enable separation of tissue compartments that are not easily resolved by other methods, we believe that histological validation and improving the quantitative accuracy of this kind of approach are important future objectives. (Note that it has been suggested that the standard nonnegativity constraint may be problematic for IR-MSE sequences in the presence of signicant exchange [119], such that extending our approach to account for exchange may require additional modications to Chap.5 RR-CSI 96 our constrained estimation framework.) More accurate modeling may also enable the use of more advanced pulse sequences that may be more ecient than IR-MSE sequences but are more sensitive to nonideal acquisition physics. Lastly, it should be noted that our estimation formulation from Eq. (3.12) uses a least- squares penalty to enforce consistency. The use of least-squares is appropriate under a Gaussian noise model, although the magnitude images we used in this work are more properly modeled with a Rician distribution. While the dierences between Gaussian noise and Rician noise are unlikely to matter very much in high-SNR regimes, it may be valuable to investigate the use of more accurate statistical noise modeling approaches [120] in future work. 5.4 Summary RR-CSI is a novel approach to imaging microstructure that uses high-dimensional contrast encoding ofT 1 andT 2 together with high-dimensional spatial-spectral image reconstruction. In this chapter, we described and evaluated the RR-CSI method, and demonstrated in vivo human brain results for the rst time. We believe that this approach reveals information that has not been easy acces- sible with traditional approaches, and has the potential to enable substantially more informative experiments in both basic research and clinical applications. Related Publications [Preprint] Daeun Kim, Jessica L. Wisnowski, Christopher T. Nguyen, Jusitn P. Haldar, \Relaxation-Relaxation Correlation Spectroscopic Imaging (RR-CSI): Leverag- ing the Blessings of Dimensionality to Map In Vivo Microstructure," 2019. doi: https://arxiv.org/abs/1806.05752 [Conference proceeding] Daeun Kim, Jessica L. Wisnowski, Christopher T. Nguyen, Jusitn P. Haldar, \Probing In Vivo Microstructure with T1-T2 Relaxation Correlation Spectroscopic Imaging," IEEE International Symposium on Biomedical Imaging, Wash- ington, D.C., 2018, p.675-678. doi: http://dx.doi.org/10.1109/ISBI.2018.8363664 Chap.5 RR-CSI 97 [Conference procedding] Daeun Kim, Jessica L. Wisnowski, Justin P. Haldar, \MR Corre- lation Spectroscopic Imaging of Multidimensional Exponential Decays: Probing Microstructure with Diusion and Relaxation," SPIE Wavelets and Sparsity XVI, San Diego, 2017, p.10394D. doi:http://dx.doi.org/10.1117/12.2272833 [Conference abstract] Daeun Kim, Jessica L. Wisnowski, Christopher T. Nguyen, Jusitn P. Haldar, \Multidimensional T1 Relaxation-T2 Relaxation Correlation Spectro- scopic Imaging (RR-CSI) for In Vivo Imaging of Microstructure,"International Society for Magnetic Resonance in Medicine, 27th Annual Meeting, Honolulu, 2018, p.783. Chapter 6 Toward Improved Eciency In previous chapters, we have demonstrated that our multidimensional MRI approaches are more informative than 1D experiments and potentially powerful. However, the multidimensional data acquisition can be slow due to the curse of dimensionality. In this chapter, we investigate the pos- sibility to make the experiments faster. In Section 6.1, we present an optimized experiment design approach that enables subsampling of high-dimensional contrast encoding to reduce the amount of data to acquire 1 . In Section 6.2, we investigate the potential to use a fast imaging sequence for our multidimensional approach. The summary of this chapter is provided in Section 6.3. 6.1 Optimized Experiment Design Approach Our previous implementations of DR-CSI/RR-CSI densely sampled the higher-dimensional contrast encoding space because we do not know exactly how much data is necessary to estimate high-quality high-dimensional correlation spectroscopic images. In this section, we investigate the possibility for substantially reducing the amount of data required to estimate high-quality high-dimensional correlation spectroscopic images. We specically focus an estimation-theoretic experiment design approach [121] that is similar approaches that have been used to optimize various other quantitative MRI experiments [79, 80, 82{84]. In the following sections, we describe our approach to optimized experiment design and present application examples in DR-CSI and RR-CSI with real MRI data. 1 Some of the text and gures are previously published in [41] and are copyright of the IEEE 98 Chap.6 Toward Improved Eciency 99 6.1.1 CRLB for the Multidimensional CSI Model Subsampling of 2D correlation spectroscopy experiments has been previously explored in previous literature [31, 32]. Our approach to optimal experiment design is based on the Cram er-Rao Lower bound (CRLB) [78]. The CRLB provides a lower bound on the variance of any unbiased estimator, and is widely used as a surrogate metric for optimizing experiments [79, 80, 82{84, 121]. We assume that CRLB-based experiment design enables accelerated MR-CSI acquisitions without a signicant loss of information, leading to higher quality MR-CSI results than experiments designed with some of the previously-proposed approaches for reducing the size of 2D correlation spectroscopy datasets. We rst dene the CRLB for RR-CSI, assuming that the signal observed from a single voxel withS distinct compartments. It is straight forward to use the similar derivation of the CRLB for other types of correlation spectroscopic imaging including DR-CSI once the necessary change to the signal model is made. Based on Eq. (3.6), the ideal noiseless signal observed from this voxel can be represented as: g(; p ) = S X s=1 f s 1 2e TIp=T 1s e TEp=T 2s ; (6.1) for p = 1; ;P . The vector 2R 3S1 concatenates the true spin density parameters (f 1 ;f S ), T 1 relaxation parameters (T 1 1 ; ;T 1 S ), and T 2 relaxation parameters (T 2 1 ; ;T 2 S ) for the voxel, while p = (TI p ;TE p ) is the pth set of experimental conditions. In the presence of additive white Gaussian noise with variance 2 , the entries of the Fisher information matrix F2R 3S3S can be calculated as: [F] ij = 1 2 P X p=1 @g ; p @ i @g ; p @ j ; (6.2) and the CRLB says that covariance( ^ ) F 1 (6.3) for any unbiased estimator ^ of . 6.1.2 Optimal Experiment Design The goal of statistical experiment design approaches is generally to choose experimental conditions p in such a way to make F 1 as small as possible, which has the eect of encouraging the set of Chap.6 Toward Improved Eciency 100 measurements to be as informative as possible (in the absence of additional information) for the unknown parameters of the given signal model [121]. In this work, we seek to nd experimental conditions that minimize the trace (\A-optimality") of the CRLB [121]: J f p g P p=1 ; = 3S X i=1 w i p [F 1 ] ii i ; (6.4) wherew i is a weighting coecient we use to place more emphasis on the more interesting parameters in . The optimal experiment design problem is nonconvex, and many approaches have been pro- posed in the literature to nd good local minima. For simplicity of implementation, we focus in this work on a preliminary retrospective investigation in which we explore subsampling of the full data we have already acquired, without the need to acquire additional new data. The greedy se- quential backward selection (SBS) algorithm [122] is convenient for this kind of retrospective data reduction. The SBS algorithm starts from a \full" set of candidate measurements, and sequentially discards the least informative measurement until reaching the desired number of measurements for a reduced experiment. 6.1.3 Application Examples Our experiment design was demonstrated from the previous DR-CSI experiments described Sec- tion 4.2.3 and RR-CSI experiments illustrated in Section 5.2.2, respectively. 6.1.3.1 DR-CSI Experiment We used the same ex-vivo mouse spinal cord datasets from previous DRCSI work (three sham controls and three with traumatic spinal cord injury). This data was sampled at every combination of 7 b-values (0, 500, 1000, 2000, 3000, 4000 and 5000 s/mm 2 ) and 4 TEs (40, 80, 120 and 160 ms) for a total of 28 images. Representative images are shown in Figure 6.1. The spectra estimated from this set of 28 images were considered as a fully sampled ground truth. The sampling schemes are displayed in Figure 6.1. For rectangular grid sampling, two dierent options (4 b-values x 3 TEs and 3 TEs x 4 b-values) were considered. For random sampling, ten dierent sampling realizations were considered, and we report results derived from Chap.6 Toward Improved Eciency 101 the best-performing option. For all sampling schemes, DR-CSI spectra were reconstructed using the same parameters illustrated in Section 4.2.3. Figure6.1: (a)The \fully sampled" diusion- and relaxation-encoded images from a control mouse spinal cord, containing 28 total images (every combination of 7 b-values and 4 TEs). (b) The sampling schemes we used with only 12 contrast encodings. Figure 6.2 shows spatially-averaged DR-CSI spectra for control and injured spinal cords. In the control cord, the spectra from the CRLB-based and 4x3 rectangular grid sampling schemes have the best performance in resolving the two distinct spectral peaks that are known to be present from the "fully sampled" ground truth. In the injured cord, the spectrum from the CRLB-based sampling has good consistency with the ground truth, and correctly demonstrates three distinct spectral peaks. Chap.6 Toward Improved Eciency 102 Figure 6.2: Estimated 2D diusion-relaxation correlation spectra (averaged across all voxels) from representative (a) control and (b) injured spinal cords. The rst column shows representative diusion- and relaxation-encoded images. The second column shows the integrated 2D spectra from the ground truth. The right four columns show the integrated 2D spectra from the four dierent accelerated scans, corresponding to the sampling schemes shown in Fig. 6.1. Chap.6 Toward Improved Eciency 103 Figure 6.3: (a) Spatially-varying DR-CSI spectra from (a) the \fully sampled" ground truth, (b) accelerated CRLB-based sampling, and (c) 4x3 Rectangular grid sampling. Spectra are shown from the transition region between injury and gray matter, corresponding to the red box drawn on the injured cord in Fig. 6.2. Figure 6.3(a) shows spatially-varying DR-CSI spectra from a region of the \fully sampled" ground truth in which the tissue transitions between injury and gray matter (GM). Fig. 6.3(b-c) show the DR-CSI spectra obtained from the CRLB-based and 4 3 rectangular grid accelerated acquisitions. The tissue transitions are still seen in the both of the accelerated acquisitions. Chap.6 Toward Improved Eciency 104 Figure 6.4: Spatial maps of the integrated spectral peaks from the DR-CSI reconstruction from representative (a) control and (b) injured cords. The correspondences between spatial maps and spectral peaks are indicated using the color-coding scheme shown in the left column (red: comp. 1, blue: comp. 2, green: comp. 3 and yellow: comp. 4). Figure 6.4 shows that the spatial maps of the integrated spectral peaks from the accelerated scans also have large similarity with the maps from the ground truth. In each case, the maps seem to be consistent with the known anatomy of the spinal cord: The rst three components appear to correspond with white matter and gray matter compartments. The last component is only present in the injured cords, and seems to be associated with a new compartment related to the injury. Chap.6 Toward Improved Eciency 105 6.1.3.2 RR-CSI Experiment For RR-CSI experiments, we used the same pumpkin datasets that had been used in sec. 5.2.2. The \full" dataset with P = 105 images is shown in Fig. 6.5. Our optimal experiment design used SBS to reduce the original dataset from its original size of P = 105 images down to subsets of P = 36 images (2.9 reduction) and P = 25 images (4.2 reduction). Figure 6.5: Illustration of the kind of 4D data acquisition we have previously used for RR-CSI, with 2D spatial encoding combined with 2D contrast encoding on a dense rectilinear grid. Figure 6.6 shows the specic subsets that were selected by the optimal experiment design. For each subsampled dataset, we applied RR-CSI reconstruction using the same parameters illustrated in sec. 5.2.2. Chap.6 Toward Improved Eciency 106 7.5 67.5 97.5 127.5 187.5 217.5 157.5 37.5 0 100 2000 700 400 200 1000 TE (ms) TI (ms) 7.5 67.5 97.5 127.5 187.5 217.5 157.5 37.5 0 100 2000 700 400 200 1000 TE (ms) TI (ms) 36 encodings 25 encodings Figure6.6: Subsampling patterns for (left)P = 36 contrast encodings and (right)P = 25 contrast encodings. Measured and unmeasured samples are marked with and symbols, respectively, and the dierence between P = 36 and P = 25 is additionally indicated with red color in the right image. Comparisons are shown in Fig. 6.7. We observe that the two subsampled datasets yield similar spectra and similar spatial maps compared to the original dataset shown in Fig. 6.5. This result conrms that P = 105 samples was not necessary to achieve good spectroscopic imaging results, and suggests that it may be possible to substantially improve the data acquisition speed of our previous RR-CSI implementation. Chap.6 Toward Improved Eciency 107 15 0 0 0 10 20 comp.1 comp.2 comp.3 composite (b) Spa!al maps 15 0 0 0 10 20 15 0 0 0 10 20 105 encodings 36 encodings 25 encodings T 2 (msec) T 1 (msec) 10 1 10 2 10 3 comp.2 comp.1 comp.3 Fully sampled (105 encodings) T 2 (msec) T 1 (msec) 10 1 10 2 10 3 comp.2 comp.1 comp.3 25 encodings T 2 (msec) T 1 (msec) 10 1 10 2 10 3 comp.2 comp.1 comp.3 36 encodings (a) Es!mated RR-CSI spectra Figure 6.7: RR-CSI subsampling results. (a): Averaged 2D RR-CSI spectra (summed over all voxels). (b): Spatial maps of the corresponding three peaks from (a) (red: component 1, green: component 2, blue: component 3). Chap.6 Toward Improved Eciency 108 6.2 RR-CSI with a MR Fingerprinting Acquisition Previous RR-CSI implementations in Chapter 5 used an inversion recovery multi-echo spin-echo (IR-MSE) sequence for data acquisition, and demonstrated good results with in vivo human brain data. However, RR-CSI is not restricted to an IR-MSE acquisition, and it is possible that other acquisition schemes may be more ecient than IR-MSE. In this section, we consider the use of a MR Fingerprinting (MRF) acquisition [42] for RR-CSI. MRF acquires data using dynamically varying pulse sequence parameters without a traditional steady state, and is potentially promising for RR-CSI because it is sensitive to bothT 1 andT 2 relaxation and is applicable to multicomponent relaxation modeling [114, 115]. This chapter investigates the theoretical and empirical potential of an inversion-recovery (IR) FISP MRF acquisition [123] for RR-CSI. 6.2.1 In Vivo Human Brain Experiments The IR-FISP MRF acquisition is rst compared against an IR-MSE theoretically, using the CRLB from estimation theory [78]. For the IR-MSE acquisition, we use the previous RR-CSI result obtained from human brain experiments in Section 5.2.4. The human brain data was acquired with TR=5000ms, and every combination of 7 TIs (0, 100, 200, 400, 700, 1000 and 2000 ms) and 15 TEs (from TE = 7.5 ms to 217.5 ms in 15 ms increments) for a total of 105 contrasts. For the IR-FISP MRF acquisition, we use the same \conventional MRF IR-FISP acquisition protocol from Ref. [84], which acquired data with TI = 12.34 ms, TE = 2 ms, and a variety of FAs and TRs for a total of 1500 contrasts. Both sequences are illustrated in Fig. 6.8. To ensure a fair comparison, we assume that data averaging is performed so that both sequences require the same amount of imaging time. In addition, we also acquire fully-sampled IR-FISP MRF data from an in vivo human brain on a 3T scanner, and apply the RR-CSI spectroscopic image estimation to this data. Chap.6 Toward Improved Eciency 109 Figure 6.8: (a) Illustration of the ip angles used by the IR-FISP MRF sequence as a function of time. (b) Example IR-FISP MRF signal evolutions corresponding to parameter values that are typical of white matter (T1=800ms, T2=70ms) and gray matter (T1=1000ms,T2=90ms). (c) Example IR-MSE signal evolutions for reference, using the same parameter values as in (b). Figure 6.9 shows in vivo human brain RR-CSI results from IR-FISP MRF. For reference, we also show IR-MSE results from the previous result (from a dierent subject and a dierent scanner). We observe that IR-FISP MRF resolves and generates reasonable spatial maps for 5 distinct spectral peaks. These 5 peaks appear to have good spatial correspondence with 5 of the 6 peaks obtained with IR-MSE, although the (T 1 , T 2 ) values are slightly dierent for the two acquisitions. IR-FISP MRF does not resolve the previously-identied sixth compartment, which is associated with short T 1 and may be associated with myelin water. Fig. 6.10 shows an illustration of the 4D spectroscopic image, demonstrating the resolution of partial voluming within each voxel. Chap.6 Toward Improved Eciency 110 Figure 6.9: RR-CSI results from (a) the IR-MSE dataset and (b) the IR-FISP MRF dataset. In each gure, the top row shows (left) the average 2D RR-CSI spectrum (integrated across all voxels) and (middle) representative individual spectra plotted from (right) dierent spatial locations indicated on the anatomical reference image (each of the peaks and the corresponding voxel is color- coded). The bottom row shows spatial maps of the integrated spectral peaks from the estimated spectroscopic images. Chap.6 Toward Improved Eciency 111 Figure6.10: Illustration of a small region (3 x 3 voxels) from the 4D spectroscopic image obtained from the IR-FISP MRF dataset. The spatially-varying spectra are plotted from the transition region between putamen (comp.2) and globus pallidus (comp.3) corresponding to the orange box drawn on the reference image in Fig. 6.9. Chap.6 Toward Improved Eciency 112 Figure 6.11 shows CRLB results assuming that 6 compartments are present within a voxel (the same 6 compartments identied from the IR-MSE result), while Fig. 6.12 shows CRLB results assuming that only 3 compartments are present within a voxel (based on the observed \typical" compartmental distributions from IR-MSE). We observe that both sequences have advantages rel- ative to each other in dierent contexts (lower CRLB is better), although we also observe that IR-FISP MRF has consistently worse CRLB for estimating components 3 and 6. The substantially worse CRLB for compartment 6 may explain the missing sixth spectral peak in the IR-FISP MRF data. Figure 6.11: CRLB results from the mixture of the six compartments. Each subgure shows (top) normalized CRLB values forT 1 ( p CRLB(T 1 )=T 1 ) and (bottom) normalized CRLB values for T 2 ( p CRLB(T 2 )=T 2 ). Chap.6 Toward Improved Eciency 113 Figure 6.12: CRLB results from three typical mixtures of three compartments. (a) The mixture of comp.1, comp.2 and comp.6. (b) The mixture of comp.2, comp.3 and comp.6. (c) The mix- ture of comp.3, comp.4 and comp.6. Each subgure shows (top) normalized CRLB values for T 1 ( p CRLB(T 1 )=T 1 ) and (bottom) normalized CRLB values for T 2 ( p CRLB(T 2 )=T 2 ). Chap.6 Toward Improved Eciency 114 6.3 Summary This chapter illustrated approaches that can make the high-dimensional experiments more ecient and faster. We proposed the optimized experiment design approach for the high-dimensional MRI approach and demonstrated it using DR-CSI data from rat spinal cord experiments and RR-CSI data from pumpkin experiments. We showed that the optimal experiment design approaches can allow substantial retrospective subsampling of these datasets, which can potentially enable much faster prospective experimental implementations of this approach in the future. In the context of improved eciency, we also performed RR-CSI using an IR-FISP MRF acquisition and theoretically evaluated the experimental eciency. Our results demonstrated that IR-FISP MRF can be used in the RR-CSI framework for ecient contrast encodings ofT 1 andT 2 . However, while IR-FISP MRF is viable, it is likely suboptimal, and the application of optimal experiment design approaches is likely benecial. Related Publications [Conference proceeding] Daeun Kim, Jessica L. Wisnowski, Jusitn P. Haldar, \Improved Eciency for Microstructure Imaging using High-Dimensional MR Correlation Spectroscopic Imaging," IEEE 51st Asilomar Conference on Signals, Systems, and Com- puters, Pacic Grove, 2017, p.1264-1268. doi: http://dx.doi.org/10.1109/ACSSC.2017.8335555 [Conference abstract] Daeun Kim, Bo Zhao, Lawrence L. Wald, Justin P. Haldar, \Multidi- mensional T1 Relaxation -T2 Relaxation Correlation Spectroscopic Imaging with a Magnetic Resonance Fingerprinting Acquisition," International Society for Mag- netic Resonance in Medicine, 27th Annual Meeting and Exhibition, Montreal, 2019, p4991. [Conference abstract] Daeun Kim, Justin P. Haldar, \Faster Diusion-Relaxation Cor- relation Spectroscopic Imaging (DR-CSI) using Optimized Experiment Design," International Society for Magnetic Resonance in Medicine, 25th Annual Meeting, Honolulu, 2017, p176. Chapter 7 Conclusion MRI has been recognized as one of the most important biomedical innovations over the past several decades, since it was rst developed in 1970s. It has provided new insights into the human body, revealing structural and functional information with a certain level of detail that would have not been accessible with other biomedical imaging technologies. Thanks to the decades of intense re- search and development, the MRI technology has advanced and matured to become a prevalent and powerful imaging tool that monitors biomedically-important tissue changes development, pathol- ogy, aging, etc. However, despite its advances, its capabilities have always been practically limited by undesirable trade-os between spatial resolution, signal-to-noise ratio (SNR), and data acquisi- tion time. Due to these limitations, modern human MRI experiments are typically performed with millimeter-scale voxels, even though many scientically- or clinically-interesting biological tissue features would only become directly visible at ner (e.g., microscopic or cellular) resolution scales. However, MRI-based study of tissue microstructure is still possible by leveraging the fact that cer- tain MR contrast mechanisms such as diusion and relaxation are sensitive to features of the local tissue microenvironment. This means that information about sub-voxel tissue features may still be accessible by formulating and solving an appropriate inverse problem. This dissertation has proposed a novel high-dimensional MRI approach to imaging microstruc- ture. The high-dimensional method takes advantage of using more than one MRI contrast mecha- nisms (e.g., diusion-relaxation or relaxation-relaxation). The method combines multidimensional data acquisition with spatially-constrained spectroscopic image reconstruction to substantially im- prove the ill-conditioned inverse problems that existing microstructure imaging methods have faced. 115 Chap.7 Conclusion 116 This method can be implemented in various ways considering which contrast mechanisms we select, which data acquisition sequences we use for contrast encoding, and whether 2D imaging or 3D imag- ing we want to perform. In this dissertation, we have implemented the method with two dierent specications: diusion-relaxation correlation spectroscopic imaging (DR-CSI) and T 1 relaxation- T 2 relaxation correlation spectroscopic imaging (RR-CSI). Through a wide range of experiments, we have shown that our high-dimensional methods improve the ill-posedness of the inverse problem and reduced compartmental ambiguities, compared to traditional 1D methods that use 1D contrast encoding. We have also shown that the high-dimensional methods provide high-quality spectrum estimation compared to traditional 2D methods that do not account for the spatial constraint. In this dissertation, we have also provided estimation-theoretic analysis to demonstrate bene- ts of high-dimensional contrast encoding and high-dimensional estimation that are key features of the proposed method. We have shown that the proposed method is more ecient than traditional methods in a statistical estimation theory perspective. Based on the estimation-theoretic analysis, we have further investigated two approaches to improve eciency of the method: an optimized ex- periment design approach and a fast imaging sequence approach. Our results have conrmed that multidimensional data acquisition with high-dimensional contrast encoding can be substantially accelerated without a loss of information by using either of the two approaches or both. While this dissertation has mainly focused on the proof-of-principle of the novel high-dimensional experiment frameworks, there are many possible extensions of the research presented in this dis- sertation. Although we have shown 2D contrast encoding, it could be possible to do even higher- dimensional contrast encoding. For example, N-D correlation spectroscopic imaging could be possible with diusion contrast encodings with N dierent diusion orientations and of course re- laxation contrast encodings could be combined to the N-D diusion contrast encodings. It could be also possible to use other fast imaging techniques in addition to the MRF technique [42] that we have used for RR-CSI implementation. For example, inversion recovery balanced steady-state free- precession [111], or a recent advanced high-dimensional contrast encoding technique [113] could be potentially used in the high-dimensional experiment frameworks. Any of these kinds of approaches could potentially be synergistically combined to make multidimensional correlation spectroscopic imaging experiments even faster. However, it would be also important to make experiments e- cient with these approaches. In this perspective, the optimized experiment design approach that Chap.7 Conclusion 117 we have presented would be benecial to these kinds of extensions. Microstructure imaging has become of great interest as it provides unique insights into the human body beyond the image with limited resolution. Although tremendous microstructure imag- ing techniques have been developed, however, its applications have been often limited due to the fact that most microstructure imaging techniques focus on using only one MRI contrast mechanism. This dissertation has suggested a novel high-dimensional MRI method to overcome the limitation by combining information from more than one contrast MRI mechanism. We believe the high- dimensional method, along with its future evolutions, would expand the role of MRI in probing important features of tissue microstructure that have previously been inaccessible via traditional MRI methods. Appendix A ADMM Based Solution of Eq. (3.11) In this section, we describe an ADMM based optimization algorithm to solve Eq. (3.11). First, we use variable splitting to \simplify" the optimization problem in Eq. (3.11) into the equivalent form: f ^ f i ; ^ x i ; ^ y i ; ^ z i g N i=1 = arg min ff i ;x i ;y i ;z i g N i=1 N X i=1 2 4 km i Kx i k 2 2 +I + (y i ) + X j2i kz j z i k 2 2 3 5 ; subject to x i =t i f i ; f i = y i ; f i = z i for i = 1;:::;N: (A.1) It is straightforward to derive an iterative ADMM algorithm for Eq. (A.1) following the derivations in [93]. Omitting the detailed derivations due to space constraints, we end up with the following algorithm (assuming the augmented Lagrangian parameter , whose value in uences the convergence speed of ADMM but theoretically has no in uence on the nal solution). Set iteration number j = 0, and initialize ^ f (0) i ,^ x (0) i , ^ y (0) i , and ^ z (0) i to arbitrary values. Also initialize corresponding vectors of Lagrange parameters g (0) xi , g (0) yi , and g (0) zi to arbitrary values. At the (j + 1)th iteration of ADMM: 118 Chap.7 Conclusion 119 1. For each voxel i = 1;:::;N, update the estimate of ^ f i according to ^ f (j+1) i = arg min f i kt i f i (^ x (j) i + g (j) xi )k 2 2 +kf i (^ y (j) i + g (j) yi )k 2 2 +kf i (^ z (j) i + g (j) zi )k 2 2 = 8 < : 1 3 h ^ x (j) i + g (j) xi + ^ y (j) i + g (j) yi + ^ z (j) i + g (j) zi i ; if t i = 1 1 2 h ^ x (j) i + g (j) xi + ^ y (j) i + g (j) yi i ; if t i = 0: (A.2) 2. For each voxel i = 1;:::;N, update the estimate of ^ x i according to ^ x (j+1) i = arg min x i km i Kx i k 2 2 +kx i (t i ^ f (j+1) i g (j) xi )k 2 2 = K H K + I 1 K H m i + (t i ^ f (j+1) i g (j) xi ) : (A.3) Note that the matrix being inverted is not very large and is always the same for every i and j, so the inverse matrix can simply be precomputed once and reused at every iteration. 3. For each voxel i = 1;:::;N, update the estimate of ^ y i according to ^ y (j+1) i = arg min y i ky i ( ^ f (j+1) i g (j) yi )k 2 2 +I + (y i ): (A.4) This solution is obtained by rst simply setting ^ y (j+1) i = ^ f (j+1) i g (j) yi , and then setting any negative entries of ^ y (j+1) i equal to zero. 4. Jointly update the estimates off ^ z i g N i=1 according to f^ z (j+1) i g N i=1 = arg min f^ z i g N i=1 N X i=1 2 4 kz i ( ^ f (j+1) i g (j) zi )k 2 2 + X j2i kz j z i k 2 2 3 5 : (A.5) This optimization problem may look complicated, but simplies in two separate ways. First, Eq. (A.5) is decoupled along the spectral dimension. If we let [z i ] q denote theqth Chap.7 Conclusion 120 entry of the length-Q vector z i , then the solution to Eq. (A.5) is equivalent to setting f[^ z (j+1) i ] q g N i=1 = arg min f^ z i g N i=1 N X i=1 2 4 k[z i ] q ([ ^ f (j+1) i ] q [g (j) zi ] q )k 2 2 + X j2i k[z j ] q [z i ] q k 2 2 3 5 (A.6) for q = 1;:::;Q. Second, the least-squares system matrices corresponding to Eq. (A.6) will have block circulant structure if the neighborhood system i is constructed using periodic boundary conditions, and thus can be inverted noniteratively using standard fast Fourier transform methods [93]. 5. For each voxel i = 1;:::;N, update the Lagrange parameters by setting g (j+1) xi = g (j) xi (t i ^ f (j+1) i ^ x (j+1) i ); g (j+1) yi = g (j) yi ( ^ f (j+1) i ^ y (j+1) i ); and g (j+1) zi = g (j) zi ( ^ f (j+1) i ^ z (j+1) i ): (A.7) 6. Set j =j + 1. Iterate steps 1-6 until convergence. Appendix B ADMM Based Solution of Eq. (3.12) In this section, we describe an aADMM based optimization algorithm to solve Eq. (3.12). Compared to the algorithm described in Appendix A (which used vector-representations of the data and the spectroscopic image within each computation step), the description below instead relies on matrix representations of the variables in Eq. (3.11) as in Eq. (3.12). This is not just a notational dierence, because computing each step using the matrix representation in Eq. (3.12) removes the need for for-loops and reduces computational complexity. The algorithm proceeds as follows: Set iteration number j = 0, and initialize ^ F (0) 2 R QN , ^ X (0) 2 R QN , ^ Y (0) 2 R QN , ^ Z (0) 2 R QN , G (0) 2 R QN , H (0) 2 R QN and R (0) 2 R QN to arbitrary values. The variables X, Y, and Z are used for variable splitting, while the variables G, H, and R correspond to Lagrange multipliers. Also choose an augmented Lagrangian parameter value > 0 (this choice in uences the convergence speed of ADMM but not the solution). At iteration (j + 1): 1. For each i = 1;:::;N, update the estimate of ^ f i (i.e., the columns of ^ F) according to ^ f (j+1) i = 8 > > > > > > < > > > > > > : 1 3 ^ x (j) i + g (j) i + ^ y (j) i + h (j) i + ^ z (j) i + r (j) i ; if t i = 1 1 2 ^ y (j) i + h (j) i + ^ z (j) i + r (j) i ; if t i = 0; (B.1) 121 Chap.7 Conclusion 122 where ^ x (j) i , g (j) i , ^ y (j) i , h (j) i , ^ z (j) i and r (j) i are respectively the ith columns of ^ X (j) , G (j) , ^ Y (j) , H (j) , ^ Z (j) , and R (j) . 2. Set ^ X (j+1) = K H K +I 1 K H M + ^ F (j+1) T G (j) ; (B.2) where I denotes the identity matrix. Note that the matrix K H K +I is small and its inverse can easily be precomputed and stored for use in every iteration. 3. Set ^ Y (j+1) = ^ F (j+1) H (j) , but replacing any negative values with zero. 4. Set ^ Z (j+1) = ^ F (j+1) R (j) I +C H C 1 : (B.3) Note that if we assume periodic boundary conditions for the spatial smoothness regular- ization, then the matrices I and C H C can both be diagonalized by the discrete Fourier transform, and we can use standard Fourier methods to quickly and analytically compute the desired matrix inversion result [? ]. 5. Set G (j+1) = G (j) ( ^ F (j+1) T ^ X (j+1) ); H (j+1) = H (j) ( ^ F (j+1) ^ Y (j+1) ); and R (j+1) = R (j) ( ^ F (j+1) ^ Z (j+1) ): (B.4) 6. Set j =j + 1. Iterate steps 1-6 until convergence. Bibliography [1] Lauterbur PC. Image formation by induced local interactions: Examples employing nuclear magnetic resonance. Nature 1973;242:190{191. [2] Jones DK, editor. Diusion MRI: Theory, Methods, and Applications. Oxford: Oxford Uni- versity Press. 2011. [3] Tournier JD, Mori S, Leemans A. Diusion tensor imaging and beyond. Magn Reson Med 2011;65:1532{1556. [4] Le Bihan D, Johansen-Berg H. Diusion MRI at 25 exploring brain tissue structure and function. NeuroImage 2012;61:324{341. [5] MacKay A, Laule C, Vavasour I, Bjarnason T, Kolind S, M adler B. Insights into brain mi- crostructure from the T2 distribution. Magn Reson Imag 2006;24:515{525. [6] Fenrich F, Beaulieu C, Allen P. Relaxation times and microstructures. NMR Biomed 2001; 14:133{139. [7] Deoni SC. Quantitative relaxometry of the brain. Top Magn Reson Imaging 2010;21:101{113. [8] MacKay A, Whittall K, Adler J, Li D, Paty D, Graeb D. In vivo visualization of myelin water in brain by magnetic resonance. Magn Reson Med 1994;31:673{677. [9] Labadie C, Lee JH, Vetek G, Springer CS. Relaxographic imaging. J Magn Reson B 1994; 105:99 { 112. [10] Whittall KP, MacKay AL, Graeb DA, Nugent RA, Li DKB, Paty DW. In vivo measurement of T2 distributions and water contents in normal human brain. Magn Reson Med 1997;37:34{43. 123 Bibliography 124 [11] Does MD, Beaulieu C, Allen PS, Snyder RE. Multi-component T1 relaxation and magneti- sation transfer in peripheral nerve. Magn Reson Med 1998;16:1033{1041. [12] Deoni SC, Rutt BK, Arun T, Pierpaoli C, Jones DK. Gleaning multicomponent T1 and T2 information from steady-state imaging data. Magn Reson Med 2008;60:1372{1387. [13] Niendorf T, Dijkhuizen RM, Norris DG, van Lookeren Campagne M, Nicolay K. Biexponential diusion attenuation in various states of brain tissue: Implications for diusion-weighted imaging. Magn Reson Med 1996;36:847{857. [14] Mulkern RV, Gudbjartsson H, Westin CF, Zengingonul HP, Gartner W, Guttmann CRG, Robertson RL, Kyriakos W, Schwartz R, Holtzman D, Jolesz FA, Maier SE. Multi-component apparent diusion coecients in human brain. NMR Biomed 1999;12:51{62. [15] Assaf Y, Freidlin RZ, Rohde GK, Basser PJ. New modeling and experimental framework to characterize hindered and restricted water diusion in brain white matter. Magn Reson Med 2004;52:965{978. [16] Assaf Y, Blumenfeld-Katzir T, Yovel Y, Basser PJ. Axcaliber: A method for measuring axon diameter distribution from diusion MRI. Magn Reson Med 2008;59:1347{1354. [17] Wang Y, Wang Q, Haldar JP, Yeh FC, Xie M, Sun P, Tu TW, Trinkaus K, Klein RS, Cross AH, Song SK. Quantication of increased cellularity during in ammatory demyelination. Brain 2011;134:3590{3601. [18] Zhang H, Schneider T, Wheeler-Kingshott CA, Alexander DC. NODDI: Practical in vivo neurite orientation dispersion and density imaging of the human brain. NeuroImage 2012; 61:1000{1016. [19] Scherrer B, Schwartzman A, Taquet M, Sahin M, Prabhu SP, Wareld SK. Characterizing brain tissue by assessment of the distribution of anisotropic microstructural environments in diusion-compartment imaging (DIAMOND). Magn Reson Med 2015;76:963{977. [20] Kaden E, Kelm ND, Carson RP, Does MD, Alexander DC. Multi-compartment microscopic diusion imaging. NeuroImage 2016;139:346 { 359. Bibliography 125 [21] Istratov AA, Vyvenko OF. Exponential analysis in physical phenomena. Rev Sci Instrum 1999;70:1233{1257. [22] H urlimann MD, Venkataramanan L. Quantitative measurement of two-dimensional distribu- tion functions of diusion and relaxation in grossly inhomogeneous elds. J Magn Reson 2002; 157:31{42. [23] Callaghan PT, Godefroy S, Ryland BN. Use of the second dimension in PGSE NMR studies of porous media. Magn Reson Imag 2003;21:243{248. [24] Galvosas P, Callaghan PT. Multi-dimensional inverse Laplace spectroscopy in the NMR of porous media. C R Physique 2010;11:172{180. [25] Freed DE, H urlimann MD. One- and two-dimensional spin correlation of complex uids and the relation to uid composition. C R Physique 2010;11:181{191. [26] Mitchell J, Chandrasekera T, Gladden L. Numerical estimation of relaxation and diusion distributions in two dimensions. Prog NMR Spect 2012;62:34{50. [27] Celik H, Bouhrara M, Reiter D, Fishbein K, Spencer R. Stabilization of the inverse Laplace transform of multiexponential decay through introduction of a second dimension. J Magn Reson 2013;236:134{139. [28] Cloninger A, Czaja W, Bai R, Basser PJ. Solving 2D Fredholm integral from incomplete measurements using compressive sensing. SIAM J Imaging Sci 2014;7:1775{1798. [29] Ahola S, Zhivonitko VV, Mankinen O, Zhang G, Kantola AM, Chen HY, Hilty C, Koptyug IV, Telkki VV. Ultrafast multidimensional Laplace NMR for a rapid and sensitive chemical analysis. Nat Commun 2015;6:8363. [30] Hatka A, Celik H, Cloninger A, Czaja W, Spencer RG. 2D sparse sampling algorithm for ND Fredholm equations with applications to NMR relaxometry. In: International Conference on Sampling Theory and Applications (SampTA). 2015; pp. 367{371. [31] Bai R, Cloninger A, Czaja W, Basser PJ. Ecient 2D MRI relaxometry using compressed sensing. J Magn Reson 2015;255:88{99. Bibliography 126 [32] Benjamini D, Basser PJ. Use of marginal distributions constrained optimization (MADCO) for accelerated 2D MRI relaxometry and diusometry. J Magn Reson 2016;271:40{45. [33] Lauterbur PC, Kramer DM, House Jr WV, Chen CN. Zeugmatographic high resolution nu- clear magnetic resonance spectroscopy. Images of chemical inhomogeneity within macroscopic objects. J Am Chem Soc 1975;97:6866{6868. [34] Kim D, Kim JH, Haldar JP. Diusion-relaxation correlation spectroscopic imaging (DR-CSI): An enhanced approach to imaging microstructure. In: Proc. Int. Soc. Magn. Reson. Med.. 2016; p. 660. [35] Kim D, Doyle EK, Wisnowski JL, Kim JH, Haldar JP. Diusion-relaxation correlation spec- troscopic imaging: A multidimensional approach for probing microstructure. Magn Reson Med 2017;78:2236{2249. [36] Kim D, Wisnowski JL, Haldar JP. MR correlation spectroscopic imaging of multidimensional exponential decays: probing microstructure with diusion and relaxation. In: SPIE Wavelets and Sparsity XVII. volume 10394. 2017; p. 103940D. [37] Kim D, Wisnowski JL, Nguyen CT. Probing in vivo microstructure with t1-t2 relaxation correlation spectroscopic imaging. In: Proc. IEEE Int. Symp. Biomed. Imag.. 2018; pp. 675{ 678. [38] Kim D, Wisnowski JL, Nguyen CT, Haldar JP. Multidimensional t1 relaxation-t2 relaxation correlation spectroscopic imaging (RR-CSI) for in vivo imaging of microstructure. In: Proc. Int. Soc. Magn. Reson. Med.. 2018; p. 783. [39] Kim D, Wisnowski JL, Nguyen CT. Relaxation-relaxation correlation spectroscopic imaging (rr-csi): Leveraging the blessings of dimensionality to map in vivo microstructure. 2019. priprint, arXiv:1806.05752. [40] Kim D, Haldar JP. Faster diusion-relaxation correlation spectroscopic imaging (DR-CSI) using optimized experiment design. In: Proc. Int. Soc. Magn. Reson. Med.. 2017; p. 176. Bibliography 127 [41] Kim D, Wisnowski JL, Haldar JP. Improved eciency for microstructure imaging using high- dimensional MR correlation spectroscopic imaging. In: Asilomar Conference on Signals, Sys- tems, and Computers. 2017; . [42] Ma D, Gulani V, Seiberlich N, Liu K, Sunshine JL, Duerk JL, Griswold MA. Magnetic resonance ngerprinting. Nature 2013;495:187{192. [43] Kim D, Bo Z, Wald LL, Haldar JP. Multidimensional T1 relaxation -T2 relaxation correlation spectroscopic imaging with a magnetic resonance ngerprinting acquisition. In: Proc. Int. Soc. Magn. Reson. Med.. 2019; p. 4991. [44] Hanson LG. Is quantum mechanics necessary for understanding magnetic resonance? Con- cepts in Magnetic Resonance Part A 2008;32A:329{340. [45] Nishimura DG. Principles of magnetic resonance imaging. Stanford University. 1996. [46] Liang ZP, Lauterbur PC. Principles of magnetic resonance imaging: a signal processing per- spective. SPIE Optical Engineering Press. 2000. [47] Haacke EM, Brown RW, Thompson MR, Venkatesan R. Magnetic resonance imaging: physical principles and sequence design. volume 82. Wiley-Liss New York:. 1999. [48] Bernstein MA, King KF, Zhou XJ. Handbook of MRI pulse sequences. Burlington: Elsevier Academic Press. 2004. [49] Abragam A, Abragam A. The principles of nuclear magnetism. 32. Oxford university press. 1961. [50] Day P. The philosopher's tree: a selection of Michael Faraday's writings. CRC Press. 1999. [51] Hoult D, Lauterbur PC. The sensitivity of the zeugmatographic experiment involving human samples. J Magn Reson 1979;34:425{433. [52] Twieg DB. The k-trajectory formulation of the nmr imaging process with applications in analysis and synthesis of imaging methods. Med Phys 1983;10:610{621. [53] Hahn EL. Spin echoes. Phys Rev 1950;80:580. Bibliography 128 [54] Carr HY. Steady-state free precession in nuclear magnetic resonance. Phys Rev 1958; 112:1693{1701. [55] Stikov N, Boudreau M, Levesque IR, Tardif CL, Barral JK, Pike GB. On the accuracy of T1 mapping: Searching for common ground. Magn Reson Med 2014;73:514{522. [56] Drain L. A direct method of measuring nuclear spin-lattice relaxation times 1949;62:301. [57] Hahn EL. An accurate nuclear magnetic resonance method for measuring spin-lattice relax- ation times. Phys Rev 1949;76:145. [58] Stejskal EO, Tanner JE. Spin diusion measurements: spin echoes in the presence of a time- dependent eld gradient. J Chem Phys 1965;42:288{292. [59] Callaghan PT. Principles of nuclear magnetic resonance microscopy. Oxford University Press on Demand. 1993. [60] Le Bihan D, Warach SJ. Diusion and perfusion magnetic resonance imaging: Applications to functional MRI. 1995. [61] Edelstein WA, Mahesh M, Carrino JA. MRI: Time is dose{and money and versatility 2010; 7:650{652. [62] Kroeker RM, Henkelman RM. Analysis of biological NMR relaxation data with continuous distributions of relaxation times. J Magn Reson 1986;69:218{235. [63] Whittall KP, MacKay AL. Quantitative interpretation of NMR relaxation data. J Magn Reson 1989;84:134{152. [64] Does MD, Gore JC. Compartmental study of T1 and T2 in rat brain and trigeminal nerve in vivo. Magn Reson Med 2002;47:274{283. [65] Lancaster JL, Andrews T, Hardies LJ, Dodd S, Fox PT. Three-pool model of white matter. J Magn Reson 2003;17:1{10. [66] Du YP, Chu R, Hwang D, Brown MS, Kleinschmidt-DeMasters BK, Singel D, Simon JH. Fast multislice mapping of the myelin water fraction using multicompartment analysis of T2* decay at 3T: a preliminary postmortem study. Magn Reson Med 2007;58:865{870. Bibliography 129 [67] Yablonskiy DA, Bretthorst GL, Ackerman JJ. Statistical model for diusion attenuated MR signal. Magn Reson Med 2003;50:664{669. [68] Prony R. Essai experimental et analytique: sur les lois de la dilatabilite des uides elastiques et sur celles de la force expansive de la vapeur de l'alkool, a dierentes temperatures. J L'Ecole Polytechnique 1795;1:24{76. [69] Peemoeller H, Shenoy R, Pintar M. Two-dimensional NMR time evolution correlation spec- troscopy in wet lysozyme. J Magn Reson 1981;45:193{204. [70] English A, Whittall K, Joy M, Henkelman R. Quantitative two-dimensional time correlation relaxometry. Magn Reson Med 1991;22:425{434. [71] Saab G, Thompson RT, Marsh GD, Picot PA, Moran GR. Two-dimensional time correlation relaxometry of skeletal muscle in vivo at 3 tesla. Magn Reson Med 2001;46:1093{1098. [72] Song Y, Venkataramanan L, H urlimann M, Flaum M, Frulla P, Straley C. T1-T2 correlation spectra obtained using a fast two-dimensional Laplace inversion. J Magn Reson 2002;154:261{ 268. [73] Bernin D, Topgaard D. NMR diusion and relaxation correlation methods: New insights in heterogeneous materials. Curr Opin Colloid Interface Sci 2013;18:166{172. [74] Slator PJ, Hutter J, Palombo M, Jackson LH, Ho A, Panagiotaki E, Chappell LC, Rutherford MA, Hajnal JV, Alexander DC. Combined diusion-relaxometry MRI to identify dysfunction in the human placenta. Magn Reson Med 2019;82:95{106. [75] Mulkern RV, Haker SJ, Maier SE. On high b diusion imaging in the human brain: rumina- tions and experimental insights. Magn Reson Imag 2009;27:1151{1162. [76] Lawson CL, Hanson RJ. Solving least squares problems. SIAM. 1995. [77] Benjamini D, Basser PJ. Magnetic resonance microdynamic imaging reveals distinct tissue microenvironments. NeuroImage 2017;163:183{196. [78] Kay SM. Fundamentals of statistical signal processing. Prentice Hall PTR. 1993. Bibliography 130 [79] Jones J, Hodgkinson P, Barker A, Hore P. Optimal sampling strategies for the measurement of spin{spin relaxation times. J Magn Reson 1996;113:25{34. [80] Cavassila S, Deval S, Huegen C, Van Ormondt D, Graveron-Demilly D. Cram er{rao bounds: an evaluation tool for quantitation. NMR Biomed 2001;14:278{283. [81] Ogg RJ, Kingsley PB. Optimized precision of inversion-recovery T 1 measurements for con- strained scan time. Magn Reson Med 2004;51:625{630. [82] Pineda AR, Reeder SB, Wen Z, Pelc NJ. Cram er-rao bounds for three-point decomposition of water and fat. Magn Reson Med 2005;54:625{635. [83] Alexander DC. A general framework for experiment design in diusion mri and its application in measuring direct tissue-microstructure features. Magn Reson Med 2008;60:439{448. [84] Zhao B, Haldar JP, Liao C, Ma D, Griswold MA, Setsompop K, Wald LL. Optimal experiment design for magnetic resonance ngerprinting: Cram er-rao bound meets spin dynamics. IEEE Trans Med Imag 2019;38:844{861. [85] Teixeira RPAG, Malik SJ, Hajnal JV. Joint system relaxometry (JSR) and Cramer-Rao lower bound optimization of sequence parameters: A framework for enhanced precision of DESPOT T 1 and T 2 estimation. Magn Reson Med 2018;79:234{245. [86] Lin Y, Haldar JP, Li Q, Conti PS, Leahy RM. Sparsity constrained mixture modeling for the estimation of kinetic parameters in dynamic PET. IEEE Trans Med Imag 2014;33:173{185. [87] Kumar D, Nguyen TD, Gauthier SA, Raj A. Bayesian algorithm using spatial priors for multiexponential T2 relaxometry from multiecho spin echo MRI. Magn Reson Med 2012; 68:1536{1543. [88] Labadie C, Lee JH, Rooney WD, Jarchow S, Aubert-Fr econ M, Springer CS, Mller HE. Myelin water mapping by spatially regularized longitudinal relaxographic imaging at high magnetic elds. Magn Reson Med 2013;71:375{387. [89] Bertero M, Boccacci P. Introduction to Inverse Problems in Imaging. London: Institute of Physics Publishing. 1998. Bibliography 131 [90] Haldar JP, Liang ZP. On MR experiment design with quadratic regularization. In: Proc. IEEE Int. Symp. Biomed. Imag.. IEEE. Proc. IEEE Int. Symp. Biomed. Imag.. 2011; pp. 1676{1679. [91] Fessler JA, Rogers WL. Spatial resolution properties of penalized-likelihood image recon- struction: Space-invariant tomographs. IEEE Trans Image Process 1996;5:1346{1358. [92] Ahn S, Leahy RM. Analysis of resolution and noise properties of nonquadratically regularized image reconstruction methods for PET. IEEE Trans Med Imag 2008;27:413{424. [93] Afonso MV, Bioucas-Dias JM, Figueiredo MA. An augmented Lagrangian approach to the constrained optimization formulation of imaging inverse problems. IEEE Trans Image Process 2011;20:681{695. [94] Gatidis S, Schmidt H, Martirosian P, Schwenzer NF. Development of an MRI phantom for diusion-weighted imaging with independent adjustment of apparent diusion coecient val- ues and T2 relaxation times. Magn Reson Med 2013;72:459{463. [95] Hwang D, Du YP. Improved myelin water quantication using spatially regularized non- negative least squares algorithm. J Magn Reson Imag 2009;30:203{208. [96] Kim JH, Tu TW, Bayly PV, Song SK. Impact speed does not determine severity of spinal cord injury in mice with xed impact displacement. J Neurotraum 2009;26:1395{1404. [97] Tu TW, Kim JH, Wang J, Song SK. Full tensor diusion imaging is not required to assess the white-matter integrity in mouse contusion spinal cord injury. J Neurotraum 2010;27:253{262. [98] Sodickson DK, Feng L, Knoll F, Cloos M, Ben-Eliezer N, Axel L, Chandarana H, Block T, Otazo R. The rapid imaging renaissance: sparser samples, denser dimensions, and glimmer- ings of a grand unied tomography. In: Medical Imaging 2015: Biomedical Applications in Molecular, Structural, and Functional Imaging, Gimi B, Molthen RC, editors. SPIE-Intl Soc Optical Eng. 2015; . [99] Meiboom S, Gill D. Modied spin-echo method for measuring nuclear relaxation times. Rev Sci Instrum 1958;29:688. Bibliography 132 [100] McGivney DF, Deshmane A, Jiang Y, Ma D, Griswold MA. The partial volume problem in MR ngerprinting from a bayesian perspective. In: Proc. Int. Soc. Magn. Reson. Med.. 2016; p. 435. [101] Hamilton JI, D A, Griswold M, Seiberlich N. MR ngerprinting with chemical exchange (MRF-X) for in vivo multi-compartment relaxation and exchange rate mapping. In: Proc. Int. Soc. Magn. Reson. Med.. 2016; p. 431. [102] Bertero M, Boccacci P, Pike ER. On the recovery and resolution of exponential relaxation rates from experimental data: a singular-value analysis of the Laplace transform inversion in the presence of noise. Proc R Soc Lond A 1982;383:15{29. [103] Zou YL, Xie RH, Arad A. Numerical estimation of the choice of the regularization parameter for NMR T2 inversion. Pet Sci 2016;13:237{246. [104] Haldar JP, Wedeen VJ, Nezamzadeh M, Dai G, Weiner MW, Schu N, Liang ZP. Im- proved diusion imaging through SNR-enhancing joint reconstruction. Magn Reson Med 2013; 69:277{289. [105] Does MD. Inferring brain tissue composition and microstructure via MR relaxometry. Neu- roImage 2018;. [106] Haldar JP, Liang ZP. Joint reconstruction of noisy high-resolution MR image sequences. In: Proc. IEEE Int. Symp. Biomed. Imag.. 2008; pp. 752{755. [107] Gong E, Huang F, Ying K, Wu W, Wang S, Yuan C. PROMISE: Parallel-imaging and compressed-sensing reconstruction of multicontrast imaging using sharable information. Magn Reson Med 2015;73:523{535. [108] Zhao B, Lu W, Hitchens TK, Lam F, Ho C, Liang ZP. Accelerated MR parameter mapping with low-rank and sparsity constraints. Magn Reson Med 2015;74:489{498. [109] Bilgic B, Kim TH, Liao C, Manhard MK, Wald LL, Haldar JP, Setsompop K. Improving parallel imaging by jointly reconstructing multi-contrast data. Magn Reson Med 2018;80:619{ 632. Bibliography 133 [110] Barth M, Breuer F, Koopmans PJ, Norris DG, Poser BA. Simultaneous multislice (SMS) imaging techniques. Magn Reson Med 2016;75:63{81. [111] Schmitt P, Griswold MA, Jakob PM, Kotas M, Gulani V, Flentje M, Haase A. Inversion recovery TrueFISP: Quantication of T1, T2, and spin density. Magn Reson Med 2004;51:661{ 667. [112] Pster J, Breuer FA, Jakob PM, Blaimer M. Fast multi-component T1 and T2 correlation measurements using steady-state free precession. In: Proc. Int. Soc. Magn. Reson. Med.. Proc. Int. Soc. Magn. Reson. Med.. 2018; p. 220. [113] Hutter J, Slator PJ, Christiaens D, Teixeira RPA, Roberts T, Jackson L, Price AN, Malik S, Hajnal JV. Integrated and ecient diusion-relaxometry using ZEBRA. Scientic Reports 2018;8:15138. [114] McGivney D, Deshmane A, Jiang Y, Ma D, Badve C, Sloan A, Gulani V, Griswold M. Bayesian estimation of multicomponent relaxation parameters in magnetic resonance nger- printing. Magn Reson Med 2018;80:159{170. [115] Tang S, Fernandez-Granda C, Lannuzel S, Bernstein B, Lattanzi R, Cloos M, Knoll F, Asslnder J. Multicompartment magnetic resonance ngerprinting. priprint, arXiv:1802.10492. [116] Parker RL, Song YQ. Assigning uncertainties in the inversion of NMR relaxation data. J Magn Reson 2005;174:314{324. [117] Prange M, Song YQ. Quantifying uncertainty in NMR spectra using monte carlo inversion. J Magn Reson 2009;196:54{60. [118] Prange M, Song YQ. Understanding NMR spectral uncertainty. J Magn Reson 2010;204:118{ 123. [119] Dortch RD, Horch RA, Does MD. Development, simulation, and validation of NMR relaxation-based exchange measurements. J Chem Phys 2009;131:164502. [120] Varadarajan D, Haldar JP. A majorize-minimize framework for Rician and non-central chi MR images. IEEE Trans Med Imag 2015;34:2191{2202. Bibliography 134 [121] Pukelsheim F. Optimal design of experiments. SIAM. 2006. [122] Reeves SJ, Zhe Z. Sequential algorithms for observation selection. IEEE Trans Signal Process 1999;47:123{132. [123] Jiang Y, Ma D, Seiberlich N, Gulani V, Griswold MA. Mr ngerprinting using fast imaging with steady state precession (FISP) with spiral readout. Magn Reson Med 2015;74:1621{1631.

Abstract (if available)

Abstract MRI is a powerful and widely used biomedical imaging technique, but it suffers from limited spatial resolution. This prevents MRI from being used to directly probe microscopic tissue information. Fortunately, indirect estimation of microscopic tissue features is still possible by modeling each voxel as a mixture of different microstructural components with distinct signal characteristics, and encoding microstructural information by manipulating MRI physics. However, the conventional indirect approaches are often associated with a highly ill-posed inverse problem, which can severely limit the ability to accurately unmix the microstructural components. This dissertation proposes and evaluates a new high-dimensional MRI approach that combines multidimensional spatial encoding with multidimensional contrast encoding (e.g., simultaneous diffusion encoding and relaxation encoding) and multidimensional regularized estimation to provide substantially improved microstructural unmixing capabilities. ❧ The new approach can be justified theoretically, and this dissertation demonstrates an estimation theoretic analysis that shows the multidimensional approach has substantial advantages over conventional lower-dimensional approaches. This estimation theoretic framework is also shown to enable optimized experiment design. The dissertation also evaluates the new approach in two different configurations (one using diffusion-relaxation contrast encoding and the other using relaxation-relaxation contrast encoding) in a wide range of experiments, including numerical simulations and real experiments with phantoms, ex vivo animals, and in vivo humans. The results demonstrate that the proposed approach provides unique information that is not available with conventional methods, and suggests that the proposed approach has the potential to have a high impact if translated to practical applications in medicine, neuroscience and biology.

Linked assets

University of Southern California Dissertations and Theses

Conceptually similar

PDF

Shift-invariant autoregressive reconstruction for MRI

PDF

Correction, coregistration and connectivity analysis of multi-contrast brain MRI

PDF

Fast flexible dynamic three-dimensional magnetic resonance imaging

PDF

Improved brain dynamic contrast enhanced MRI using model-based reconstruction

PDF

Estimating liver iron non-invasively with high-field MRI

PDF

Functional real-time MRI of the upper airway

PDF

Improving sensitivity and spatial coverage of myocardial arterial spin labeling

PDF

Optimization methods and algorithms for constrained magnetic resonance imaging

PDF

Novel theoretical characterization and optimization of experimental efficiency for diffusion MRI

PDF

New theory and methods for accelerated MRI reconstruction

PDF

Dynamic cardiovascular magnetic resonance imaging for improved assessment of ischemic heart disease

PDF

Improving the sensitivity and spatial coverage of cardiac arterial spin labeling for assessment of coronary artery disease

PDF

Fast upper airway MRI of speech

PDF

Mapping water exchange rate across the blood-brain barrier

PDF

Measuring functional connectivity of the brain

PDF

Seeing sleep: real-time MRI methods for the evaluation of sleep apnea

PDF

Quantitative MRI for oxygenation imaging in cerebrovascular diseases

PDF

Characterization of lenticulostriate arteries using high-resolution black blood MRI as an early imaging biomarker for vascular cognitive impairment and dementia

PDF

Dynamic functional magnetic resonance imaging to localize activated neurons

PDF

Unveiling the white matter microstructure in 22q11.2 deletion syndrome with diffusion magnetic resonance imaging

Asset Metadata

Creator Kim, Daeun (author)

Core Title High-dimensional magnetic resonance imaging of microstructure

School Viterbi School of Engineering

Degree Doctor of Philosophy

Degree Program Electrical Engineering

Publication Date 08/05/2019

Defense Date 06/04/2018

Publisher University of Southern California (original), University of Southern California. Libraries (digital)

Tag constrained reconstruction,correlation spectrum,magnetic resonance imaging,microstructure,multicomponent modeling,OAI-PMH Harvest,relaxometry and diffusometry

Format application/pdf (imt)

Language English

Contributor Electronically uploaded by the author (provenance)

Advisor Haldar, Justin P. (committee chair), Nayak, Krishna S. (committee member), Wood, John C. (committee member)

Creator Email daeunk@usc.edu,dankim108@gmail.com

Permanent Link (DOI) https://doi.org/10.25549/usctheses-c89-206711

Unique identifier UC11663241

Identifier etd-KimDaeun-7700.pdf (filename),usctheses-c89-206711 (legacy record id)

Legacy Identifier etd-KimDaeun-7700.pdf

Dmrecord 206711

Document Type Dissertation

Format application/pdf (imt)

Rights Kim, Daeun

Type texts

Source University of Southern California (contributing entity), University of Southern California Dissertations and Theses (collection)

Access Conditions The author retains rights to his/her dissertation, thesis or other graduate work according to U.S. copyright law. Electronic access is being provided by the USC Libraries in agreement with the a...

Repository Name University of Southern California Digital Library

Repository Location USC Digital Library, University of Southern California, University Park Campus MC 2810, 3434 South Grand Avenue, 2nd Floor, Los Angeles, California 90089-2810, USA