University of Southern California Dissertations and Theses

Theory and simulation of diffusion magnetic resonance imaging on brain's white matter

(USC Thesis Other)

Theory and simulation of diffusion magnetic resonance imaging on brain's white matter

PDF

Download a page range

Download transcript

Copy asset link

Request this asset

Transcript (if available)

Content THEORY AND SIMULATION OF DIFFUSION MAGNETIC RESONANCE IMAGING ON BRAIN'S WHITE MATTER by Shahryar Karimi-Ashtiani A Dissertation Presented to the FACULTY OF THE USC GRADUATE SCHOOL UNIVERSITY OF SOUTHERN CALIFORNIA In Partial Fulllment of the Requirements for the Degree DOCTOR OF PHILOSOPHY (ELECTRICAL ENGINEERING) August 2010 Copyright 2010 Shahryar Karimi-Ashtiani Dedication To my lovely wife Shirin, and to my parents ii Acknowledgements This work would not have been possible without the help of many family members, friends, advisors, and collaborators. During the time of my PhD which was quite an era, I have enjoyed the help and blessing of many people and my enormous debt and gratitude can hardly be repaid to them. I would like to thank my PhD dissertation committee members and my advisor Professor Jay Kuo, at USC's Viterbi School of Engineering, for their guidance. I am truly grateful for all moral and medical supports I received from Dr. and Mrs. Craig Morris and treating me like a family member which result in severe impact on my health condition improvement. There is no doubt in my mind that I would not have come this far without having them on my side. My thanks go to Dr. Antonio Desalles at UCLA neurosurgery for making great improvement in my health condition, making the continuation of my work possible. My thanks go out to Dr. Majeed Ahmadi for his ongoing pure advices, inspirations, and unreserved helps. I also wish to thank my grandparents Mrs. and Mr. Nikzad. My thoughts go out to my late grandfather Mr. Ali Karimi-Ashtiani whose blessing touched my life in very profound ways. iii I would like to thank my brother and sister, Shahram and Shahrzad. I am grateful for my wife, Dr. Shirin Ebrahimi for her love and sacrices she made towards my well being and surpassing dicult life's hurdles. Finally, I would like to thank my mother Parvaneh Nikzad, and specially my father, Mr. Hossein Karimi-Ashtiani for being such magnicent parents. I always enjoy their unconditional love and supports. Shahryar Karimi-Ashtiani Summer 2010 iv Table of Contents Dedication ii Acknowledgements iii List of Figures vii Abstract ix Chapter 1 : Introduction 1 1.1 Signicance of the Research . . . . . . . . . . . . . . . . . . . . . . 1 1.2 Brain Imaging Techniques and Open Issues . . . . . . . . . . . . . . 3 1.2.1 Diusion Magnetic Resonance Imaging (D-MRI) . . . . . . . 3 1.2.2 Diusion Tensor Imaging (DTI) . . . . . . . . . . . . . . . . 4 1.2.3 High Angular Resolution Diusion Imaging (HARDI) . . . . 5 1.2.4 Challenges and Open Issues . . . . . . . . . . . . . . . . . . 5 1.3 Review of Previous Research . . . . . . . . . . . . . . . . . . . . . . 7 1.4 Contributions of the Research . . . . . . . . . . . . . . . . . . . . . 8 1.5 Organization of Dissertation . . . . . . . . . . . . . . . . . . . . . . 9 Chapter 2 : Research Background 10 2.1 Central Nervous System . . . . . . . . . . . . . . . . . . . . . . . . 10 2.2 Diusion Magnetic Resonance Imaging . . . . . . . . . . . . . . . . 15 2.3 Physical Modeling of Self-diusion Process . . . . . . . . . . . . . . 19 2.4 Variational Formulation of Self-diusion PDE . . . . . . . . . . . . 21 2.5 Numerical Solution of Self-diusion PDE . . . . . . . . . . . . . . . 25 2.5.1 The Galerkin Method and Method of Lines . . . . . . . . . . 25 2.5.2 Numerical Solution of ODEs . . . . . . . . . . . . . . . . . . 26 2.5.3 Finite Element Methods . . . . . . . . . . . . . . . . . . . . 31 Chapter 3 : Anisotropic Self-diusion Simulation 35 3.1 Introduction . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 35 3.2 Theory and Methods . . . . . . . . . . . . . . . . . . . . . . . . . . 40 3.2.1 Variational Formulation and FEM . . . . . . . . . . . . . . . 45 3.2.2 Validation of the Solution . . . . . . . . . . . . . . . . . . . 53 v 3.2.3 Membrane Wall Eect . . . . . . . . . . . . . . . . . . . . . 54 3.2.4 Voxel Average Propagator . . . . . . . . . . . . . . . . . . . 56 3.3 The Finite Element Method (FEM) . . . . . . . . . . . . . . . . . . 58 3.4 The Minimum Rule . . . . . . . . . . . . . . . . . . . . . . . . . . . 61 3.5 Variational Formulation . . . . . . . . . . . . . . . . . . . . . . . . 62 3.6 Initial Unrestricted Self-diusion Regime . . . . . . . . . . . . . . . 65 3.7 Results . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 66 3.7.1 Simulation Parameters . . . . . . . . . . . . . . . . . . . . . 66 3.7.2 Contributions of the FEM Basis Functions in the Global Ap- proximate Solution . . . . . . . . . . . . . . . . . . . . . . . 67 3.7.3 Equivalent membrane model . . . . . . . . . . . . . . . . . . 68 3.7.4 Hexagonal Bundles of Healthy Axons . . . . . . . . . . . . . 70 3.7.5 Voxel Average Propagators for Single Tract Scenarios . . . . 72 3.7.6 Equivalent Tensor Model for Propagator . . . . . . . . . . . 75 3.7.7 Voxel Average Propagator for Two Crossing Tracts . . . . . 80 3.7.8 Demyelinated Tracts . . . . . . . . . . . . . . . . . . . . . . 81 3.8 Discussion . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 84 3.9 Conclusion . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 87 Chapter 4 : Theory for Practical MRI and Diusion MRI Simulation 89 4.1 Physics of MRI and Diusion MRI . . . . . . . . . . . . . . . . . . 89 4.2 Generic Imaging Framework in Practical MRI/Diusion MRI . . . . 102 4.2.1 Practical MRI Signal Generation . . . . . . . . . . . . . . . 102 4.2.2 Practical Diusion MRI Signal Generation . . . . . . . . . . 111 4.3 Conclusion and Major Contributions . . . . . . . . . . . . . . . . . 122 Chapter 5 : Conclusion and Future Work 125 5.1 Conclusion . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 125 5.2 Future Work . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 127 Bibliography 129 vi List of Figures 2.1 : Illustration of the human central nervous system . . . . . . . 11 2.2 : Illustration of the action potential . . . . . . . . . . . . . . . 11 2.3 : Illustration of a neuron cell . . . . . . . . . . . . . . . . . . . 12 2.4 : The cross section of an axon and the myelin sheath . . . . . 13 2.5 : The white and the gray matters in the human brain . . . . . 14 2.6 : White tracts . . . . . . . . . . . . . . . . . . . . . . . . . . . 15 2.7 : Motion of water molecules in axons . . . . . . . . . . . . . . 16 2.8 : 3D Elements . . . . . . . . . . . . . . . . . . . . . . . . . . . 33 3.1 : Axons model . . . . . . . . . . . . . . . . . . . . . . . . . . 45 3.2 : Tetrahedralization . . . . . . . . . . . . . . . . . . . . . . . . 48 3.3 : Membrane equivalent model . . . . . . . . . . . . . . . . . . 55 3.4 : Tetrahedral element . . . . . . . . . . . . . . . . . . . . . . . 59 3.5 : Three compartment Environment . . . . . . . . . . . . . . . 66 3.6 : Distribution of solution energy . . . . . . . . . . . . . . . . . 69 3.7 : Taner-Powel equivalent models . . . . . . . . . . . . . . . . . 70 3.8 : Transverse propagator . . . . . . . . . . . . . . . . . . . . . 71 3.9 : Parallel and transverse propagators . . . . . . . . . . . . . . 73 vii 3.10 : Validation of the FEM solution . . . . . . . . . . . . . . . . 74 3.11 : Average propagator . . . . . . . . . . . . . . . . . . . . . . . 76 3.12 : Time evolution of propagator for the single tract case . . . . 77 3.13 : Equivalent tensor model . . . . . . . . . . . . . . . . . . . . 80 3.14 : Voxel average propagator for two crossing tracts . . . . . . . 81 3.15 : Eects of demyelination on propagator . . . . . . . . . . . . 82 3.16 : Time evolution of propagator for the demyelination case . . 83 3.17 : Decomposing components . . . . . . . . . . . . . . . . . . 87 4.1 : A spin magnet and its magnetic eld . . . . . . . . . . . . . 90 4.2 : The coordinate system adopted in an MRI machine . . . . . 90 4.3 : Spin precession in the presence of an external magnetic eld 91 4.4 : Parallel and transverse components of local magnetization . 92 4.5 : RF pulse application . . . . . . . . . . . . . . . . . . . . . . 93 4.6 : De-phasing . . . . . . . . . . . . . . . . . . . . . . . . . . . . 95 4.7 : A linear gradient along the x-axis . . . . . . . . . . . . . . . 97 4.8 : The spin echo pulse sequence . . . . . . . . . . . . . . . . . 98 4.9 : The spin echo pulse sequence with diusion gradients . . . . 101 4.10 : 1-D MR image reconstruction . . . . . . . . . . . . . . . . . 109 4.11 : Traveling spin magnet . . . . . . . . . . . . . . . . . . . . . 113 4.12 : A spin traveling in dierent materials . . . . . . . . . . . . . 116 4.13 : Relative error . . . . . . . . . . . . . . . . . . . . . . . . . . 122 4.14 : Relative error as a function of T 2 . . . . . . . . . . . . . . . 123 viii Abstract Diusion MRI (D-MRI) has opened a new front for uncovering the convoluted structure of the central nervous system by providing capability for the non-invasive identication of white matter tract geometries in the brain. One of the major open issues in fully extending this technology to the clinical domain is the lack of in vivo validation of the results, which makes it dicult to have objective comparisons of dierent D-MRI techniques. To this end, the application of simulated data from known ground truths appears to be the second best choice. It is well understood that, this imaging modality is characterized by the shape of the self-diusion (SD) prole within the brain bers. The previous methods for the quantication of the SD process in the white matter environments suer from the lack of enough gen- erality or solution precision, and often impose excessive computational complexity, which limits their full extent of applicability. The contributions of this research are two folds: 1) the development of a new numerical method to compute the self- diusion SD prole and 2) the provision of a generic framework for reconstruction of diusion MR images under imperfect imaging conditions. In the rst part of this thesis, a numerical paradigm based on the nite element methods (FEM) for nding the solution of SD partial dierential equation (PDE) ix in the white matter environment is proposed. The standard FEM is incapable of accommodating the boundary conditions of the PDE on the interfaces of dierent white matter materials. Theoretical constraints to modify the FEM are investi- gated such that it becomes applicable to the problem of interest. Our method is not conned to any special geometry and virtually can handle any microscopic models of white matters. One of the highlights of the developed method is address- ing the eect of partial permeability of the cell membrane. Also, a self-validation technique is proposed to guarantee the accuracy of the solution. In the meantime, the aggregate SD prole of the MRI voxel is analytically derived to show the de- pendency of the macroscopic behavior of the diusion at the macroscopic level, as a function of the contained tissue microscopic parameters. Several simulation results are provided to showcase the eectiveness of our method. In the second part of this thesis, a generic theoretical framework for reconstruc- tion of diusion MR images, which factors most of the imaging artifacts due to true values of imaging conditions, is developed. It is worthwhile to point out that, on the simulation front, the previous D-MRI reconstruction formulations were derived under ideal imaging conditions. However, in reality, the actual values of imaging parameters leave the available reconstruction techniques inapplicable. The severe impacts of those parameters on the quality of the white matter mapping from the D-MRI data necessitate the extension of existing methods to include them. For example, the existing methods cannot accommodate spatially varying T1 and T2 parameters, a condition very common in the white matter tissues. We provide an x analytical derivation for such conditions, and show that, under these situations, the k-space signal and the eective spin magnet density will no longer form a Fourier transform pair. Hence, the arising question is the interpretation of the image under spatially varying relaxation parameters, which is the time average of the spin den- sity during the sampling period. Besides, under long data acquisition periods, the spin magnets experience dierent material by having diusion motion, which urges to nding eective values for those parameters. We provide a formulation for those quantities as functions of the parameters of the composing materials. In reality, the D-MRI data acquisitions are performed with a long diusion gradient, where spin phase encodings also occur during the wide pulses. This requires the quanti- cation of the spin magnet motions in the white matter environment. This problem is addressed analytically, and it is shown that the change in the reconstruction formulation appears in the form of modication in the q-space vector. xi Chapter 1 : Introduction 1.1 Signicance of the Research The central nervous system (CNS) plays a signicant role in dierent low and high level aspects of human life. It contributes substantially to maintaining a vast range of vital activities of dierent organs and provides a fast global communicating tools for message passings and coordinations between dierent parts of the body. The extent of the CNS in uence ranges from subconscious activities such as generation and controlling signals for heart beats, and responses to outside stimuli to more perceptual high level and deliberate functions such as thinking, walking, etc. Considering the signicance of the CNS, failure of the system could impose dramatic consequences which make the continuation of the subject's life dicult or even impos- sible. Especially in the past decades, as a downside of industrialization and growth of large societies, physical and psychological interactions among people have signicantly increased, which cause neurological and neuropsychological diseases to soar noticeably. Also, regional and international military con icts leave societies with many victims of 1 wars, who suer from a variety of medical problems and injuries which urge to nd better treatments to enhance the quality of their lives. Naturally, the rst step towards nding treatments and cures for the disorders is to have thorough understanding of the CNS anatomies and functions. In the past two cen- turies, as a result of the endeavors of neuroscientists, valuable structures and operations of the CNS have been unraveled. However, due to the complexity of the system, which includes microscopic processes, such as the molecular interactions of neuron synapses, as well as the macroscopic functions, a lot of questions remain unanswered, and many issues are still open-ended. This is a serious problem when we consider the anatomy of the CNS, specially in the brain. The CNS is subject-specic and there are many across- subject variations. For example, trajectories of the white tracts in the brain undergo variations from one person to another. In this research, we attempt to develop a software tool for precise computation of the self-diusion process in the white matter tracts of the human brain. We would like to handle the general white matter tract and axon conditions including the complex geome- tries observed in real ber tracts. The solution of the partial dierential equation (PDE) governing the self-diusion process in 3D, is explored using the nite element method (FEM), which is believed to be among the best numerical methods in terms of the prob- lem formulation and the solution accuracy. Our study is expected to provide a numerical foundation for the future exploration of the relation between the microscopic axons pa- rameters and the average voxel self-diusion prole. A prominent feature of this work is that the spatial and temporal evolutions of the self-diusion prole are quantitatively provided, as opposed to the parametric or the non-parametric ad-hoc methods such as 2 the tensor model that have been proposed before. In addition, to simulate more realistic data, a theoretical framework, where MRI and diusion MRI formulations are extended to include real world imaging parameters, is provided. 1.2 Brain Imaging Techniques and Open Issues 1.2.1 Diusion Magnetic Resonance Imaging (D-MRI) The diculty of the CNS study necessitates the employment of multidisciplinary ap- proaches. The invention of the Magnetic Resonance Imaging (MRI) [16, 31] in late 80s was a manifestation of such perception. During the past three decades, the applications of many engineering, physics, and mathematical methodologies and tools in medicine have become prevalent, which set up a new ground for the treatment of the CNS abnormalities. In early 90s, the development of Diusion MRI (D-MRI) opened a new front for the assessment of the brain white matter tract anatomy and better understanding of the connectivity of dierent regions of the brain. This led to procedures to demonstrate neural tracts in the CNS, which are referred to as tractography. It remains the only technology known for the non-invasive in vivo study of the geometry of the white matter tracts in the brain and the spinal cord. The D-MRI tractography is characterized by the simple principle that the self-diusion of water molecules in well-organized tissues has strong directional preference along the tissues orientation. Therefore, by measuring the direction of the self-diusion process, it is possible to estimate the orientation of the underlying tissue in the 3D space. The D-MRI is the only technology available for such measurements. 3 Since the invention of D-MRI, this imaging modality has found many applications in the diagnoses and studies of many diseases and abnormalities associated with the CNS [47]. Examples include brain ischemia [38, 75, 74], autism [3, 63, 5], oncology [55, 42, 69], epilepsy [32, 35], multiple sclerosis [26, 17, 4], Parkinson's [76], and Alzheimer's [12, 54] diseases, and Huntington disease [53]. D-MRI has also been employed for studies of psychiatric disorders and diseases such as schizophrenia [52, 36, 40], and behavior disorders [1, 39]. 1.2.2 Diusion Tensor Imaging (DTI) In early 90s, the invention of Diusion Tensor Imaging (DTI) [7, 10, 6] was a major breakthrough to the understanding of geometries of white tracts, which was characterized by modeling the self-diusion within the MRI voxels by the tensor model. The simplicity of parameter estimation in this model and the relatively short data acquisition time in MRI scanners, made DTI the most popular method in the eld of tractography, which retains its signicance to the present date. However, this method is prone to several shortcomings which severely limit its practical applications. The tensor model is not capable of capturing the properties of the self-diusion prole under the partial volume condition (i.e. when two or more white tracts cross each other in a voxel). Besides, it is dicult to relate microscopic parameters of white tracts (e.g., axon diameters and myelin sheath thicknesses) to macroscopic integral self-diusion of the MRI voxel. 4 1.2.3 High Angular Resolution Diusion Imaging (HARDI) In the last decade, the family of High Angular Resolution Diusion Imaging (HARDI) [70, 71, 49], has received a lot of attention for surpassing the partial volume problem imposed by DTI, in order to unveil the structure of the white tracts beyond the voxel resolution of MRI. In these techniques, imaging is performed with diusion gradients in many directions to have better assessments of the 3D prole of the self-diusion process. However, this could substantially increase the imaging time, which causes a number of undesired practical and economical consequences. The HARDI method, does not provide any link between the microstructure of axons and their aggregate eects on the average propagator of the voxel, either. 1.2.4 Challenges and Open Issues Realization of the white matter tract micro-geometry in a voxel from nite resolution D-MRI signals remains an open issue. The signicance of this problem becomes more transparent since a number of degenerations and nerve abnormalities, such as axon de- myelination, occur in the microscopic level. For accurate diagnostic purposes, it is crucial to address this issue in most general cases. The natural direction to follow to answer this complex problem is to study the eect of the white tracts microstructure composition on the self-diusion prole and D-MRI signals, with the perception of leveraging the lessons learned from this study to the reverse problem of assessments of the micro-geometry of axons from the D-MRI data. However, the large degree of variability of axon constituting and organizational parameters makes it dicult to perform such studies for general cases. In the literature, there are a number 5 of analytical and numerical attempts in this regard, but a majority of them consider only simplied models of axons, which limit the extents of their applicabilities to practical and clinical problems [56, 11, 59, 61, 51, 62, 27, 28, 34]. There is considerable variability in anatomical and structural morphology of human brains. This includes inconsistency in sizes, shapes and locations of dierent regions of the brain across dierent subjects. This makes it extremely dicult and inaccurate to establish statistical models for the study of the brain and the diagnoses of abnormalities in individual subjects. As a result of large variations in dierent brains, any attempt for setting up a statistical model would lead to a large model variances which increases the degree of uncertainty in decision making according to the model. This is also valid for the DMRI tractography. We observe a lot of dierences between the results of two subjects. Thus, tractography results should be interpreted on the subject by subject basis. This raises the serious question of in vivo validation of tractography results. That is, a D-MRI tractography algorithm might manifest a number of the white tract trajectories of the subject. The question is that whether they all represent real white tracts and, if yes, how precise the identied tracts are with respect to the ground truth. On the other hand, it is desired to be assured that all existing white tracts in the region of interest are identied. The answers to these questions for in vivo tissues are challenging, since it is hardly possible to identify the ground truth of white tracts for an alive subject with the current available technology. It might be possible to validate results by scanning terminal patients and perform postmortem studies, [72]. However, this is often a tedious, expensive and erroneous procedure. Besides, results of such studies, which are usually performed for special cases, are dicult to generalize. 6 1.3 Review of Previous Research As an alternative validation solution, application ofD-MRI phantoms has become popular in recent years [8, 45]. The general methodology is to simulate the D-MRI signal for known geometries of white tracts and apply the data as the gold standard for validation of dierent tractography methods [2, 64, 15, 18, 25]. All previous methods for the simulation of D-MRI signals are built on the assumption that the self-diusion proles in MRI voxels are based on the Gaussian distribution characterized by the voxel tensor model or the mixture of tensors for the case of several tracts crossing in a voxel [67, 9, 50, 30, 41, 68]. Despite the simplicity of such an approach, there are a number of issues associated with them. For example, it is dicult to nd relations between parameters of the white tract microstructure and the tensor characteristics. To evaluate the performance of tractography, it is of special interest to experiment the eects of these parameters on the result. Generally speaking, a conclusive evaluation should include the assessments of tractography outcomes under dierent tract conditions. Even for the simple case of a single white tract, the tensor model is decient as a generative model. Also, it is crucial that the simulated signal resembles real D-MRI data, which requires the generative algorithm to take into account the underlying physical equations governing the self-diusion process and D-MRI signal generation. Obviously, such capacity is not available with the tensor model. Furthermore, it has been shown that the performance of tractography is strongly impacted by dierent D-MRI artifacts and noise [43, 20, 44]. The origin of many of these phenomena is in the physical processes 7 involved in theD-MRI signal formation which cannot be well characterized by the tensor model. 1.4 Contributions of the Research The following specic contributions are made in this research. We develop a software tool for precise computation of the self-diusion prole in which many constraints imposed by the previous work have been relaxed. The implemented method provides the capacity to handle general white tract and axon conditions including complex geometries observed in real ber tracts. The numerical solution of the partial dierential equation (PDE), which governs the self-diusion process in 3D, is performed using the nite element method (FEM). The FEM is among the best numerical methods in problem formulation and the so- lution accuracy for boundaries of complex geometries. This framework is expected to provide a numerical foundation for future exploration of the relation between mi- croscopic axon parameters and the average voxel self-diusion prole. A prominent feature of our work is that the spatial and temporal evolutions of the self-diusion prole can be quantitatively provided. We provide extensions to MRI andD-MRI signal formulations by including realistic imaging and intrinsic parameters. Most existing equations are only valid under ideal or near ideal conditions. Our new derivation incorporates more challenging and realistic situations, e.g., the eects of the T 1 and T 2 relaxations and wide diusion gradients for D-MRI data acquisition. 8 1.5 Organization of Dissertation The rest of this dissertation is organized as follow. In Chapter 2, some research back- ground information is provided. Specically, we show that the D-MRI signal formation equations are characterized by the self-diusion prole of the voxels. In Chapter 3, the FEM formulation of the self-diusion problem in a multi-compartment environment is presented. The eects of dierent white tract components on the trend of propagation of water molecules as well as the average prole of the voxel under unusual tract conditions such as ber crossings and demyelination are illustrated. In Chapter 4, theoretical exten- sions of MRI and D-MRI formulations are provided in order to include a larger range of parameters involved in image reconstruction. In Chapter 5, concluding remarks are given and possible future extensions are pointed out. Since our method provides the precise values of the prole for the voxel, it is possible to account for all the above mentioned requirements in the simulation of fair gold standards for validation of tractography. 9 Chapter 2 : Research Background 2.1 Central Nervous System The human Central Nervous System (CNS) consists of the brain and the spinal cord. As illustrated in Fig. 2.1, the communication between the CNS and other parts of the human body is through peripheral nerves. The CNS interferes in and controls many vital activities of the body and serves as a fast communicating tool between dierent organs. There are many electric and chemical subsystems in the CNS which make the task of message passing possible. The message passing throughout the CNS is performed via electric impulsive signals which are referred to as action potentials. Fig. 2.2 shows a sample action potential. Usually a message includes a train of action potentials. For example, when the in- tention is to ex a muscle, the region of the brain in charge of that muscle generates a sequence of action potentials which are carried through the CNS transmission lines to the muscle. The more action potentials are in the train, the stronger the muscle exes. 10 Figure 2.1: Illustration of the human central nervous system:http://www.shands.org Figure 2.2: Illustration of the action potential: http://piratesandrevolutionaries 11 In order to generate and transmit such electric signals for communication, the CNS bio- logical structure has been adapted to special anatomical characteristics to accommodate such necessities. Neuron cells are the constituting elements of the CNS. Their special shapes provide them with the capability to receive action potentials from other neurons or sensors and amplify and transmit them to subsequent neurons or muscles. A neuron is illustrated in Fig. 2.3. Figure 2.3: Illustration of a neuron cell: http://www.enchantedlearning.com A neuron anatomy consists of three major components: dendrites, cell body, and axon. The role of the cell body, which includes the cell nucleus, is to keep the cell alive. The dendrites act as receivers and the neuron acquires electric messages from the upwind nervous system, which are either the axons of other neurons or the sensory signals of receptor cells (in the end nerve cases). It is evident from the gure that a neuron can receive the action potential from a number of dendrites. A neuron could possess up to millions of dendrites, which shows the complexity of the signaling activities in the CNS. The axon serves as the output mechanism to transmit the action potential to the downwind nervous system, which could include dendrites of other neuron or muscles (in the case of end nerves). It is worth mentioning that the relation between the received 12 and output action potentials of a neuron is often more complicated than a simple linear equation and there are many microscopic and molecular processes involved. For example, it is well understood that neurotransmiters play a substantial role in dening the input- output relation of a neuron. The mono-polar structure of a neuron makes it possible to have a reliable electric communication path over the CNS. The conveyance of messages to distant locations in the CNS is achieved through long axons which act analogously to long electric wires. This could range from the microscopic distance between two adjacent neurons to a very large distance from the thalamus in the brain to the feet. Similar to wires, axons are surrounded by isolating layers to prevent the action potential charges from leaking as they travel along axons. These isolating layers are referred to as myelin sheaths, which are of protein material and whitish color. In Fig. 2.4 the cross section of a single axon and the myelin sheath layers are illustrated. As can be seen from the gure, the myelin sheath thickness could be quite large, which gives this region a comparable volume fraction to the axon region. It is a good electric isolator as well as a good water insulator. Figure 2.4: The cross section of an axon and the myelin sheath Often, axons of similar functions are grouped together to form bundles of nerve bers which are referred to as the white tracts (WT), which is a justied term since the color 13 of the tracts is actually white due to the white color of myelin sheaths. This is contrary to the gray matter of the CNS which are mostly of cell bodies and dendrites. From the macroscopic point of view, a major dierence between the white and the gray matters of the CNS lie in their structural organizations. The tissues belonging to the white matters often have a well-parallel organization, while such a structural order is rarely seen in the gray matters. By the macroscopic point of view, the scale of the MRI voxels is in the order of millimeters. Fig. 2.5 shows the white and gray matter regions of the brain. Figure 2.5: The white and the gray matters in the human brain :http://www.cic- caracas.org The gray matter region includes the surface layer of the brain sphere which is referred to as the cerebral cortex. Most of the brain high-level/concous processing occurs in dierent subregions of the cortex. The communication between the regions of the cortex as well as the regions of the cortex and other parts of the body is achieved through WTs. In Fig. 2.6, two dierent WTs are illustrated. They are the visual tracts, which covey the visual signals between the eyes and the visual cortex of the brain, and the lateral tract of the corticospinal tract, which connects the motor cortex region to the spinal cord. They provide ne motor control of limbs. 14 Figure 2.6: White tracts: The lateral part of corticospinal tract (left) and the visual tract (right): http://www.lrn.org 2.2 Diusion Magnetic Resonance Imaging Diusion Magnetic Resonance Imaging (D-MRI) is the only available imaging technology for non-invasive uncovering the microstructure of human brain's WTs. It is character- ized by detection of dierent tendencies in the self-diusive Brownian motion of water molecules in dierent spatial orientations. Therefor, if the underlying tissue within an MRI voxel has a well-organized microstructure which causes water molecules to have directional preference to move along the tissue orientation, tracing the molecule motion could unveil the structural orientation of the tissue. Such tissues can be observed in the white matters of the brain, the spinal cord and in muscles. In the CNS white mat- ters, the partial permeability of axon myelin sheaths hinders the penetration of water molecules in transverse directions to the WTs orientation and, consequently, the average long term self-diusion of molecules has directional preference along the orientation of WTs as illustrated in Fig. 2.7. 15 Figure 2.7: Motion of water molecules in axons: A cross section of motion of a water molecule within an axon, where the myelin sheath directs the motion along the axon orientation. The long-term self-diusion time is characterized by the required time during which the microstructure of axons is well-probed by the water molecules. Typical values for such a time period is in the range of 20ms in real MRI scanners. As a result, the self- diusion of water molecules in WTs is anisotropic, which makes them a good candidate to be studied with the D-MRI technology. The signal in a D-MRI scanner is formed according to the following two equations [16]. First, we have the k-space signal given by S(k;q)/ Z R 0 (r 0 )e 2ik:r 0 Z R P r 0(r; )e 2iq:(r 0 r) drdr 0 ; (2.1) where S is the k-space signal, k is a sample extracted from the k-space, q is a q-space sample which could be proportional to the applied diusion gradient if the applied gra- dient is constant as the time evolves, (r 0 ) is the local density of water molecules, the P r 0(r; ) is the self-diusion propagator. Then, the D-MRI and the k-space signals are 16 related together through a Fourier transform pair. That is, the point-wise D-MRI signal observed at point r 0 as a result of an applied diusion gradient can be rewritten as E (r 0 ;q)/F[S(k;q)]; (2.2) whereF denotes the Fourier transform operation. The above equations concern the ideal signal acquisition scenario where the resolution of an image is assumed to be innite. In practice, the signal formation is performed in a voxel-wise fashion where the observed signal is the average representation of the point-wise signal values over the voxel. This is mostly due to the limited scanning time and the eect of the T 2 weighting, which cause imperfect k-space sampling. As a result of limited scan time, the k-space sampling is often truncated which causes blurring in the output D-MRI data. In the course of data acquisition in realD-MRI, all identities characterizing Eqs. (2.1)- (2.2) except for the self-diusion propagator, P r 0(r; ), are known. It is characterized by the the probability of displacement of a water molecule from the position r 0 at time t = 0 to the position r at time t = . It is worthwhile to point out that the eect of nonhomogeneous molecules density (r 0 ) is alleviated by normalization against the zero diusion gradient image. Discovery of the structure of WTs in the MRI voxel from a limited amount of D- MRI data necessitates the models for the self-diusion propagator. This leads to the emergence of several parametric and non-parametric models to describe the self-diusion behavior in WTs. The tensor model found a lot of popularity for the characterization of the self-diusion prole in the WTs. The underlying assumption for the validity of such a 17 model is that the water molecules self-diusion process within D-MRI voxels is governed by the Gaussian distribution. The covariance of the Gaussian distribution is attributed by a second-order tensor model which is a symmetric positive denite (SPD) matrix in form of P r 0(r; ) = N (r; 0; ); (2.3) = 2D; (2.4) D = 0 B B B B B B B @ D xx D xy D xz D xy D yy D yz D xz D yz D zz 1 C C C C C C C A : (2.5) The simplicity of the tensor model is of signicant advantage in the study of self-diusion in single uniform and straight tracts, which emerged to be the Diusion Tensor Imaging (DTI) technology [7, 10, 6]. It cab be shown that the voxel signal in DTI is governed by the following realtion: S/e b T Db ; (2.6) where b is proportional to the applied diusion gradient q. It can be deduced from Eq. (2.6) that at least seven independent measurements are requited for the parameter estimation of the model. The downside of the DTI method is that, when two or more tracts pass each other in a voxel, this model is no longer valid due to the partial volume eect. 18 2.3 Physical Modeling of Self-diusion Process The governing equations of the self-diusion propagator in the 3D space involve partial derivatives of spatial and temporal variables, which leads to the Partial Dierential Equa- tions (PDE). Finding the solution to the self-diusion prole in a diusive environment is not an easy task. In this section, an overview of the numerical methods for the solu- tion of PDEs is provided. Before discussing the numerical techniques employed for our study, it should be mentioned that there is a rich literature for the analytical solution to PDEs based on many sophisticated techniques. Although, these techniques provide deep theoretical insights into PDEs, they are mostly intractable in solving complicated PDEs encountered in practice. Thus, numerical techniques have to be used to nd approximate solutions for the cases which fall outside the extent of the applicability of the analytical approach. We use P r 0(r;t) to denote the self-diusion propagator in a 3D spatial environment. The conditional probability ux at the position r is given by J =DrP r 0(r;t); (2.7) whereD is the self-diusion tensor at pointr. Its physical interpretation of this quantity is the local rate of self-diusion directional exchanges of the environment. Such a tensor is dened on the point-wise basis and is in general dierent from the average self-diusion tensor of an MRI voxel. On the other hand, the self-diusion propagator is essentially a probability density function and its integral values at all time instances must be equal to 19 unity. It can be shown that this property is accounted by the continuity theorem [16], which leads to r:J =@P r 0(r;t)=@t: (2.8) Plugging the value of the probability ux from Eq. (2.7) into Eq. (2.8), we obtain the fundamental equation for the self-diusion propagator as @P r 0(r;t)=@t =r:D(r)rP r 0(r;t); (2.9) which is also known as Fick's second law. Eq. (2.9) is a linear PDE that involves the rst order temporal and the second order spatial derivatives, which belongs to the family of parabolic PDEs. For a single compartment problem including the entire space R 3 , in order to have uniquely dened solutions, it is necessary to accompany Eq. (2.9) with additional initial and boundary conditions given by P r 0(r; 0) =(rr 0 ); (2.10) P r 0(1;t)! 0: (2.11) The shape of the propagator prole and its temporal regime are strongly characterized by the self-diusion tensor D. If the environment exhibits uniform diusivity in all ori- entations without any directional preference, the self-diusion process becomes isotropic and the 3D locations of all points with the same propagation probabilities form a sphere. The isotropic self-diusion process is one of the rare cases in which a closed-form solu- tion for the self-diusion PDE exists. If the isotropic self-diusion tensor is provided 20 by D = 2D is I where I is a 3-by-3 identity matrix and D is is an isotropic self-diusion coecient of the environment, it can be shown that the isotropic self-diusion propagator is given by [16] P r 0(r;t) = (4D is t) 3=2 exp((r 0 r) 2 =4D is t): (2.12) In a multi-compartment environment such as the biological tissues, water molecules do not experience uniform resistance and they travel along dierent orientations with dif- ferent speeds. This is largely due to local variability in the material and composition of tissue microstructure. As a result, the self-diusion process is anisotropic. That is, water molecules have a faster average traveling speed in directions with less opposition. Unfor- tunately, due to the structural complexities of biological tissues, it is often intractable to nd analytical solutions to the PDE in Eq. (2.9) for such cases. Inevitably, numerical methods are used to nd approximate solutions. 2.4 Variational Formulation of Self-diusion PDE For time dependent PDEs such as the parabolic self-diusion process, the space of spatial variables must be specied at all time instances. This could be a complicated problem for the cases where the space is a function of time. Fortunately, for the cases where the geometry of the environment is not a variable of time as in the case of WTs, the spatial domain is representable by compact forms. Assume that, at each time instance, the space 21 of spatial variable r is denoted by . For x geometries, the global space of spatial and temporal variables for time interval (0;T ) is denoted by T = (0;T ): (2.13) The solution of the exact parabolic PDE given by Eq. (2.9) is attributed by strong regularity constraints. Let the regions of the interior and the boundary of T be denoted by T andC 2 ( T ), respectively,C( T ) be the space of second continuously dierentiable and continuous functions. Then, the space of the solution can be written as P r 0(r;t)2C 2 ( T )\C( T ); (2.14) which must satisfy the initial and the boundary conditions as well. These requirements are often too strong in practical applications. Thus, they are often reduced to weaker function spaces such as the Sobolev space. Before attending to denition of the Sobolev space, we introduce several related con- cepts and denitions. Denition 2.1. IfV is a normed space with normk:k V ,g2V , and 0<r2R, the open ball with center g and radius r is the set B(g;r) =fv2V ;kvgk V <rg. Denition 2.2. Suppose V is a normed space and S V . S is an open set if, for any g2S, there exists an r> 0 such that B(g;r)S. 22 Denition 2.3. Let R d be an open set, f2L p ( ), and D be a partial derivative operator of degree . The function D w f is said to be the weak th order derivative of f if Z D ! f'dx = (1) Z fD 'dx 8'2C 1 ( ): Denition 2.4. Let R d be an open set, k be an integer greater than one, and p2 [1;1]. Then, the Sobolev space is dened as the space of functions with the following characteristics. W k;p ( ) =ff2L p ( );D w f exists and belongs to L p ( ) for all ;jj kg. Also the norm on the Sobolev space is dened as kfk k;p = 0 @ Z X jjk jD w fj p dx 1 A 1 p = 0 @ X jjk kD w fk p p 1 A 1 p : As mentioned before, the regularity constraint in Eq. (2.14) is too strong to satisfy in practical numerical approaches. This is mostly due to the involvement of the second- order derivatives of functions. Thus, if this constraint is relaxed to a larger space which contains the space dened in Eq. (2.14), it will be possible to reformulate the original problem in Eq. (2.9) in a new integral (variational) form as [57] Z @P r 0(r;t) @t u(r)dr = Z u(r)rD(r)rP r 0(r;t)dr; u(r)2V; (2.15) where V is the space of the solution, which is the Sobolev space. The practical choices for k and p in denition 2.4 is k = 1 and p = 2, and the resultant space is denoted by H 1 = W 1;2 . The function u(r), which is also called the test function, could be any function extracted from space V . 23 The variational formulation in Eq. (2.15) involves the second-order derivatives. Since the computation of high order derivatives in numerical methods is complicated and also susceptible to errors, it is often desired to reduce the derivative order as much as possible. To achieve this goal, the following Divergence theorem and its corollaries [57] can be applied. Divergence Theorem 2.1. For a continuously dierentiable vector eld, F (r), over volume and the smooth and closed surface @ , the following integral equality holds Z rF (r)dr = Z @ F (r) ! ds; (2.16) where ! ds denotes the local normal to surface @ . There are a number of corollaries which can be directly deduced from the divergence theorem. One of such results directly concerns the variational formulation in Eq. (2.15). Corollary of the Divergence Theorem 2.2. Suppose that (r) and (r) are two scalar-valued functions dened on region with the closed surface@ . Also, it is assumed that the rst derivatives of(r) and the rst and second derivatives of(r) exist and are continuous over the region . Then, as a direct consequence of the divergence theorem, we have Z (r)r 2 (r)dr = Z r(r)r(r)dr + Z @ (r)r(r) ! ds: (2.17) The order of spatial derivatives in the variational formulation of the self-diusion PDE in Eq. (2.15) can be reduced using Eq. (2.17) and accounting for boundary condition in Eq. (2.11). It can be shown that, for region over which the self-diusion process 24 and its coecients are dierentiable, the weak formulation can be reduced to a simplied version as Z @P r 0(r;t) @t u(r)dr = Z D(r)ru(r)rP r 0(r;t)dr u(r)2V: (2.18) 2.5 Numerical Solution of Self-diusion PDE 2.5.1 The Galerkin Method and Method of Lines The Sobolev space in which the solution of Eq. (2.18) is sought has an innite dimension, which makes it inappropriate for numerical methods. For numerical solution, we have to limit the space of solutions to a nite-dimensional space. The Galerkin method [57, 58] provides a systematic approach that uses a sequence of embedded nite dimensional subspaces,fV n g 1 n=1 , V n V n+1 , to ll the space V as n goes to innity. Generally, Eq. (2.18) can be solved uniquely for each nite dimensional subspace. Furthermore, under a suitable assumption, the sequences of the solutionfP r 0(r;t) n g 1 n=1 whereP r 0(r;t) n V n will converge to the exact solution of equation 2.18. The self-diusion equation solution is a function of spatial and temporal variables. For numerical approximation, one common practice is to perform the discretization on the spatial domain and consider the innite domain of the temporal variable, which is known as the Method of Lines (MOL) semi-discretization of the solution. The approximate solution for each subspace V n is then provided by P r 0(r;t) n X i=1 i (t)' i (r); (2.19) 25 where ' i (r) is the basis function of the nite dimension space V n and i (t) is the time dependent coecient. By substituting Eq. (2.19) in Eq. (2.18) and examining all basis functions foru(r), the original PDE is reduced to a simple system of ordinary dierential equations (ODEs) in the form of M (t) =S (t); (2.20) where M and S are called the mass and the stiness matrices, respectively, which are direct consequences of the integral equation in Eq. (2.18), and (t) is a column vector whose entries are coecients i (t) in Eq. (2.19). 2.5.2 Numerical Solution of ODEs As discussed in the previous subsection, the simplied variational formulation of the original PDE in Eq. (2.18) can be reduced to a simple rst-order ODE, if the PDE solution is approximated by the MOL semi-discretization. This is a linear and rst order multi-dimensional ODE with constant matrix coecients. It can be uniquely solved if it is accompanied by the initial value of the vector . Hence, the solution of the original PDE in Eq. (2.15) is equivalent to the solution of following system of equations: M (t) 0 = S (t) (2.21) (0) = 0 ; (2.22) 26 where 0 is the vector consisting of the coecients of the approximate representation of the initial value function in Eq. (2.10) by basis functions of the nite dimensional space V n as provided by P r 0(r; 0) = n X i=0 i (0)' i (r): (2.23) Note that Eq. (2.21) can be rewritten as (t) 0 =M 1 S(t) = ((t);t): (2.24) For the estimation of basis function coecients three major approaches have been pro- posed [58]: 1) the best interpolant method, 2) the projection-based interpolant method, and 3) the Lagrange nodal interpolant method. The rst approach is employed in this study. It provides the most accurate yet computationally complex results so that it is appropriate for spatial domains of small to moderate sizes. The best interpolant is characterized by minimizing the following norm kP r 0(r; 0) n X i=0 i (0)' i (r)k V ; (2.25) which is achieved through global orthogonal projection of the initial value functionP r 0(r; 0) into the spaceV n . To perform such a projection, it is required to solve a linear system of n equations. Being similar to the PDE case, it is analytically intractable to nd the exact solution for ODEs. Thus, numerical techniques are employed to nd approximate solutions. For 27 ODEs arising from the semi-discretization of time-evolving PDEs by MOL, the family of one-step ODE solvers is of special interest. In the one-step method, the solution of Eq. (2.20) at time step t +4t, denoted by t+4t , can be approximated by t+4t ="( t ;t;4t): (2.26) where t is the approximate solution at time t and4t is the time step. Explicit and Implicit Euler Methods The Euler methods are among the oldest and simplest numerical ODE solvers. Their most popular feature is the simple implementation burden. However, especially in the case of the explicit Euler method, its numerical behavior may become unstable. Also, the solution precision is relatively poor as compared with higher-order methods. The explicit Euler formulation is achieved by approximation of the rst order derivatives by the forward time dierence, (t k ) = k+1 k 4t k ; (2.27) where(t k +4t k ) = k+1 . With the combination of Eqs. (2.27) and (2.24), the one-step formula for the explicit Euler method is given by " ((t);t;4t) =(t) +4t ((t);t); (2.28) or, in the discrete form, k+1 = k +4t k k ;t k : (2.29) 28 If the rst derivative is approximated by the backward time dierence; namely, (t k+1 ) = k+1 k 4t k ; (2.30) we obtain the implicit Euler method as k+1 = k +4t k k+1 ;t k +4t k : (2.31) The implicit Euler method is numerically more stable than the explicit one. However, it involves the solution of a system of linear equations. In practice, we usually have to face the stiness of ODEs, a characteristic that makes the Euler method less attractive. Stiness is a complicated phenomenon that involves the characteristics of the under- lying ODE, the numerical solution, and the time step. Its exact denition and properties are beyond the scope of our discussion. The simple interpretation of stiness is the co- existence of several terms with dierent temporal rates in the solution, which makes the choice of time step more dicult. If the time step is too small, it makes the convergence very slow, which will result in accumulated errors due to more time steps. On the other hand, if the time step is too large, some abrupt variations of the solution might not be captured by the numerical method, which consequently causes deviations or even diver- gence in the solution. The best known one-step scheme to address this issue is the family of higher-order implicit Runge-Kutta (IRK) methods. 29 Higher-order Runge-Kutta Methods For practical situations, there are cases in which stiness does not exist. Thus, it makes the employment of explicit methods possi- ble. The advantage of this approach is that no linear system of equations is to be solved. One of the commonly utilized explicit methods is the higher-order explicit Runge-Kutta (ERK) [14] schemes. The benet of the ERK method instead of lower-order explicit Euler method is the high degree of solution accuracy achieved at the cost of more algorithmic computations. The s-stage ERK algorithm is provided in the following, z i = 0 @ k +4t k i1 X j=1 a ij z j ; t k +c i 4t k 1 A ; c 1 = 0; (2.32) k+1 = k +4t k s X i=1 b i z i ; (2.33) where the parameters c i , a ij , and b i are given at the stage s. Formulation of implicit Runge-Kutta (IRK) is similar to ERK, but it is required to solve a system of linear equations for the solution z i = 0 @ k +4t k s X j=1 a ij z j ; t k +c i 4t k 1 A ; c 1 = 0; (2.34) k+1 = k +4t k s X i=1 b i z i : (2.35) The major dierence between ERK and IRK is that z i is provided implicitly in IRK. Although the solution ofz i involves higher numerical burden, the fact that the solution is achieved through contributions of allz i values makes the method robust against stiness. 30 2.5.3 Finite Element Methods It was shown in the last subsection that, to implement IRK methods, we are engaged in the solution of a system of linear equations. Considering the large spatial dimen- sions of practical problems and the number of basis functions involved, those systems of equations are usually large, which could easily drain the memory and computational resources of the system. Besides, the performance of numerical solvers highly depend on the characteristics of the stiness and mass matrices, introduced in Eq. (2.20). To improve the performances of the linear solvers, it is desired that those matrices be sparse with small condition numbers [57]. These matrix properties are strongly attributed by the integral terms in Eq. (2.18), which are highly in uenced by the characteristics of the basis functions. For approximate solutions with a large number of basis functions, it is dicult to fulll these requirements with basis functions of unlimited spatial support. As a solution, a partitioned space with basis functions of limit spatial support can be employed. The most famous method towards this objective is the Finite Elements Methods (FEM). The FEM is among the superior numerical methods for the solution of PDEs in terms of solution accuracy and its ability to handle irregular boundary conditions. The methodology is based on dividing the spatial domain to non-overlapping small volumes called elements. The space of functions within each element is approximated by the linear summation of basis functions. The basis functions are designed such that their spatial supports are limitted to the element dimension. 31 One common choice for the basis functions is the piece-wise polynomial functions such as the Lagrange and Lobatto shape functions. The local renements of approximate solutions are controlled by the size of elements as well as the degree of the polynomial basis, which are termed as h and p renements, respectively. The FEM with high order polynomial basis functions is referred to as the higher-order FEM. It is known that, for higher-order FEM, the family of the Lobatto shape functions are preferable in terms of implementation burden and their consequence on the behavior of stiness and mass matrices. Having the spatial domain partitioned into nite element regions, it is appropriate to reformulate the simplied variational formulation in Eq. (2.18) in terms of nite volume integrals as M X i=1 Z K i @P r 0(r;t) @t ' j (r)dr = M X i=1 Z K i D(r)r' j (r)rP r 0(r;t)dr; ' j 2V n ; j = 1; 2;;n; (2.36) where M is the total number of elements and K i is the region associated with element i. In the 3D space, there are a number of elements proposed for partitioning of the global space [58]. The most common ones are brick, tetrahedral and prismatic elements as illustrated in Fig.2.8. The basis functions of an element is associated with its vertices, edges, faces and interior regions, and named accordingly. For example, the basis functions associated with vertices are referred to as vertex basis functions. Due to implementational simplicity, early FEM methods employed vertex basis functions of lower order polynomials. There is no mechanism to increase the degree of polynomial basis functions in lower order FEM. 32 Figure 2.8: Three common 3D elements: the brick (left), the tetrahedral (middle) and the prismatic elements (right). Thus, the only way to increase the spatial resolution of the solution is to decrease the size of elements. This is an inecient approach and often constrains the scalability of the method to large and inhomogeneous environments such as the biological tissues. In recent years, the extension towards the higher-order FEMs by denition of edge, face, and interior region polynomials of larger degrees has received a lot of attention. Having the capability to adjust the order of polynomials for each individual element alongside ecient partitioning techniques (i.e. mesh generation) provides ecient allocation of resources and makes more complicated problem solvable. Dening basis functions and the interior region for every individual element is not ecient and even might be impractical for large systems. To this end, it is customary to dene those identities for a reference element, ^ K, and their counterparts in any arbitrary element K are achieved through a mapping system X K : ^ K!K. By incorporating this mapping relation into the integral terms in Eq. (2.36), the integral region can be trans- formed from the region K to the region ^ K. The general method for construction of such a reference mapping is through a procedure called transnite interpolations. However, 33 it is usually unnecessary to go that far in practice and a simple isoparametric mapping suces for many problems. An example of such a mapping will be provided in the next chapter for the tetrahedral elements. 34 Chapter 3 : Anisotropic Self-diusion Simulation 3.1 Introduction Recent advancements in medical imaging and computer science have signicantly revo- lutionized the eld of Diusion Magnetic Resonance Imaging (D-MRI) and tractography in the past two decades. It is the only noninvasive technique known for unraveling the microstructures of the white tract tissues in the central nervous system (CNS). TheD-MRI signal is characterized by the measurements of the water molecules Brow- nian motions within the biological tissues which is a physical process referred to as self- diusion. The random motion of an individual molecule is too complicated to study, but it is possible to predict the average behavior of all molecules during relatively large time period. This time period which is referred to in the literature as the long term diusion time [16] is characterized by the time interval after which the motion of a water molecule is uncorrelated with the other molecule motions. The applicability of D-MRI method is rooted in an important assumption that the directions of the water molecule motions are identical to the orientation of the tissues. Well-organized tissues exhibit anisotropic 35 self-diusion, the property which makes them potential candidates to study by D-MRI. Such tissues are observed in the brain white matters, spinal cord, and muscles. D-MRI principles were developed in mid 80s, by incorporating the early work of Stejskal and Tanner [60], for encoding molecular diusion eects, into the Nuclear Mag- netic Resonance (NMR) imaging principles [37, 46, 66]. The D-MRI signal acquired from a biological tissue is a measure of the underlying self-diusion processes within the microstructure environment of the sample. Therefore, the general problem in the applica- tion of D-MRI methodology to the study of the microgeometries of material is retrieving the self-diusion prole from the observed D-MRI signal. This is an ill-posed problem and often additional constrains are required for identifying unique solutions. For the successful assessments of underlying tissue microstructure, it is necessary to have thor- ough understanding of the relations between tissue structural parameters such as their geometry packs, morphologies, and constituting material, and the self-diusion regime. Quantication of the self-diusion propagator within the brain microstructure is not a trivial task, and a common trend is to apply parametric or non-parametric models, such as the tensor model. Despite the simplicity of such strategies, their extents of validations are up to satisfaction of some criteria. Considering the domain of variability of the brain white matter microstructure, those assumptions could easily be violated and the models will no longer be valid. In early 90s, the tensor model for characterization of self-diusion anisotropy in the brain white matters revolutionized the D-MRI methodology which emerged to Diusion Tensor Imaging (DTI), a method which continued to retain its signicance to the present date [7, 10, 6]. The directional preferences of self-diusion process which is due to the 36 structure of the tissues of study, are characterized by the tensor eigenvectors and inten- sity of self-diusion in each direction is represented by the corresponding eigenvalues. Naturally, such a simple model is suitable for capturing the behavior of the Gaussian self- diusion propagation, which makes them prone to inaccuracy for the non-Gaussian cases such as crossing-bers. The mixture of tensors is another parametric model proposed to address the problem of partial volume eect for more intricate cases of crossing white tract in the MRI voxel [15]. The major shortcomings of such approaches are their inadequate exibilities to accommodate variability of self-diusion proles, as well as lack of direct relations between the white tract microstructure parameters and their characteristics. As an alternative view to this problem, a dierent family of techniques is proposed, trying to account for the microstructure parameters in quantication of the self-diusion prole, which fall under two main subcategories. The rst subgroup involves the methods which attempt to analytically uncover the exact solution of the self-diusion prole in the brain tissues. Given the complicated geometries of the white tracts, the scope of these techniques is limited to over-simplied models of axons and expanding their applicability to more general geometries will extremely be dicult, if not impossible. However, the results of such attempt could be benecial for providing ground truths for validation of numerical methods. The reader can learn about dierent visions of this class of methods by referring to [56, 11, 59, 61, 51]. To overcome the intractability of analytical methods for the general geometries, nu- merical methods and simulations have been utilized. Monte Carlo (MC) simulations of water molecule self-diusion in simple models of axon microstructure were introduced by Szafer et al. [62]; Ford and Hackney [27]; and Ford et al. [28], for the purpose of D-MRI 37 signal simulations. In these approaches, the random walk motions for a limited number of particles within the models of axons are simulated, as the time evolves. One major con- cern associated with MC simulations is that the actual eect of the self-diusion prole in D-MRI signal reconstruction might be underrepresented, since a relatively limited num- ber of particles for simulation are considered, which are far less than the actual amount of water molecules contribute in signal formation. Secondly, it is not straightforwardly clear how to scale the results of computer simulations which are usually performed for a few axons, to average self-diusion propagation of the voxel which involves thousands of bers. Hwang et al. [34] exploit the Finite Dierence (FD) method to tackle the solution of self-diusion approximate Partial Dierential Equation (PDE) in simple cylindrical models for axons. They include the eect of thick myelin sheaths into the calculation in terms of eective permeability of the axon membrane layers. Similar to MC methods, they provide no answer to scaling the derived microscopic results for bundle of axons to computation of average self-diusion proles of voxels, which are in the macroscopic range. Besides, their methods do not provide mechanisms to accommodate crossing con- ditions for the white tracts. And, nally, it is not an easy task to accommodate the complex tracts geometries with the FD methods. In our view, the solution of the PDE of self-diusion propagation, in its most general form has not properly been addressed in the literature. Bearing in mind the signicance of this problem to dierent aspects of D-MRI technology, necessitates the development of numerical tools for addressing this issue. To this end, the purpose of this study is to target the solution of the PDE by relaxing most of the constraints imposed by the previous work and include a vast range of variability in geometries of the white tracts 38 and composition of the constituting material. To achieve such ambitious objective, the higher-order Finite Element methodology (FEM) is employed [58, 21, 29, 19]. The FEMs are believed to be among the high performance numerical methods for solution of PDEs with complicated boundary conditions. In spite of their superior exibility and precisions, sucient care must be paid in the design of FEMs for the cases of multi-compartment geometry packs, such as in the white tract scenario. In the method section, the FEM formulation of the problem is presented alongside our especial design to make the methodology tractable for the problem of self-diusion in multi-compartment environments. Furthermore, we provide a numerical assessment scheme to validate the accuracy of the FEM results. The D-MRI signal is sensitive to average self-diusion proles of the voxels, rather than their point-wise values. As a distinguished highlight of this work, the average self-diusion propagations of voxels are formulated. This is an essential step forward, towards understanding the relation between microscopic parameters of white tracts which are generally unobservable in practice, and the average voxel quantities which are in the macroscopic range and are measurable withD-MRI. In the result section, the eects of dierent white tract components on the trend of propagation of water molecules are illustrated, as well as the average prole of a voxel under unusual tract conditions such as ber crossings and demyelinations. In the discussion section, the signicance of the results and the corresponding likely applications to dierent aspects of D-MRI is interpreted. Finally, we conclude our discussion by providing the potential directions for the future explorations. 39 3.2 Theory and Methods It is theoretically well-understood that the observed signal of a D-MRI voxel is strongly tied with the self-diusion prole of the spin magnets (e.g. water protons) within the voxel. Such relation constitutes the foundations of many applications of this imaging technique for exploring the structures of white matters in the CNS. The D-MRI signal, as a result of the self-diusion process, is governed by the following principal equalities [16] S(k;q) / Z R 0 (r 0 )e 2ik:r 0 Z R P r 0(r; )e 2iq:(r 0 r) drdr 0 (3.1) E (r 0 ;q) / F [S(k;q)]; (3.2) where (r 0 ) is the concentration of the spin magnets which is a function of the voxel tissue material, k is a sample extracted from the k-space, q is proportional to the ap- plied diusion gradient vector, is the eective diusion time during the course of signal acquisition, and S is the k-space signal. The value of k is a function of the applied magnetic eld gradients of conventional - non-diusion- MRI [31] along the x-, y-, and possibly z- axes( usually D-MRI signal is acquired in 2D, so the k-space sampling is per- formed in x- and y- imaging axes). Also,E (r 0 ;q) represents the D-MRI signal observed at the position r 0 , as s result of applying diusion gradient q and eective diusion time . The intensity of this signal is inversely proportional to intensity of the self-diusion along the diusion gradient at position r 0 . The generally unknown physical quantity in equation 3.1, which plays a substantial role in unraveling the structures of white matters, is the self-diusion propagator,P r 0(r; ).It is dened as the probability density function 40 of transition of a water molecule from the coordinate position r 0 at time t = 0, to the position r, at time t = . The shape and evolution of propagator is severely impacted my the underlying tissue structure, especially in the case of the white tracts which the organization and geometry pack of axonal bers play predominant role in shaping the 3D prole of the propagator. Under controlled imaging conditions inD-MRI scanners, all the participating quanti- ties in Eqs. 3.1 and 3.2, except for the propagator, are known- the eect of inhomogeneous (r 0 ) is alleviated as the result of normalization by the signal acquired with zero diusion gradient, i.e q = 0. A major endeavor of this work will be devoted to quantication of the propagator as a function of spatial and temporal variables. The self-diusion propagator is governed by a parabolic PDE [57] as given below (For the purpose of notation brevity, the notation P (r;t) to represent the cases where the origin is transformed to r 0 , @P (r;t) @t =rD(r)rP (r;t); (3.3) whereD(r) is the self-diusion coecient at the position r. The parabolic term species that the PDE is a function of spatial and temporal partial derivatives and the solution to such a problem is a time evolving function. The solution of the equation is uniquely specied, if it is accompanied by additional initial and boundary conditions. For single 41 compartment environments (i.e. diusive media without inner layers), they are provided as simple initial and boundary conditions as given below, P (r; 0) = (r); (3.4) P (1;t) = 0: (3.5) Here, a compartment is referred to a spatial 3D region in which all the variables in- volved are continues and all the participating derivatives exist. The D-MRI principal signal generation equations concern the diusion prole within the brain white matter microstructure. Having the geometry of the microstructure, it will, in principle, be pos- sible to compute the propagator by solving the system of equations 3.3-3.5 within the given geometry (throughout this article, the true solution of equation 3.3 is referred to as the exact solution, and the solution acquired by numerical methods are referred to as the approximate solution). The brain white matters consist of bundles of axons which exhibit dierent self- diusion speeds at dierent locations. Consequently, the volume of the tissue of interest is partitioned into a number of compartments such that the self-diusion coecient is continues within a compartment, and it has discontinuities on the surface boundaries between neighboring compartments. Hence, the solution of the original PDE in equation 3.3 has to be examined in a multi-compartments scenario. Unfortunately, the additional conditions provided in equations 3.4 and 3.5 will not be adequate to uniquely identify the solution. 42 For a white matter environment consisting of the cytoplasm, the myelin the mem- brane, and the extracellular compartments, two additional boundary conditions must be satised. For transition from the myelin to the extracellular regions, the only condition to satisfy by the solution is the continuity of the normal ow[13, 70] which is dened as, J (r) = D (r)rP (r):~ n(r) (3.6) J + (r) = J (r); (3.7) whereJ is referred to as the self-diusion probability ow, and~ n is the normal vector to the boundary surface. This equality is of special signicance for coupling the quantities of propagators of two adjacent compartments on the points located on their shared bound- ary, though its validity within a compartment is a trivial fact. The + and subscripts are to discriminate between the quantities on two opposite sides of a layer. As a result of partial permeability of the thin membrane wall, for boundaries between the cytoplasm and myelin regions, the boundary condition becomes more complicated. Tanner [65] argued that the eect of a thin partially permeable boundary can be captured by enforcing the following additional boundary condition [65, 70], P + (r;t)P (r;t) = 1 J (r); (3.8) where is the membrane wall permeability constant. As a necessary condition for the solution of Eq. 3.3, the environment in which self- diusion occurs must be specied rst. That is, it is required to have the geometry of 43 microstructure, and the values of self-diusion coecients for all locations. In general, such geometries are very complicated and nding exact solutions will become analytically intractable. Inevitably, it is necessary to employ numerical methods to approximately dis- cover the answers. Often, there is a compromise between the complexity of the numerical method utilized, and the precision of the approximate solution. It is a common belief that the FEMs are among of the superior methods, in terms of the representation of problems and precision of the results, for tackling the solutions of PDEs with irregular boundaries/layers conditions [57]. Their special characteristics provide the exibility to accommodate a variety of problem constraints with minimal eorts. Subsequently, con- sidering the structural complications of the brain white matters, the FEMs appear to be potential numerical tools. For the purpose of modeling the geometries and structures of the white tracts, the hexagonal array of cylinders will be adopted for experiments of this study. The advantage of using this model is that the volume fraction of axons in this model is very similar to the one from the actual organization of axons in real white tracts. It is worth mentioning that the extent of FEM's applicability goes beyond this simple model, and it is capable of handling far more complex geometries including the actual white tract structures. Without loss of generality, the cylindrical model is adopted to dene the geometries of self-diusion environments. In the model of hexagonal arrays of cylinders [56], a myelinated axon is represented by an inner cylinder as the cell cytoplasm region, and a coaxial cylindrical sheath that models the thick myelin layer around the axon cell. The eect of the axon membrane partial permeability is lumped into the permeability of the myelin sheath. Also, the organization of myelinated axons within a white tract is assumed 44 to be hexagonal, as illustrated in Fig.3.1. In this model, the myelin sheath thickness is considered and its magnitude order is comparable with the axon diameter. Also, it is further assumed that the self-diusion coecient is constant within a compartment (There are three dierent compartments identied in the provided model which are cytoplasm, myelin sheath, and extracellular regions). Figure 3.1: Axons model: (Left) A cylindrical model for a myelinated axon, where self- diusion coecients are denoted by D m , and D c for the inner cell and myelin regions, respectively, and (Right) the hexagonal array organization of axonal bers 3.2.1 Variational Formulation and FEM For the time being, we assume that the membrane layers are not present and the white matter tissue simply consists of the organization of the cytoplasm, the myelin, and the extracellular materials. However, we will show how to incorporate the eect of the mem- brane walls into the equations in Section 3.2.3. Under this assumption, the only boundary condition to enforce is the continuity of the propagator ow in Eq.3.7. The FEMs are among advanced methods for the solution of PDEs in their most general forms. In spite of their powerful potentials, they could not be directly applied to 45 the problem of self-diusion in the brain for a number of reasons, which will be claried below. Fortunately, the FEMs come with a host of nice properties that control the characteristics of the solutions. Those properties can be exploited in formulation of the problem such that the resultant numerical solutions become tractable. Since our intention is to nd the solution of equation 3.3 by means of numerical meth- ods, the partial derivatives involved have to be computed numerically. Such numerical approximations are usually susceptible to errors especially in the case of second order derivatives which are very sensitive to the system noise level. In order to reduce the order of second order derivatives, the common practice is to transform the original PDE in 3.3 into the variational domain. As the rst step towards the numerical solution, the space of the original PDE is transformed to the variational domain by multiplication of the both sides of Eq. 3.3 by an arbitrary function u(r) and integrating over the whole single compartment volume . Selection of the function space from which u(r) is extracted, depends on the problem characteristics. The common choice, which covers a vast range of practical problems, is the space of square-integrable functions L 2 , Z @P (r;t) @t u(r)dr = Z u(r)rD(r)rP (r;t)dr: (3.9) The divergence theorem (see chapter 2.) can be utilized to reduce the order of spatial derivatives on the right hand side of Eq. 3.9. As mentioned in the previous chapter, for application of the divergence theorem, some conditions must be satised. Due to discon- tinuities of self-diusion coecients on compartment boundaries, the integrant function 46 is not dierentiable at those regions and the theorem is not applicable everywhere. This problem can be eliminated by partitioning the domain of integration to volumes with continues self-diusion coecients, which are the compartments associated with the cy- toplasm, myelin, and extracellular regions. Hence, Eq. 3.9 can be rewritten as X i Z C i @P (r;t) @t u(r)dr = X i Z C i u(r)rD(r)rP (r;t)dr; (3.10) and application of the divergence theorem can reduce the order of spatial derivatives involved as X i Z C i @P (r;t) @t u(r)dr = X i Z C i D(r)ru(r):rP (r;t)dr + Z S i u(r)D(r)rP (r;t): ~ ds ; (3.11) whereC i species the region associated with compartment i, andS i is the outer surface ofC i . For the exact solution, the application of layer boundary conditions in Eq. 3.7 eliminates the second integral (the surface integral) on the right hand side of equation 3.11. However, for approximate solutions, this property is not automatically held, and the design of the solution must be adapted to contain the property. To understand how this constraint is imposed on the approximate solution, it is important to understand the properties of the FEMs and design the approximate solution such that it accounts for the layer boundary conditions. In the FEMs, the spatial domain is partitioned into non-overlapping small 3D volumes which are referred to as elements. The hexagonal array of cylinders used to model axonal bers and the corresponding FEM partitioning by tetrahedral elements (known as the 47 tetrahedralization mesh) are illustrated in Fig. 3.2. The solution within an element is approximated by the linear participation of nite number of the element basis functions. For further information about the tetrahedral FEM and piecewise polynomial basis func- tions, which are utilized for the implementation of this study, please refer to section 3.3 and discussions in [58, 33, 23]. Figure 3.2: Tetrahedralization: (Left) Hexagonal array of cylinders model for ber tract, (Right) 3D mesh as the result of tetrahedralization To have the approximate solution fulll the boundary conditions, special considera- tions must be accounted in the design of the tetrahedralization and basis functions. We elaborate our approach to achieve the desired characteristics and simplications of the variational formulation in Eq. 3.11 in section 3.5. It is shown in section 3.5 that If the tetrahedralization is performed such that, no boundary passes through the interior regions of elements and the Minimum Rule of conformity to the Sobolev space,(see section 3.4) is violated for all elements, the approximate solution will have the exibility to accom- modate layer boundary conditions. It is also shown that, the variational formulation is 48 transformed to a simple rst order ordinary dierential equation (ODE) of the following form, M (t) = S(t) (3.12) (t) = [ 1 (t); 2 (t);); n (t)]; (3.13) where vector (t) consists of the coecients of the basis functions in the approximate solution (see section 3.4) at time instance t, and S and M are square stiness and mass matrices, whose dimensions are equivalent to the total number of basis functions from all elements. There are quite a few numerical methods for the solution of the system of linear ODEs as the one in Eq. 3.12. Considering the fact that the self-diusion environment is multi-compartment, the temporal variation of the solution at dierent locations could exhibit dierent rates especially at the time instances when the diusion wave arrives at a location. As a result, the ODE in Eq. 3.12 could become sti, which necessitates the employment of implicit numerical methods, such as the family of implicit Rung-Kutta (IRK) methodologies [14]. Application of simple explicit methods such as the Euler explicit solution requires the selection of very small time steps which makes the solution to converge slowly and also increases the accumulated errors. The value of vector at the time t = 0, is exploited for the initialization of the ODE solution. This vector can be calculated by the projection of the initial condition function provided in Eq. 3.4 into the space of basis functions. This is accomplished by 49 minimization of a norm distance between the actual initial function and the approximate representation. P (r;t) X i i (t = 0)' i (r) H ; (3.14) where ' i (r)s are the basis functions of the FEM, and k:k H is the appropriate norm dened on function space H. The choice of the function space is made according to the nature of the problem and the degree of computational complexities it imposes(e.g. the Sobolev space). Taking into the account the fact that the initial value function is in form of the Dirac delta function, its approximation by a linear combination of piecewise polynomials is inaccurate, irrespective of the number of contributing basis functions. On the other hand, the self-diusion propagator is essentially a probability density function whose integral value over the global spatial domain must be united by one at all time instances. This conservation rule for the propagator is implied in the original PDE in Eq. 3.3 and there is no extra explicit mechanism to enforce this property. As a result, any error induced in the calculation of the initial value vector, will be preserved and propagated as the time evolves. To this end, it is crucial to suppress the error of the initial value function approximation, by minimizing the estimated error of Eq. 3.14. To do this there are two options. The rst choice is to use dierent (non-polynomial) basis functions which are capable of modeling the dirac function. This could lead to alteration of the well-established theories of FEM which could engage very dicult situations. The more practical solution is to use initial value functions which are better representable by the polynomial basis function. 50 The following observation is of special importance in handling this problem. Let us narrow our attention on the quantity of the propagator of a point within a compartment, during a very small time interval [0;t " ], wheret " is a small time interval right after initial timet = 0. If the value oft " is suciently small such that the water molecules do not nd enough time to arrive at any restricting layer boundaries, the process of self-diusion will be free (unrestricted) for the mentioned time interval (see section 3.6 for formal proof). A closed form solution for unrestricted self-diusion process exists [16] which can be exploited for evaluation of the global behavior of the propagator during the mentioned time interval, P free (r;t) = (4Dt) 3 2 e [r 2 =(4Dt)] ; t2 [0;t " ]: (3.15) It can easily be examined that by letting t!1, Eq. 3.15 converges to the Dirac delta function. As a result, if the system is initialized at time t = t " with the initial value function P free (r;t " ), the solution must be unchanged. This statement is supported by the fact that the ODE in Eq. 3.12 is linear with time invariant matrix coecients. The advantage of this time shift is that the initial value function is now an exponential, whose approximation with linear combination of piecewise polynomials is far more accurate that the original delta function. For assessing reasonable values for t " the Einstein equation [16] for calculating the self-diusion length can be employed. That is, l = p 2DT " (3.16) 51 where l is the distance of the point of interest within the compartment, to the nearest boundary. The Einstein equation is derived by assuming unrestricted self-diusion and extracting the variance of the Gaussian propagator. Eq. 3.12 involves the contributions of the basis functions from all elements of the sample mesh. For a large sample, there is a possibility that water molecules do not have time to reach some of the mesh elements within the maximum diusion time of interest . Therefore, their propagator values will be united by zero for the whole time of study, and the coecients of their basis functions will be zero. Such zero coecients do not contribute to the solution, and merely add undesired dimensionality to the ODE in 3.12. Theoretically, this could impose no harm to the solution, but in practice, the system memory and computation resources could dramatically be drained, and the performance of the numerical linear solver for solution of the ODE, could severely be degraded. For exclusion of those ineective elements, we introduce the concept of eective radius of propagation for multi-compartment self-diusion media. The eective radius of propagation, is characterized as the upper bound for the max- imum range of water molecules self-diusion propagation in a multi-compartment envi- ronment during the time interval . If the largest self-diusion coecient exposed by the environment is denoted byD max , the maximum range of propagation at any location within the compartments is upper-bounded by the free self-diusion length with that maximum self-diusion coecient. This is, due to the fact that free-self diusion occurs 52 faster than restricted one, since there are no boundaries to hinder the propagation of wa- ter molecules. This free self-diusion length denes the eective radius of propagation, which could be calculated by the Einstein relation via R eff = p 2D max : (3.17) For implementation, it is sucient to include the elements which fall within the sphere of radius R eff , only. This could dramatically reduce the size of stiness and mass matrices in Eq. 3.12. 3.2.2 Validation of the Solution To condently extend the results of solution of the self-diusion propagator to clinical applications, it is essential to validate the obtained results. This could be quite challenging since the ground truth is not available for evaluation of the FEM's results for general geometries of microstructure. For simple geometries, it might be possible to nd closed- form solutions to validate the results, but this is hardly extendable to general complicated geometries. A valid solution is identied if it simultaneously satises the original PDE Eq. 3.3, the initial and boundary conditions provided in Eqs. 3.4 and 3.5, and the layer boundary constraint of Eq. 3.7. According to the above discussions on the time shift of the sys- tem initiation time, the initial condition function is approximated with the high level of accuracy, provided that the applied norm is carefully chosen for the problem. If the for- mulization of FEMs is conducted for the spatial domain conned to the regions specied 53 by the eective radius of propagation, the boundary condition at innity is automatically satised, too, since it is assumed that the value of the propagator beyond that radius is zero. It was previously stated that if the provisions provided in section 3.5 are accounted for, the boundary conditions are accommodated by design. Therefore, for solution val- idation, it just remains to conrm it in the original PDE. The FEM solution provides the ability to evaluate the evolution of the propagator in time and space. This makes it possible to numerically evaluate spatial and temporal partial derivatives of the propaga- tor function. The solution is validated by conrmation of the partial derivatives and the self-diusion coecient in the original PDE. For the cylindrical axonal model for which the self-diusion coecient is constant within a compartment, the validation test reduces to the following simple rule, D C i = @P (r;t) @t r 2 P (r;t) ; (3.18) whereD C i is the self-diusion coecient of the compartment of the study. Naturally, the sanitation of this method is up to accuracy of estimation of partial derivatives. To reduce the error due to numerical computations of derivatives, the Eq. 3.18 can be examined at the several samples of a compartment and the the results are deduced through averaging. 3.2.3 Membrane Wall Eect The membrane is an extremely thin layer containing the cytoplasm region. Despite its negligible volume, it is a porous medium that exhibits partial permeability. As a result, the restricting eect of the membrane walls on the self-diusion process is substantial and should be eectively included in the derivation of the propagator solution. 54 The boundary condition on an innitesimally thin permeable boundary is given in Eq.3.8 [65] which was inspired by the analogy to the thermal diusion problem [51]. In its current implementation, our method is incapable of direct accommodation of Eq.3.8. However, Powels et al. [51] provided an equivalent nite thickness boundary model for the membrane wall, which could be integrated into our method straightforwardly. Figure3.3 illustrates two equivalent models for the membrane wall. Figure.3.3.a shows an extremely Figure 3.3: Membrane equivalent model:a) An innitesimally thin porous boundary. The boundary condition is governed by Eq.3.8. b) Equivalent nite thickness layer with simple condition in Eq.3.7 on the boundaries. The two models are equivalent as long as =D b =b thin partially permeable boundary with the permeability constant , and Figure.3.3.b presents the equivalent nite thickness layer with the thickness b and self-diusion co- ecient D b . It can be proven [51] that, for a small value of b, as long as the condition =D b =b holds, the realizations of the self-diusion propagator in the neighboring cyto- plasm and myelin regions are the same. That is, the boundary condition in Eq.3.8 can be safely replaced by inserting a nite thickness compartment and the new boundary condition on the interfaces of the new added layer reduces to Eq.3.7. This is a perfect t for our developed method in Section 3.2.1. 55 3.2.4 Voxel Average Propagator The source ofD-MRI signal is the microscopic motions of water protons within the axon microstructure. The spatial scale of such phenomena is in the order of [m] which is beyond the resolution of D-MRI measurement. On the other hand, by D-MRI data acquisitions, the average behavior of all water protons within the voxel is manifested. Naturally, understanding the microscopic structure of axons from D-MRI is desirable too. To achieve such a goal, it is crucial to understand the relationship between these microscopic and macroscopic processes. In this section, we provide a derivation to relate the self-diusion within a set of axons in a hexagonal array of cylinders model to the integral average values over the MRI voxel. The application of the Fourier transform to theK-space signal provides theD-MRI signal at a point of interest, r 0 , as E(r 0 ; ;q)/(r 0 ) Z R 3 P r 0(r; )e 2iq:(rr 0 ) dr: (3.19) However, due to practical limitations such as T 2 relaxation [31], the sampling of the K- space is often imperfect and the reconstruction of the signal is prone to local averaging over voxels, rather than point-wise acquisitions as in Eq. 3.19. Hence, an estimate of the observedD-MRI signal in the voxelV j could be obtained by averaging this equation over all spin magnets of the voxel. 56 E av V j (;q) / 1 N j Z V j Z R 3 (r 0 )P r 0(r; )e 2iq:(rr 0 ) drdr 0 = Z R 3 e 2iq:R ( 1 N j Z V j (r 0 )P r 0(r 0 +R; )dr 0 ) | {z } P V j (R;) dR; (3.20) where N j = R V j (r 0 )dr 0 denotes the total number of spin magnets in the voxel, and P V j (R) is the weighted average of the propagator over the voxel V j . It is known that [56] for long term self-diusion processes (by long term it is meant the time interval after which the motions of water protons are uncorrelated), the concentration parameter (r 0 ) is uniform and constant for the entire regions of a compartment. It is, therefore, possible to evaluate the average voxel propagator by partitioning the voxel volume into its constituting compartments. P V j (R; ) = 1 N j X k ( k) Z C k P r 0(r 0 +R; )dr 0 = X k f k P k (r; ); (3.21) where k is the constant concentration parameter, f k = k V k P m mVm (V k is the volume of compartment k conned to the voxel region) is the spin magnet fraction, and P k (R; ) is the average propagator of compartment k. In the subsequent sections, the signicance of equation 3.21 in derivation of average voxel propagator for dierent white tract condition will be more clear. 57 3.3 The Finite Element Method (FEM) In the FEM, the spatial domain is partitioned into non-overlapping small volumes, which are referred to as element domains. The space of functions within an element is ap- proximated by the nite number of basis functions. The common choice for the basis functions is piecewise polynomial functions, the polynomials whose supports are limited to the dimensions of the element. There are a number of options for the choice of elements in 3D, such as bricks, tetrahedral, and prismatic nite elements. Since the isoparamet- ric mapping (dened below), has constant Jacobean in the case of tetrahedral elements, the computation of the stiness and mass matrices, as introduced in section 3.5, can be performed on the reference domain. This will dramatically reduce the computational complexity of the FEM. To this end, we will conne ourselves to the tetrahedral elements in the current study. The procedure of representing the spatial domain with tetrahedral elements is referred to as tetrahedralization. This goal is achieved by mesh generating algorithms such as Delaunay triangulation [24](refer to Fig. 3.2 for a sample tetrahe- dralization mesh). For meshes with pure tetrahedral elements, all elements and their basis functions, which are referred to as physical elements, are represented by a reference tetrahedron and its basis functions, respectively, through the isoparametric mapping. Fig. 3.4shows the reference tetrahedral element. The spatial domain of the reference tetrahedron is dened by the following relation, ^ K = n 2 R 3 ;1 1 ; 2 ; 3 ; 1 + 2 + 3 1 o (3.22) 58 Figure 3.4: Tetrahedral element: Reference tetrahedral element and the reference coor- dinate The space of functions within the reference element is approximated by a linear space, spanned by the basis functions of the reference domain. As mentioned previously, the common choice is the piecewise polynomial functions. The family of Lobatto polynomial functions [33, 58], is of special interest for establishing the basis polynomial functions, due to their well-behaved numerical characteristics. Associated with the reference tetrahedron vertices, edges, faces, and interior regions, four types of basis functions can be dened. Having the following ane coordinates, 1 = 2 + 1 2 ; 2 = 1 + 2 + 3 + 1 2 ; 3 = 1 + 1 2 ; 4 = 3 + 1 2 ; (3.23) we can dene the following four categories of basis functions. 59 Vertex Functions: there is one basis function associated with each vertex, ^ ' v i = j ; i = 1; 2;:::; 4 (3.24) where j is the index of the element face, opposite to vertex i. Edge Functions: For every edge, according to the maximum allowed degree of edge polynomial functions, p e j , k basis functions are dened as, ^ ' e j k = j1 j2 k2 ( j1 j2 ); 2kp e j (3.25) where j1 and j2 are the indices of faces which share only one vertex with the edge, and k (x) is the kernel polynomial of the Lobatto shape functions [33]. Face Functions: For every face of the reference element, according to maximum al- lowed degree of face polynomial functions,p s i , (p s i2)(p s i1) 2 basis functions can be dened. ^ ' s i n1;n2 = A B C n11 ( B A ) n21 ( A C ); 1n1;n2; n1+n2p s i 1 (3.26) where A;B, and C are the indices of the face vertices, and A, and C correspond to the smallest and largest indices, respectively. 60 Bubble Functions: It is possible to dene (p b 3)(p b 2)(p b 1) 6 bubble functions associ- ated with the reference element interior region, wherep b is the maximum degree of bubble polynomial functions. ^ ' b n1:n2;n3 = n11 ( 1 2 ) n21 ( 3 2 ) n31 ( 4 2 ) 4 Y i=1 i : (3.27) The projection from the reference tetrahedral element ^ K to an arbitrary element within the mesh, denoted by K, is achieved through the isoparametric mapping, presented in equation 3.28, x K () : ^ K!K : x K () = 4 X i=1 ! V i ^ ' v i (); (3.28) where ! V i is the 3D coordinate vector of vertexi of elementK. As a result of this mapping, the basis functions of element K can be represented according to basis functions of the reference element. That is, '(K) = ^ 'x 1 K (x): (3.29) 3.4 The Minimum Rule In order to nd the solution of PDEs by the FEMs, the desired solution is approximated by linear contributions of basis functions of mesh elements. Such a representation for the PDE of the self-diusion propagator has the form specied in Eq. 3.30, which is also recognized as discretization by the Method of Lines [58]. P (r;t) = n X i=1 i (t)' i (r): (3.30) 61 Often, it is necessary to constrain the basis functions to guarantee the smoothness of the solutions at the boundaries of elements. The common criteria for smoothness are the existence and square-integrability of spatial partial derivatives of the solution. This is also known as H 1 conformity or conformity to the Sobolev space. In the case of tetrahedral elements, it can be shown [58] that the sucient condition for the conformity to the Sobolev space is the fulllment of the following inequalities, p e j p s i p b : (3.31) 3.5 Variational Formulation By taking the boundary condition in Eq. 3.7 into the account, we can eliminate the surface integral terms of Eq. 3.11. Recall that the boundary conditions are the characteristics of the physical quantity of interest, which are satised by the exact solution. It must be further guaranteed that the FEM approximate solution in Eq. 3.30 fullls the boundary condition, too, a property which is not held in general. For example, in the case of polynomial basis functions, for all boundaries passing through an element, the following relation holds for points positioned on the layers, 8 > > > < > > > : D + rP + :~ n = D rP :~ n rP + = rP )D + =D Contradiction! (3.32) which could cause contradictions if the self-diusion coecient is not continues during transition from one side of the boundary to the opposite side. To address this diculty, 62 we propose the following design which allows the partial derivatives of the solution at the compartment boundaries to be discontinues and, in the meantime, maintains the smooth- ness within the compartment regions. Satisfying the boundary condition necessitates that at points located on the compartments layers which are characterized by separating two regions with dierent self-coecient values, the rst derivatives of the solution to be dis- continuous as given in Eq. 3.32. This can be achieved by violation of the Minimum Rule, provided in section 3.4, for the elements which have intersection or shared faces with com- partment boundaries. This could place serious burden on the implementation since we have to treat some elements dierently. This problem can be eliminated by forcing com- partment boundaries to have incident with the FEM mesh only at the element faces. In other words, the compartment boundaries never split any elements. To this end, violating the Minimum Rule for all elements gives the approximate solution the exibility to have jumps in its rst order derivatives at the compartment boundaries. Meanwhile, for the element faces, with no shared regions with the compartment boundaries, the boundary condition is translated to, 8 > > > < > > > : D + rP + :~ n = D rP :~ n D + = D ) rP + :~ n =rP :~ n: (3.33) Since this equation is applicable to layers in all direction, it can be concluded that the partial derivatives exist, which guarantees the smoothness of the solution. It should be noted that, with the proposed setting, the self-diusion coecients within an element must be continues which is a trivial requirement, since it is usually assumed that its 63 value is constant within a compartment. As a consequence of the above discussion, the surface integral part in Eq. 3.11 can be safely canceled out. X i Z C i P t udr = X i Z C i D i ru:rPdr; (3.34) Z P t udr = Z D(r)ru:rPdr; u(r)2V n (3.35) X j Z K j P t udr = X j Z K j D j ru:rPdr: (3.36) It should be emphasized that Eq. 3.36 could not directly be derived from Eq. 3.3, since the self-diusion coecients have discontinuities at compartment boundaries. Also the original space of function u is conned to a new space, spanned by the basis functions of all elements, which is a subspace of the function space L 2 . Substituting the approximate FEM solution into Eq. 3.36 and taking into the account the isoparametric map, we obtain, n X l=1 l (t) X j jJ K j j Z ^ K ~ j q ~ j l d = n X l=1 l (t) X j jJ K j jD j 3 X m=1 3 X s=1 3 X r=1 @ s @x m : @ r @x m : Z ^ K @ ~ j l @ s @ ~ j q @ r d; q = 1; 2;:::;n (3.37) where J K j is the Jacobean, and @m @x l the entries of the inverse of the Jacobi matrix of the isoparametric map associated with element K j , and ~ j l = ' n x 1 K j . Eq. 3.37 can 64 be represented in the form of a simple ordinary dierential equation as presented in the following relations, M = S = [ 1 (t); 2 (t);:::; n (t)] (3.38) whereS andM aren byn matrices which are referred to as stiness and mass matrices, respectively. 3.6 Initial Unrestricted Self-diusion Regime In this section, it is mathematically proven that for a pointr 0 , within a multi-compartment environment, for the small time interval [0;t ], self-diusion will have an unrestricted (free) regime. Let us consider a three-compartment environment with self-diusion coef- cients D 1 ;D 2 , and D 3 which belong to compartments 1; 2; and 3 respectively, and the goal is to evaluate the self-diusion evolution at a point within region 1, as illustrated in Fig. 3.5. Let's assume that the time interval is suciently small such that the water molecules do not have time to travel beyond distancel, whichl is the distance of the point r 0 from the nearest boundary. By denition, the propagator values will be united by zero for all locations whose distances from r 0 are larger than l, which include the regions 2 and 3. As a result, we have the following equalities, @P r 0(r;t) @t = Dr 2 P r 0(r;t); P r 0(r;t) = 0 =) 0 = 0:D: r>l; t2 [0;t " ] (3.39) 65 Figure 3.5: Three compartment Environment: A three-compartment self-diusion envi- ronment Eq. 3.39 is valid irrespective of all nite values ofD, which implies that, during the small time interval, the original PDE is independent of the self-diusion coecient values of regions 2 and 3. In other words,D 2 andD 3 can pick any arbitrary nite values including the D 1 value. That is, the self-diusion can be considered unrestricted. 3.7 Results 3.7.1 Simulation Parameters For experiments of this study, unless otherwise specied, the following parameters setting will be employed. The sample dimensions (i.e., the volume over which the PDE solution is evaluated) are set to 15 by 15 by 15 [m] along the x-, y-, and z-axes, respectively, as shown in Fig.3.8. The radii of inner cylinders corresponding to the cytoplasm regions (see Fig.3.1) are set to R in = 3 [m] and we consider thick myelin sheath layers, which results in the myelinated axons, to be with total radii of R out = 4 [m]. For hexago- nal organization of axons within the white tract, the distance between the axels of two 66 neighboring cylinders is chosen to be L = 10 [m]. The self-diusion coecients for the cytoplasm, myelin, and extracellular regions are selected asD c = 1:0 [m 2 =ms],D m = 0:5 [m 2 =ms], and D e = 1:9 [m 2 =ms], respectively. Following the same parameterization for water concentrations in [56], the water molecules concentrations are set to c = 0:85, m = 0:5, e = 0:9, and g = 0:9 for the cytoplasm, myelin, extracellular, and o-tract regions, respectively. Tetrahedral elements are utilized for construction of the sample mesh of the FEM. Each element is equipped with vertex, edge, and face basis functions. In order to relax the minimum rule of conformity to the Sobolev space, the maximum degree of polynomial functions for edge, face, and bubble basis functions are set top e = 3, p f = 4 and p b = 0, respectively. Such a setting provides each element with a degree of freedom of 14 basis functions. Approximately 40,000 tetrahedral elements were employed for discretization of the sample of interest. For solution of the rst order ODE system in Eq.3.12, the rst-order IRK method is applied. 3.7.2 Contributions of the FEM Basis Functions in the Global Approximate Solution Dening the type and number of basis functions in the approximate solution in Eq. 3.30 is of special signicance for precision of the solution, as well as the ecient allocations of memory and computational resources. Observing the distribution of the solution signal energy over dierent basis functions at dierent time instances could guide us towards better design of the FEM setups. In Fig. 3.6 typical energy distributions of experiments of this study at two dierent time instances are provided. It shows that at the timet = 1:0 [ms] the approximate solution receives contributions from all kinds of tetrahedron basis 67 functions, whereas after t = 15 [ms] the energy is compacted in the lower-order (Vertex) basis functions. This is an indication of signicance of higher-order FEMs for the transient regime, which naturally has a strong indirect impact on the steady state solution, since the one-step method for solution of the ODE is employed. 3.7.3 Equivalent membrane model The experiments conducted in this section are used to validate the accuracy of discussions in Section 3.2.3 regarding the equivalent model for a thin membrane wall, where it was emphasized that as long as the ratio D b =b is kept constant, the realizations of the self- diusion propagator outside the equivalent layer are unchanged. Two cases are presented in Figure.3.7, where two large compartments with self- diusion coecients 0.6 and 0.1 m 2 =ms for the left and right regions, respectively, are separated by a thin layer of widthb and self-diusion coecientD b . The experiments were performed with two sets of parameters as depicted in the gure. The key concept is that the ratio of D b =b is the same in both experiments. Figure.3.7 shows that the realizations of the self-diusion propagator outside the equivalent thin layers are very similar as predicted in Section 3.2.3. Then, as the width of the layer becomes very small in the limit (b! 0), the approximate layer should converge to an actual membrane wall by keepingD b =b the same and Eq.3.8 will hold. In practice, an innitesimally thin mem- brane wall with permeability constant , can be replaced with a thin layer, whose width and self-diusion coecient satisfy the =D b =b condition. 68 Figure 3.6: Distribution of solution energy: Distribution of energy over the basis functions of the solution of self-diusion PDE for a point located in the cytoplasm region. The low- level (vertex) basis functions are indexed from 1 to 7035, (Top) The distribution for t = 1:0 [ms], (Bottm) The distribution for t = 15 [ms] 69 Figure 3.7: Taner-Powel equivalent models: The equivalent model for the membrane wall. The propagator is evaluated along the direction normal to the boundaries. The origin is translated to the self-diusion start point. In the two scenarios of study, the ration D b =b is kept constant to 10. It can be seen that, outside the thin equivalent layers, the realizations of propagators stay very close together. 3.7.4 Hexagonal Bundles of Healthy Axons The results of the FEM solutions of the self-diusion propagators at dierent locations denoted by labels 1; 2;::: and 5 are illustrated in Figs.3.8 and 3.9 . Fig.3.8 species the transverse propagators -the probability of propagation in the plane perpendicular to the axon orientation- at dierent location after the time elapses by 15 [ms]. The eects of the myelin sheaths around the axon cells in characterization of the transverse propagators are evident. In other words, despite the fact that the concentration of water molecules in the myelin regions is very low and those regions have negligible direct contributions to the average propagator according to Eq. (3.21), they form the shape of the propagators in the neighboring cytoplasm and extracellular regions due to their partial permeability. In Fig. 70 3.8, the time evolutions of transverse and parallel (along the axons orientation) propaga- tors at the position 1 are provided, as well as transverse and parallel propagators at the labeled points after time t = 15 [ms]. The reason for choosing this self-diusion time is to provide the water molecules adequate time to probe the axon model microstructure. Figure 3.8: Transverse propagator: (Top row) (Left) Hexagonal organization of myeli- nated axons and the labeled points for evaluation of propagators, (Middle) Transverse propagator at position 1, (Right) Transverse propagator at position 2, (Bottom Row) (Left) Transverse propagator at position 3, (Middle) Transverse propagator at position 4, (Right) Transverse propagator at position 5 Since the propagation along the axon orientation encounters no limitation, the evolu- tion of parallel self-diusion prole always stays close to the Gaussian function. On the other hand, as a result of the partial permeability of the myelin sheaths, there are drastic deviations from the Gaussianity in the transverse propagators, as shown in Fig3.9. For validation of the PDE solution, following the discussions which led to derivation of Eq.3.18, the estimate of the self- diusion coecient for each compartment is achieved 71 by averaging the results of evaluation of Eq. 3.18 over several sample points of the com- partment region. Fig.3.10 shows the results of estimated self-diusion coecients for the cytoplasm and myelin regions, provided that their true values are D c = 0:6 [m 2 =ms], and D m = 0:015 [m 2 =ms], respectively. It can be seen that, for the myelin region, it takes more than 5 [ms] for the estimated values to fall in the valid range. The origin of this delay for getting accurate estimates from Eq. 3.18 is the time required for the water molecules to depart from the point 1 and surpass the myelin sheath. In other words, for the time intervals less than 5 [ms], the numerator and denominator of Eq. 3.18 are not well-established -almost zero- and the application of the equality in estimation of self-diusion coecients is meaningless. Also, the standard deviations of the estimates represent the precision of the computed parameters over the sampled points of the cor- responding compartment. It can be observed from the gures that the quality of the estimates become superior as the time elapses. 3.7.5 Voxel Average Propagators for Single Tract Scenarios D-MRI signals are sensitive to average values of self-diusion propagators computed over the voxels, rather than their point-wise values. To this end, we formulate the result of such averaging for the case of a single tract passing through the MRI voxel, as illustrated in Fig. 3.11. The key relation for the calculation of the average propagator over the voxel is Eq. 3.21. The equation is characterized by average propagators as well as spin magnet fractions of compartments. For computation of the compartments average propagators, discrete integration techniques are employed.That is, the PDE of self-diusion propagator 72 Figure 3.9: Parallel and transverse propagators: (Top row) (Left)Time evolution of the parallel propagator at position 1,(Right) Time evolution of transverse propagator at po- sition 1, (Bottom row) (Left)Parallel propagators at time t=15 at labeled positions in Fig.3.8 ,(Right)Transverse propagators at time t=15 at the labeled positions 73 Figure 3.10: Validation of the FEM solution: Validation of the propagator at point 1 with known true self-diusion coecient values of D c = 0:6 [m 2 =ms], and D m = 0:015 [m 2 =ms]. The results are achieved by averaging of the evaluations of Eq. 3.18 at sample points of each compartments, (Left) the cytoplasm region estimated self-diusion coecient and the corresponding standard deviation, (Right) the same quantities for the mayelin region. is solved at dierent sample point of the compartment and the results are integrated by discrete techniques. The spin magnet fractions are characterized by partial volumes and concentrations of water molecules over the compartments. To evaluate the propagators in the regions of the voxel that are not parts of the passing tract(s), an isotropic propagator model as provided in Eq. 3.15, is considered. It can be concluded from Fig. 3.11 that the anisotropic self-diusion of the myelinated ax- ons are mostly due to the cytoplasm and extracellular regions. The average propagator of the myelin sheaths exhibits weak anisotropy, which alongside the fact that the concentra- tion of the water molecules in those regions are relatively lower than other compartments of the white tracts, severely downplay the direct role of the myelin sheaths in determina- tion of the total average anisotropy of the voxel propagator. The noticeable observation 74 is a large amount of transverse propagation provided by the extracellular regions which largely contributes to increasing the average transverse propagators of the voxel. The source of this incident is the open space among the axons. The 2D average parallel and transverse propagators indicate substantial degree of anisotropy within the voxel. This is an indication of the required time for manifestation of anisotropic propagation in the white tract. In Fig.3.12 the time evolution regime for self-diusion propagator for the case of a single axon is provided. It can be seen that the prole stays symmetric at all temporal moments and it exhibits severe degree of anisotropy with the preference along the axon orientation. Once again, the eect of the myelin sheaths in hindering the transverse self-diusion process in evident. 3.7.6 Equivalent Tensor Model for Propagator Because the required time for data acquisitions in DTI is relatively small, compared to other D-MRI methods, and the computational complexity for the post processing algorithms is low, it still remains very popular for tractography of the brain white tracts, even after two decades of the DTI invention. So it is of special interest to understand and realize the extent of the tensor model performance and identify its shortcomings and situations it fails to perform well. In this section, the ideal performance of a tensor model in characterization of a voxel average self-diusion propagator for a single tract, which was computed by our developed FEM method, will be evaluated. The reason we used the term "ideal" is that in real DTI technology, the self-diusion prole is modeled by a tensor. However, in estimation of the 75 Figure 3.11: Average propagator: (Top row)(Left)The average propagator of the cyto- plasm region in 3D, (Middle) The average propagator of the myelin region in 3D ,(Right) The average propagator of the extracellular region in 3D, (Middle row)(Left) A single tract passing through a voxel and the reference coordinate system, (Middle) The aver- age transverse (in x-y plane) propagator of the voxel, (Right) The average parallel and transverse propagator of the voxel in 2D along the z- and y-axes, (Bottom row) The average propagator in 3D with tract volume fraction T vf = 0:75. All 3D propagators are evaluated on a shell of radius 8 [m]. The times elapsed for all gure are t = 15 [ms]. The isotropic self-diusion coecient for o-tract regions of the voxel is set to D g = 2:0 [m 2 =ms]. 76 Figure 3.12: Time evolution of propagator for the single tract case: Time evolution of the self-diusion prole at center of an axon, (top row) from left to right, propagator proles at t = 4; 5; 6 [ms],(second row) propagator proles at t = 7; 8; 9 [ms], (third row) propagator proles at t = 10; 11; 12 [ms], (bottom row) propagator proles at t = 13; 14; 15 [ms] 77 tensor model parameters, the exact value of the propagator are not directly measurable, and they have to be estimated from the DTI signal which is subject to errors and imper- fectness. Therefor, in real DTI technology, the tensor model characterizes a noisy version of the underlying self-diusion prole. As opposed to the real DTI measurements, our method makes it possible to provide clean propagator data to evaluate the performance of a tensor model under ideal imaging conditions. To make the tensor parameter esti- mation task easier, we consider the single tract scenario in Fig.3.11 in which the average propagator is symmetric around the z- axis and it has strong directional preference along that axis. Therefore, the equivalent tensor model principal eigenvector would be aligned with the z- axis and the eigenvalues corresponding to minor eigenvectors would be equal. Note that, this proper choice of coordinate reduces the total tensor parameters to be estimated, from six to three parameters. That is, D eq = 0 B B B B B B B @ D ? 0 0 0 D ? 0 0 0 D k 1 C C C C C C C A (3.40) 78 In DTI, it is assumed that the self-diusion prole within the D-MRI voxel is governed by a single equivalent tensor model, and consequently the propagator prole within a voxel is provided by the following closed-form equation, P (r;t) = (4t) 3=2 (D k D 2 ? )exp 0 B B B B B B B B B B B B B B B B B B B B B @ r T 0 B B B B B B B @ 1=D ? 0 0 0 1=D ? 0 0 0 1=D k 1 C C C C C C C A r 4t 1 C C C C C C C C C C C C C C C C C C C C C A (3.41) The tensor characterizing parametersD ? andD k are estimated by tting the true values of the voxel average propagator acquired from our FEM method, into equation 3.41. For this study, we employed the minimum mean square estimation technique (MMSE) for parameter estimations. This will involve non-linear estimation which was optimized by the gradient descent method. The resultant estimated parameters of the tensor and its prole of the voxel average self-diusion are provided in Fig.3.13. Fig.3.13 shows that even for a simple case of the voxel with a single tract passing through, the tensor model may not be an appropriate choice. Its precision might be adequate for detection purposes such as the estimation of the ber orientation but its capability for simulation of the diusion process in the voxel is questionable. For example, in the literature, there are a number of work in which the tensor model was employed to simulateD-MRI phantom for validation of tractography methods. Fig.3.13 illustrates an apparent error which could easily propagate into the simulated signals by such methods. 79 Figure 3.13: Equivalent tensor model: (left) Voxel average self-diusion prole for a single tract in Fig.3.11 provided by our FEM method, (right) Equivalent estimated tensor with estimated parameters of D ? = 0:1 [m 2 =ms] and D k = 2:0 [m 2 =ms] This raises a serious concern about the validity of the tensor model, especially in the applications in which the tensor is treated as a generative model. 3.7.7 Voxel Average Propagator for Two Crossing Tracts Eq. 3.21 is scalable to accommodate conditions in which two or more white tracts intersect in the MRI voxel. Comparing to the single tract case, the extra information needed is the volume fraction associated with each tract. Fig.3.14 shows the outcomes of the studies for two scenarios with dierent tract volume fractions. . Quantication of the eects of individual tract volume fractions on the net voxel propagator, is one of the highlights of the developed method, which is illustrated in the mentioned illustration. It will be illustrated that, the average propagator of two crossing tracts, is not simply a linear weighted average of individual propagators, which negates the constitutions of the mixture models for representation of ber crossing situations [15]. In Fig.3.14 it 80 can be seen that the isotropic regions have reduce the degree of anisotropy of the voxel propagator. Figure 3.14: Voxel average propagator for two crossing tracts: (Top) The average prop- agator for crossing tracts with volume fractions T vfz = 0:45 and T vfy = 0:35, (Bottom) the same analysis for T vfz = 0:25 and T vfy = 0:50 3.7.8 Demyelinated Tracts The exibility of the FEM in representation of dierent complex boundary conditions make it possible to study the eects of nerve degeneration and morphological abnormali- ties on the self-diusion prole through inducing them in the axon model and observe the results on the simulated outcomes. Among such degenerations is the demyelination case 81 in which the myelin sheath layers around the axon cells are partially or totally destroyed. Such phenomenon can be seen in neurodegenrative diseases like Multiple Sclerosis in which the demyelination hinders or blocks the transmission of action potentials along the axons due to electric leakage of the message signals. For the current work, the partial demyelination of axons in the hexagonal model of the white tracts, is considered. As illustrated in Fig.3.15, in the scenario under study, some part of the myelin sheath along the z- axis is destroyed. As a result of this nerve degeneration, the propagator prole became severely lateral towards the broken sheath region. For comparison, the prole of the healthy nerve is also provided. It can be seen that, the parallel propagation along the z- axis underwent rigorous reduction as compared to the healthy case which could severely reduce the degree of apparent anisotropy in the observed D-MRI signal. Figure 3.15: Eects of demyelination on propagator: Eect of axon demyelination on the self-diusion prole (Left) A hexagonal organization with partial demyelination in the central axon, (Middle) the 3D self-diusion prole in the point marked by the red arrow, (Right) the self-diusion prole for the same location for healthy axons In Fig.3.16, the time evolution of the self-diusion propagator at a point located in the center of a demyelinated axon, as illustrated in Fig. 3.15, is provided. 82 Figure 3.16: Time evolution of propagator for the demyelination case: Time evolution of the self-diusion propagator in the center of the demyelinated axon in Fig.3.15. (Top row) from left to right, propagator proles for time instances, t = 4; 5; 6 [ms], (Second row) propagator proles for time instances,t = 7; 8; 9 [ms], (Third row) propagator proles for time instances, t = 10; 11; 12 [ms], (Bottom row) propagator proles for time instances, t = 13; 14; 15 [ms] 83 Fig. 3.16 shows that at the beginning, due to the eect of large self-diusion coecient of the extracellular region, the lateral component of the propagator has a very fast growth, since there is no myelin sheaths barrier to slow it down. However, as the time passes by, as a result of partial permeability of the surrounding axon myelin sheaths, the fast evolution pace of the lateral component is hindered, but it still has severe asymmetry, after t = 15 [ms]. This knowledge could be of special interest for designing of D-MRI pulse sequence for diagnosis of demyelination. It can be seen that, the more time elapses, the less lateral the self-diusion prole becomes. Therefor, if the self-diusion time in imaging is too long, it may become dicult to observed the lateral prole as a signature of demyelination. 3.8 Discussion One of the highlights of the developed method is its capability to quantitatively link the white tract microstructure parameters such as axon diameters, inter axon distances, and myelin sheath thicknesses, to self-diusion proles. In the computations, the modular eects of constituting components of the white tracts are factored. It is possible to accommodate a vast range of variability in the white tract components irrespective of the level of geometrical complexities they exhibit. Thus, it is possible to study the behavior of individual components and their contributions to the aggregate self-diusion propagators. For example, in Fig. 3.11, it is evident that for long self-diusion time, the myelin sheaths substantially contribute to parallel and transverse propagation. However, since the number of water molecules in those regions is relatively less than the other component of the white tracts, their eect on the average propagator is insignicant 84 by Eq. 3.21. However, their roles are prevalent in formations of other compartment self-diusion. InD-MRI data acquisition, it is well understood [16] that for unraveling the underlying microstructure of the brain white matters, it is necessary to give the self-diusion process enough time to probe the white tracts structure dimensions. One of the prominent aspects of the developed technique is quantication of this phenomenon. For example, by examining the time evolution of parallel and transverse propagators at point 1 in Fig. 3.9, it can be observed that the propagator shows a signicant amount of anisotropy at timet = 15 [ms], contrary to thet = 4 [ms] case. This is more apparent from the average parallel and transverse propagators of the voxel in Fig. 3.11. This could potentially lead to more eective designs for D-MRI pulse sequences by customization to the prior knowledge of characteristics of microstructure of interests. It is also promising to evaluate the consequential eects of tract structural parameters on the self-diusion proles, as illustrated in Fig. 3.15 for the case of demyelinated axons. Inducing variations in parameter values such as myelin sheath thicknesses, axon diameters, and axon spacing, and evaluating their impacts on the propagator, could help uncover distinguishing D-MRI signatures of neuro-degenerative diseases. This might be superior over direct measurements of the propagators from real D-MRI signals of the patients, since the ground truth is already available. In direct methods, the signature for a specic nerve abnormality is sought through data mining inD-MRI signals acquired from the subjects diagnosed with that abnormality. There are a number of issues associated with this methods including the lack of adequate local resolution and also validations. 85 One application of the presented work is validity assessments of propagator parametric and non-parametric models for representations of dicult tract conditions. It can be observed from Fig. 3.11 that the average prole of the self-diusion of a single tract over the voxel closely resembles the Gaussian function. Hence, the famous tensor model, commonly utilized in DTI methodology, is of adequate precision. However, the unilateral extension of the propagator prole towards the demyelinated region in Fig. 3.15 raises serious concerns over the accuracy of such models for studies of similar situations such as multiple sclerosis. This observation is alarming for being more cautious about applications of the anisotropy measures such as Fractional Anisotropy (FA) which are based on the validity of the tensor model, in diagnoses of demyelination related illnesses [22, 73]. The partial volume problem due to limited spatial resolution is an eminent issue in the D-MRI technology, especially in the brain regions where two or more tracts cross one another. There have been a great amount of eorts in the literature [15] to account for such situations through a variety of modeling schemes. For the two-tract situation, it can be concluded from Eq. 3.21 that the voxel average propagator consists of a linear combination of three decomposing components, P v (R; ) =f An1 P An1 (R; ) +f An2 P An2 (R; ) +f ISO P ISO (R; ); (3.42) where P An1 (R; ) and P An2 (R; ) are average anisotropic propagators due to the pass- ing tracts through the voxel, and P ISO (R; ) is the average isotropic component, cor- responding to the net eect of the o-tract isotropic regions. Fig. 3.17 illustrates such 86 a decomposition for the two crossing-tract cases as presented in Fig. 3.14. This phe- nomenon can be examined for three or more crossing tracts and, as a general rule, for representing the average propagator of the voxel, with n crossing tracts, n + 1 indepen- dent component are required. This observation emphasizes the signicance of precisely considering the isotropic component, especially for the cases in which tracts have small volume fractions. It could be a perceptible source of errors for modeling schemes, should we ignore or downplay the contribution of the isotropic term. Figure 3.17: Decomposing components: Decomposing components of the voxel average propagator of the crossing tracts in Fig. 3.14. In addition to two anisotropic propagator proles for each tract, one more isotropic prole is provided to account for the o-tract region of the voxel. 3.9 Conclusion In this chapter, we developed and implemented the tetrahedral FEM solver to address the solution of the self-diusion propagator within the hexagonal array of cylinder model for the brain white tracts. The extent of the developed method encompasses far more complex geometries for the white tracts. It provides the foundation for numerical simulations for studies of the white tracts under a variety of conditions. In our view, the current work can potentially be further explored in three major dimensions. The rst application is 87 to leverage the information learned from simulation of self-diusion in the brain tissues, towards better pulse sequence design of D-MRI scanners, provided that the time and space evolutions of self-diusion is now quantitatively available. An example of such possibility is assessment of the required time to observe anisotropic propagation within a voxel, which can straightforwardly be determined by the developed tool. Reverse engineering of neuro-degenerative diseases could be regarded as another likely domain for the employment of our method. As exemplied by the case of demyelination, such a situation can simply be simulated in order to nd discriminating signatures for diagnostic purposes, instead of data mining in real D-MRI signals. Finally, simulating the D-MRI data from the calculated self-diusion propagators of known white tract geometries, could serve as the ground truth for the validation of dierent tractography methods. Through these procedures, many tracts conditions and D-MRI artifacts can be incorporated into the simulated signal, to evaluate their eects on tractography. As previously stated, having an ad-hoc approach for modeling the behavior of self-diusion prole (e.g. tensor model), lacks the required precision for simulating the self-diusion process and the generated signal would not resemble real D-MRI data. By having a more realistic model for self-diusion process, more realistic D-MRI signal can be simulated, which gives the resultant ground truth more merit, in evaluation of performance of tractography methods. 88 Chapter 4 : Theory for Practical MRI and Diusion MRI Simulation Most existing MRI and diusion MRI signal modeling and simulation formulations are derived under perfect imaging conditions. However, these conditions are often violated in real world MRI scanners which makes those relations prone to errors. The impact of the formulation inaccuracy becomes more serious when we consider their consequences on the performance of many diagnostic techniques. In this work, we extend those theories to accommodate real imaging and intrinsic parameters. 4.1 Physics of MRI and Diusion MRI The source of MRI signals is the variation of the ensemble eld of magnetic particles in the scanner receiver coil. For medical imaging applications, these particles are often protons of the hydrogen molecules in the tissue water. As a result of the revolution of a particle (also known as the spin magnet) around its axis, a small magnetic eld is induced in the surrounding space with the direction along the revolution axis as illustrated in Fig. 4.1. With the precise coordinations of dierent magnetic elds in MRI scanners, the small 89 elds of the spins can be added to establish a sensible magnetic eld in the MRI machine receiver coil. Figure 4.1: A spin magnet and its magnetic eld An MRI scanner equipped with a main strong magnetic eld which is static, and a number of temporally changing magnets. The conventional coordinate system in MRI scanners is shown in Fig. 4.2. Often, the direction of the main magnetic eld, denoted byB 0 , lies in the direction of the z-coordinate. Figure 4.2: The coordinate system adopted in an MRI machine. 90 If a spin magnet is under a strong external magnetic eld, its revolution axis starts a precession motion around the external eld as shown in Fig. 4.3. The angular velocity of the precession, denoted by ! r is proportional to the magnitude of the external eld in form of [31] ! r = B 0 ; (4.1) where is the gyromagnetic constant coecient, which is the characteristic of the particle. This normalized angular velocity is also known as the Larmor frequency. Figure 4.3: Spin precession in the presence of an external magnetic eld. In the presence of an external main eld, it is easier to study the average magnetic behavior of spin magnets in the volume unit, which is referred to as local magnetization vector and denoted by ~ M(~ r). It represents the average magnetic eld of the spin magnets 91 in an innitesimally small volume at position ~ r. It can be shown that the magnitude of this vector is proportional to the external eld [31] as M 0 (~ r) =C 0 (~ r)B 0 (~ r) T ; (4.2) where C is the particle constant, 0 (~ r) is the local density of particles, and T is the absolute temperature. It is useful to decompose local magnetization vectors into two components which are parallel and perpendicular to the main eld orientation. The two components are known as parallel and transverse components as shown in Fig. 4.4 and denoted by M k and ~ M ? , respectively. Figure 4.4: Parallel and transverse components of local magnetization. The transverse magnetization rotates in the x-y plane with the Larmor frequency. In MR experiments, the application of a short RF pulse, which rotate in the transverse plane with the Larmor frequency, tips the local magnetization vectors into the transverse plane and causes them to have equal phases right after its application, as shown in Fig.4.5. After the RF pulse application, the magnetization vectors start rotating in the x-y plane with 92 the Larmor frequencies. The motion of the local magnetization vector can be predicted by the solution of the Bloch Equation given below [31]: d ~ M(~ r) dt = ~ M(~ r) ~ B ext (~ r) + 1 T 1 (~ r) (M 0 (~ r)M k (~ r))^ z 1 T 2 (~ r) ~ M ? (~ r); (4.3) where ~ B ext (~ r) is the local value of the net external magnetic eld. Figure 4.5: RF pulse application: Application of a rotating RF pulse with the Larmor frequency tips the spins into the transverse plane. It can be shown that if the coecients T 1 and T 2 in the Bloch equation become very large, its solution is a rotating vector in the transverse plane with the Larmor frequency. That is, for such ideal conditions, the local magnetization vectors remain in phase and lie in the x-y plane rotating with the Larmor frequency. However, in practice, the nite values of those coecients violate such ideal conditions and the phenomena are called relaxations. The rst phenomenon of interest, known as theT 1 relaxation, represents the tendencies of spin magnets to be aligned with the external magnetic eld, after being tipped into the transverse plan by the RF pulse. As a result, the parallel components of the magnetization 93 vectors, which are zero right after the RF pulse, have exponential growths to their original M 0 value; namely, M k (t)(~ r) =M 0 (~ r)(1 exp t T 1 (~ r) ): (4.4) Another type of relaxation is due to the interaction of the magnetic elds of spin magnets with each other which causes gradual de-phasing. As illustrated in Fig. 4.6, such incidents weaken the net local transverse magnetization vector. As a result, the local transverse magnetization decays exponentially at a rate of 1 T 2 [31]: M ? (~ r;t) =M ? (~ r; 0) exp t T 2 (~ r) ; (4.5) whereM ? is the magnitude of transverse magnetization vector ~ M ? , and the time reference is the application of the RF pulse unset. The T 2 relaxation eect is an intrinsic property of spin magnets and cannot be compensated by imaging parameters. There is another type of de-phasing, which is the consequence of inhomogeneity in the main external eld. This is known as the T 0 2 relaxation. Even in modern scanners, establishing a strong uniform external eld is very dicult and spatial variations always exist. Spin magnets at dierent positions are under a dierent external eld, which causes them to have unequal Larmor frequencies and gradually de-phase. As a result, there is another exponential decay in the magnetization vector at a rate of 1 T 0 2 . Fortunately, this kind of relaxation can be compensated by the application of the Spin Echo (SE) pulse 94 sequences, which will be explained later on. The combined de-phasing types of T 2 and T 0 2 relaxations can be represented by a single exponential with an unied rate of 1=T 2 = 1=T 2 + 1=T 0 2 : In practice, theT 1 value is much larger than theT 2 value for a lot of tissue material and, if the experiment is not designed for T 1 measurement, theT 1 eect can often be ignored. In contrast, the situation is very dierent for the T 2 parameter. Since its magnitude is usually smaller or in the same range of the MR data acquisition time, it could result in substantial errors in the experiment results if its eect is not accounted for properly. Figure 4.6: De-phasing: Illustration of spins' de-phasing due to the T 2 relaxation. As stated before, the ensemble eect of all magnetization vectors creates a magnetic ux in the MRI scanner receiver coil. There are three sources of temporal variation in the ensemble eld; namely, the rotating transverse components of the local magnetization vectors, the T 1 relaxation, and the T 2 relaxation. The temporal variation of the ux 95 induces an electromotive force (emf) in the receiver coil. Mathematically, we have the following two equations [31]: M (t) = Z ~ rec (~ r) ~ M(~ r;t)d~ r; (4.6) emf = d dt M (t); (4.7) where ~ rec (~ r) = ? (~ r) h cos ~ rec (~ r) ^ x +sin ~ rec (~ r) ^ y i + k (~ r) ^ z; (4.8) is the magnetic eld induced by the receiver coil if it is powered by one unit of electric current. It is a characteristic of the coil and temporally constant during MR experiments. Eq.(4.6) suggests that the magnetic ux is a function of magnetization vector phases at dierent positions. The magnetization vector phase value is determined by the angular velocity of the magnetization vector (the local Larmor frequency) and its initial phase value. Since the RF pulse application sets all spin magnets to have the same initial phase, the value of the ux is dened by the Larmor frequency according to Eq.(4.1) and, hence, the local value of the external eld. This property, if appropriately employed, can be used for the frequency encoding of dierent positions of the samples under study. The methodology for utilizing such encoding includes the application of dierent magnetic eld gradients due to the spatial variation in the main magnetic eld. In MRI scanners, the variation often happens linearly along a direction. Fig. 4.7 illustrates a gradient 96 which occurs along the x-axis. One or more gradients may exist in a scanner which may or may not temporally coincide. We can select the appropriate gradient parameter that is specic to local characteristics of the tissue in the MR images. Figure 4.7: A linear gradient along the x-axis. Depending on the gradient pulse parameters such as their magnitudes, durations, and temporal coordinations and the application of the RF pulses, dierent pulse sequences can be dened for MR data acquisition. The Spin Echo pulse sequences are widely utilized in dierent MR experiments, since they can eectively eliminate the T 0 2 relaxation and recover the signal loss. There are a number of variants for such pulse sequences. One of the common designs is shown in Fig. 4.8. In the spin echo sequences, the second RF pulse reverses the phases of spin magnets, causing them to re-phase at echo time T E . The subsequent derivations in this study are based on the assumption that the spin echo sequences are used in data acquisition. 97 Figure 4.8: The spin echo pulse sequence. For perfect data acquisition with no relaxation eects, the time domain signal received in the MRI machine coil can be written as [31] Signal(t)/ Z M 0 ! r (~ r) ? (~ r) sin ~ B rec (~ r) (~ r;t) T E (~ r) d~ r; (4.9) (~ r(x;y;z);t) = Z t T E (B 0 +G x x)dt; (4.10) where T E (~ r) is the local magnetization phase at the echo time T E , and the application of the rst RF pulse occurs at t = 0. The received signal in Eq.(4.9) is multiplied by the cosine and sine demodulators to form the complex demodulated signal. The following equation is the complex signal after simplication, S(t)/ Z M 0 ! r (~ r) ? (~ r) exp i !r (~ r)t+(~ r;t)+ T E (~ r) ~ rec (~ r) ; (4.11) / Z eff (~ r) exp i h R t T E Gxxdt+4Gyyy i d~ r: (4.12) 98 In Eq.(4.12), eff (~ r) is called the eective spin density which is dened as eff (~ r) = M 0 ! r (~ r) ? (~ r) exp i ~ rec (~ r) : (4.13) The k-space variables for 2D MRI reconstructions are dened as k x = 2 Z t T E G x (t)dt; (4.14) k y = 2 Z 4G y (t)dt; (4.15) ~ k =k x ^ x +k y ^ y: (4.16) Plugging Eqs. (4.14) and (4.16) into Eq. (4.12), we obtain the fundamental formula for the MR complex signal, which also is known as the k-space signal, in form of S( ~ k) / Z eff (~ r) exp i2 ~ k~ r d~ r: (4.17) Eq. (4.17) shows that the k-space signal is the Fourier transform of the eective local spin density. The application of the inverse transform results in the local characteristics of the sample tissue. That is, we have eff (~ r) / Z S( ~ k) exp +i2 ~ k~ r d ~ k: (4.18) Eq. (4.18) oers the main relation for magnetic resonance image reconstruction. The common practice is to sample thek-space signal received in the scanner coil, and the MR image is reconstructed by applying the inverse Fourier transform to the sampled data. It 99 should be re-emphasized that, under ideal data acquisition condition with no relaxations, thek-space signal and the MR image are perfect Fourier pair, a property which does not hold under sever relaxation eects. In practice, the sampling and transforming procedures are performed in the discrete domain through the fast Fourier transform (FFT). Besides, as innite sampling of the k-space signal is impractical, windowing often causes spatial blurring in the MR images. In diusion MRI, in addition to the gradients that reveal the local spin density within the tissue, a pair of additional strong diusion gradients are employed in a specic direc- tion, in order to phase-encode the moving spin magnets during the diusion time. Fig. 4.9 shows the spin echo pulse sequence for diusion MRI data acquisition. After the application of the rst gradient, we can spatially phase-encode the sample spin magnets, and after the second pulse application, those phase shifts are perfectly recovered if the spin stays stationary during diusion time4. Any spin motions during this time period cause residual phase shifts which induce a signal decay in the scanner coil and is treated as a contrast for quantifying the local motions in the direction of the applied gradient orientation. If~ g + is a vector representing the rst diusion gradient magnitude and direction, the q-space variable is dened as ~ q = 2 ~ g; (4.19) 100 Figure 4.9: The spin echo pulse sequence with diusion gradients. Alternatively, the k-space signal in the scanner coil for innitesimally narrow diusion gradient pulses ! 0 is given by S( ~ k;~ q) = Z ~ r eff (~ r) exp i2 ~ k~ r Z ~ R P ~ r (~ r + ~ R;4) exp i2(~ q+ ~ k) ~ R d ~ R d~ r: (4.20) In practice, the phase shifts caused by the diusion gradient is much larger that those due to other gradients. That is, ~ q + ~ k~ q: This fact makes Eq. (4.20) a perfect Fourier transform. The MR image can be recon- structed by applying the inverse Fourier transform to S( ~ k;~ q). The result is the local eective density weighted by the local tendencies of spin magnets to move along the diusion gradient. 101 4.2 Generic Imaging Framework in Practical MRI/Diusion MRI The MRI and Diusion MRI signal generation formulations presented in the last sec- tion are valid only under ideal imaging conditions. For example, there is no relaxation with unlimited k-space sampling , which result in a perfect Fourier pair between the k- space signal and the eective density of spin magnets. In the case of Diusion MRI, the derivation was under the assumption that the diusion gradient was very narrow, and the Fourier relations hold if the magnitude of ~ k was negligible as compared to that of ~ q, which is dicult to fulll under the narrow diusion gradient condition. In practice, all the above simplifying conditions can be easily violated, and the results of the last section might be partially or completely false. In this section, we intend to relax these constraints and extend the formulations to include more realistic imaging and intrinsic parameters such as wide diusion gradients and the T 1 and theT 2 relaxations. It should be emphasized that our subsequent formulations are derived under the spin echo pulse sequence assumption, but the same methodology can be easily extended to other pulse designs. 4.2.1 Practical MRI Signal Generation Inclusion of relaxation parameters have already been achieved for constant T 1 and T 2 parameters in [31]. However, complex tissues are composed of dierent material exhibiting dierent relaxation eects. As a result, the relaxation parameters become functions of the position so that methods described in [31] do not provide accurate results. 102 We should point out that the fundamental equations, i.e., Eqs. (4.6) and (4.7), for the received signal in the scanner coil are unchanged regardless of the values of relaxation parameters. Since the derivation is done for the spin echo sequences, we derive the formu- lation of the modiedk-space signal during the data acquisition phase as depicted in Fig. 4.8. Often, in order to reduce the noise eect and suppress the T 2 relaxation, the entire k-space is not sampled in one round yet multiple measurements are performed through repetitions in MRI data acquisition experiments. The repetition period is denoted by T R . It is straightforward to show that the magnitudes of transverse local magnetization vectors right after the application of the rst RF pulse is given by [31], M ? (~ r;t = 0) = M 0 (~ r) 1e T R =T 1 (~ r) : (4.21) In the spin echo sequences, at echo time t =T E , the eect of T 0 2 relaxation is completely compensated and in the absence of external eld gradients, the magnitudes of transverse local magnetization vectors suer exponential decays as given below: M T E (~ r) =M ? (~ r;t = 0) exp T E =T 2 (~ r) : (4.22) It is known that, before and after the echo time, the de-phasing exponential decay is governed by the following time constant: 1 T y 2 (~ r;k x ) = 8 > > > < > > > : 1 T 2 (~ r) 1 T 0 2 (~ r) t<T E or k x < 0; 1 T 2 (~ r) + 1 T 0 2 (~ r) t>T E or k x > 0: (4.23) 103 Now, if the data acquisition is performed under the eects of G x and G y gradients as shown in Fig. 4.9, the transverse magnetization vectors undergo phase shifts, and the governing equation around the echo time is given by: ~ M ? (~ r;t) = exp (tT E ) T y 2 (~ r;kx) M T E (~ r) exp i((~ r;t)+ T E (~ r)) ; (4.24) (~ r;t) = Z t T E (B 0 +G x x)dt; (4.25) T E (~ r) = (G y y y +B 0 T E ); (4.26) where the exponential decay due to the T y 2 relaxation is considered in Eq. (4.24). By substituting Eq. (4.24) into Eq. (4.6), the magnetic ux can be computed as M (t) = Z ? (~ r)M T E (~ r) exp (tT E ) T y 2 (~ r;kx) cos (~ r;t) + T E (~ r) ~ rec (~ r) d~ r; (4.27) and, according to Eq.(4.7), its time derivative is the signal received in the scanner coil in form of Signal(t) / Z ? (~ r)M T E (~ r) exp (tT E ) T y 2 (~ r;kx) = 1 T y 2 (~ r;kx) cos (~ r;t)+ T E (~ r) ~ rec (~ r) +!r sin (~ r;t)+ T E (~ r) ~ rec (~ r) : (4.28) Applying the following changes of variables, 1 T y 2 (~ r;k x ) = (~ r;k x )! r (4.29) sin(~ r;k x ) = (~ r;k x ) p 1 +(~ r;k x ) 2 (4.30) 104 along with employment of trigonometric properties, we can obtain the receiver coil signal as Signal(t) / Z ? (~ r)M T E (~ r) exp (tT E ) T y 2 (~ r;kx) ! r q 1 +(~ r;k x ) 2 sin (~ r;k x ) (~ r;t) T E (~ r) + ~ rec (~ r d~ r: (4.31) Multiplying with the sine and cosine demodulators with the Larmor frequencies and performing low-pass ltering, we can express the complex MR signal received in the scanner coil as S(t)/ Z ? (~ r)M T E (~ r)e (tT E ) T y 2 (~ r;kx) ! r q 1 +(~ r;k x ) 2 (4.32) exp i h R t Te G x xdt +G y y y +(~ r;kx)+ ~ rec (~ r i d~ r: If the proportionality coecient in Eq. (4.33) is , and the modied eective spin density is dened as (~ r;k x ) = ? (~ r)M T E (~ r) exp (tT E ) T y 2 (~ r;kx) ! r q 1 +(~ r;k x ) 2 exp i (~ r;kx)+ ~ rec (~ r) ; (4.33) we end up with the following compact form for the k-space signal: S( ~ k) = Z (~ r;k x ) exp i2 ~ k~ r d~ r: (4.34) which accounts for the relaxation eects as well. 105 Eq. (4.34) is very similar to its counterpart in Eq. (4.17). However, the signicant dierence is that the Fourier relationships no longer hold between S( ~ k) and the modied spin density (~ r;k x ). In other words, applying the inverse Fourier transform to the k- space signal will no longer result in the modied spin density. The problem arises when the MRI community still reconstructs the MR image with the inverse transform. That is, the MR image, denoted by math (~ r), is reconstructed by math (~ r) =F 1 h S( ~ k) i : (4.35) Consequently, we have ^ S( ~ k) =F math (~ r) ; (4.36) where ^ S( ~ k) may dier from S( ~ k) in general according to the Fourier convergence theo- rem. The "math" subscript for the reconstructed MR image is to specify that this is a mathematical concept and has no clear relation with the modied spin density which is a physical quantity. Thus, it is interesting to nd the link between them. To realize the physical interpretation of the MR image under the relaxation eects, we multiply Eqs. (4.34) and (4.36) by 1 and integrate over the entire space, which yields Z S( ~ k)d ~ k = Z Z (~ r;k x ) exp i2 ~ k(~ r) d~ rd ~ k; (4.37) Z ^ S( ~ k)d ~ k = Z Z math (~ r) exp i2 ~ k(~ r) d~ rd ~ k: (4.38) 106 The left-hand-sides of Eqs. (4.37) and (4.38) are equal to 1. Subtracting both sides of those equations and simplifying the outcome, we get the following relations: Z Z ~ r (~ r;k x ) math (~ r) exp i2 ~ k(~ r) d~ r d ~ k; (4.39) = Z Z ~ k (~ r;k x ) exp i2 ~ k(~ r) d ~ k math (~ r)(~ r ~ r 0 ) d~ r; (4.40) Z ~ r h 4 ~ k (~ r) + 12C 1 (~ r) + 6C 2 (~ r) math (~ r) i d~ r: (4.41) Eq. (4.41) is the consequence of Eq.(4.40) by replacing (~ r;k x ) with its rst order piece- wise polynomial approximation, (~ r;k x ) 4 ~ k (~ r) + 12C 1 (~ r)jk x j + 6C 2 (~ r); (4.42) C 2 (~ r)/C 1 (~ r): (4.43) As an immediate result of Eq. (4.41), we have math ( ~ r 0 ) 4 ~ k ( ~ r 0 ) + 12C 1 ( ~ r 0 ) + 6C 2 ( ~ r 0 ): (4.44) For the spin echo sequences with strong signal recovery, we have 4 ~ k ( ~ r 0 ) 12C 1 ( ~ r 0 ) + 6C 2 ( ~ r 0 ); (4.45) and, consequently, math ( ~ r 0 ) 4 ~ k ( ~ r 0 ): (4.46) 107 Since the spin density at position ~ r 0 , changes with time, the averaging operation causes blurring in the nal image. It should be noted that this is dierent from the spatial blurring due to the limited k-space sampling. Under no relaxations, the spin density is not a function of time and no blurring occurs in this regard. In the data acquisition process using real scanners, we only receive the k-space signal in the scanner coil, which is not concerned with the way it is generated. However, in the simulation, it is necessary to quantify S( ~ k) (or alternatively ^ S( ~ k)) through the computa- tion of Eq. (4.34) (or Eq.(4.36)) and applying the inverse Fourier transform in the MR image. The domain of integration includes the entire or part of the sample of study with the dimension in the range of several inches. On the other hand, a biological tissue often consists of a lot of compartments with complicated micro-geometries in the range as small as a micron. As a result, performing a direct spatial integration such as the one in Eq. (4.34) is practically intractable. We will show that Eq. (4.46) is of great importance to compartmentalize the integration and reduce the complexity tremendously. In MRI technology, the k-space is evaluated in the discrete domain. The k-space signal is sampled at the 1=L interval where L is the image length or the eld of view (FOV). The sampling is not innite and limited to a window with width W . As a result, the image is convolved with a sinc function r 1=W with the rst zero crossing at 1=W . An example of such an image reconstruction process for the 1-D scenario is illustrated in Fig. 4.10. Existence of the spatial convolution reduces the resolution of the reconstructed image to 1=W along each dimension. Therefore, it is sucient to reconstruct the image in a discrete fashion with one sample within an MRI voxel whose dimension length is 108 1=W . Adding more samples to the image will not include any extra visual information to the image. Figure 4.10: 1-D MR image reconstruction: Illustration of the discrete 1-D MR image reconstruction. If the spatial domain is partitioned to N small voxels of length 1=W , denoted by v i , Eq. (4.36) can be reformulated by the following approximate relation, ^ S( ~ k) N X i=1 Z v i math (~ r) r 1=W (~ r) exp i2 ~ k~ r d~ r; (4.47) N X i=1 Z v i avr (~ r) exp i2 ~ k~ r d~ r; (4.48) avg (~ r) = math (~ r) r 1=W (~ r); (4.49) where is the convolution operator. By plugging Eq. (4.46) into Eq. (4.49), we have avg (~ r) ~ k (~ r); (4.50) (~ r) = avg (~ r) r 1=W (~ r): (4.51) 109 By substituting Eq.(4.50) into Eq.(4.47), we can get ^ S( ~ k) N X i=1 Z v i avg (~ r) exp i2 ~ k~ r d~ r: (4.52) In its analog form, Eq. (4.52) does not provide any computational saving. As stated previously, MRI signal processing is performed in the discrete domain. To conform to the Nyquist theorem, certain criteria must hold for the k-space signal and image sampling [31]. For the 2D cases, if the dimensions of the image along the x- and y-axes are L X and L y , and the image has N x and N y samples along those directions, respectively, the measured discrete k-space signal is given by ^ S m (p4k x ;q4k y ) = 4x4y Nx X m=1 Ny X n=1 ( mL x N x ; nL y N y ) exp i2( mp Nx + nq Ny ) ; (4.53) p = 1; 2;;N x and q = 1; 2;;N y ; where4x = 1=W x and4y = 1=W y are dimensions of the MRI voxels (W x and W y are k-space sampling intervals for the x- and y-axes, respectively). In Eq. (4.53), (~ r) is sampled at the center of each voxel, which represents the result of weighted averaging of ~ k (~ r) with the sinc function located at the center of the voxel. Since the voxel size is very small and the tissue region conned to it is approximately uniform, the weighted averaging with thesinc function is equivalent to uniform averaging over the voxel volume. In other words, an alternative relation for the computation of (~ r) at the center of voxel v i can be derived as (~ r c i ) = avg ~ r2v i ~ k (~ r): (4.54) 110 Eq. (4.54) is the key relation for reducing the k-space computation complexity. A careful inspection on Eq. (4.33) reveals that, for an uniform receiver coil magnetic eld (which is an accurate practical assumption), the value of ~ k (~ r) is almost constant for a specic composing material (e.g. fat) of the tissue. Therefor, if the spatial average of ~ k (~ r) in a volume unit of material j is denoted by j , the average spin density of MRI voxel i can be computed by (~ r c i ) = N c;i X j=1 f j;i j ; (4.55) wherec j is the compartment volume of material j conned to voxeli,N c;i is the number of compartments in the voxel, and f j;i = c j v i is the volume fraction of material j. 4.2.2 Practical Diusion MRI Signal Generation In this section, a theoretical framework for the formulation of diusion MRI signals with wide diusion gradients is presented. The derived relations in the previous subsection for the inclusion of relaxation eects, will be incorporated in the diusion MRI equations to provide an inclusive signal generation scheme. In real diusion MRI data acquisition experiments, the diusion gradients take a large portion of imaging time. In other words, despite the previous formulation assumptions [16, 48] in Fig. 4.9 is not innitesimally small, and diusion time 2 + comprises most of the imaging period. For the rest of our derivations, the following practical assumptions are made. First, the spin diusive motions occur only during the diusion time introduced before and 111 the translational motions of particles during the application of x- and y- gradients are negligible. Second, since the magnitudes and application periods of diusion gradients for wide pulses are much larger than other existing gradients, a spin phase shift due to its diusive translation is predominantly caused by the diusion gradient. With the above assumptions, we can proceed with the modication of Eq. (4.20) with two major goals; namely, inclusion of the T 1 and the T 2 parameters, and factoring the existence of very wide diusion gradients. To begin with, it is necessary to modify theT 1 andT 2 values in Eqs. (4.21) and (4.22), respectively. Originally, these equations are the established parallel local magnetization component after time T R , and the decayed transverse component of local magnetization after time T E , respectively, both calculated at position ~ r. For a stationary spin magnet at~ r, it is imposed to the same rates of theT 1 and theT 2 relaxations for the whole period of imaging since those parameters are the material characteristics. However, for the case of diusion MRI, since spins are moving across dierent tissue materials from position ~ r to ~ r 0 , they experience dierent rates of relaxations. Under spin diusive translations, if the relaxations still have exponential regimes, it will be possible to nd analytical modications to Eqs. (4.21) and (4.22) as presented below. Fig. 4.11 illustrates a spin magnet starting from position ~ r 1 traveling to ~ r 2 during large time interval T . There are always N T small equal time steps, t, to partition the time interval into, such that the spin remains in one material during the entire time of 112 Figure 4.11: Traveling spin magnet: A spin magnet traveling across dierent materials (color coded) with dierent T 1 parameters. any step. The parallel local magnetization component in Eq. (4.4) for a small patch of spins traveling from ~ r 1 to ~ r 2 during time interval T R can be modied as M k (~ r 1 ;~ r 2 ;t =T R ) = M 0 0 @ 1 N T Y i=1 exp t=T 1;i 1 A ; (4.56) = M 0 1 exp t P N T i=1 (1=T 1;i ) ; (4.57) = M 0 1 exp T R 1 N T P Nc c=1 (nc=T 1;c ) ; (4.58) = M 0 1 exp T R P Nc c=1 (ft;c=T 1;c ) ; (4.59) = M 0 1 exp T R =T e 1 (~ r 1 ;~ r 2 ) ; (4.60) 113 wherec is the material index,f t;c is the time fraction the spin magnets travels in material c, T 1;c is the relaxation rate for material c, and the eective T 1 relaxation rate is dened as T e 1 (~ r 1 ;~ r 2 ) = Nc X c=1 (f t;c =T 1;c ): (4.61) As a result of all spins traveling to that position after time T R , the total value of parallel magnetization at ~ r 2 is derived by integration over the entire sample, M k (~ r 2 ) = Z ~ r 1 P ~ r 1 (~ r 2 ;T R )M k (~ r 1 ;~ r 2 ;t =T R )d~ r 1 ; (4.62) and the value of transverse magnetization right after the RF pulse at time t = 0 (next repetition cycle) is equivalent to Eq.(4.62). That is, we have M ? (~ r;t = 0) =M k (~ r): (4.63) With a similar procedure, the eective value for the T 2 relaxation parameter can be derived and the nal result is T e 2 (~ r 1 ;~ r 2 ) = Nc X c=1 (f t;c =T 2;c ); (4.64) 114 and, for a small spin patch moving from ~ r to ~ r 0 during the time interval T E , Eq. (4.22) can be modied by M T E (~ r; ~ r 0 ) = M ? (~ r;t = 0)P ~r ( ~ r 0 ;T D ) exp T E =T e 2 (~ r; ~ r 0 ) : (4.65) The derived eective relaxation coecients in Eqs. (4.61) and (4.64) might provide insights for understanding the relatively low value of relaxation coecients in white mat- ters. The white matter volume is mostly comprised of the cytoplasm and extracellular uids which exhibit large values ofT 1 andT 2 , and a small portion occupied by the myelin sheath which is mostly of fatty material with much smaller relaxation coecients com- pared to those of uids. The white matters, however, exhibit relaxation parameters which are close to those of the myelin, which is unexpected given the small volume fraction of those regions. Having a closer look at the eective relaxation parameter equations reveals that, for the white matters, since the larger a material compartment is, the faster the spins travel, the time fractions f t;c for all composing material are in the same range, and the eective parameters are dominated by the material with small relaxation parameters, which is the myelin. For a wide diusion gradient, the phase encoding occurs during the entire pulse in- terval, as opposed to the innitesimally narrow pulse where it happens instantaneously. In order to quantify the phase shifts, it is necessary to estimate the spin motion equation during the applied diusion gradient. In the following, we show that the the motion with a constant speed provides an accurate model for our purpose. 115 Figure 4.12: A spin traveling in dierent materials. Biological tissues often have structured and repetitive geometries. For example, in the brain white matters, the structure consists of parallel organization of dierent material with periodic geometry. Fig. 4.12 illustrates a simple 1D motion in an environment consisting of three dierent materials which repeat periodically for n times (the number three for the number of material is arbitrary and can be any number). The material compartments are assumed to be homogeneous and, as a result, a spin travels with constant speeds of v 1 , v 2 , and v 3 in each compartment. The compartment widths are denoted byd 1 ,d 2 andd 3 and, consequently, the times it takes the spin magnet to traverse the width of a compartment are t 1 = d 1 v 1 , t 2 = d 2 v 2 , and t 3 = d 3 v 3 . Also a block width is dened as D =d 1 +d 2 +d 3 and T =t 1 +t 2 +t 3 is the time the spin travels in a block. The average speed after a spin passes n blocks can be computed as 116 V = n 2 (d 1 +d 2 +d 3 ) n (t 1 +t 2 +t 3 ) ; = nD 1 : (4.66) Now, let us consider a case when a spin travels fort =n 0 T +t seconds, wheret<T . The exact value of the distance the spin travels is given by d actual =D +f(t); (4.67) where f(:) is a function specifying the exact position of the spin within the D + 1 block. On the other hand, assuming a constant speed motion for the spin with average speed given in Eq. (4.66), the distance traversed by the spin can be estimated by ^ d = t V; nD +t 0 D T : (4.68) There are two scenarios to evaluate. First, for large diusion time, the nD terms in the left-hand-sides of Eqs. (4.67) and (4.68) become dominant and the estimated distance converges toward the true value. In contrast, for small diusion time, there might be a discrepancy between those values (although still small), but that error becomes insignicant. Recall that for small diusion time, the spin traverses over a small distance and the corresponding phase shift is negligible. Thus, the constant speed model can be safely employed in the subsequent derivations. 117 For a spin moving from position~ r to ~ r 0 and experiencing G x , G y , and g gradients as shown in Fig. 4.9, the total phase shift (namely, the phase of the complex exponential in Eq. (4.24)) at time t near echo time T E can be expressed as total (t) =! r t + x (t) + y + g : (4.69) We start the derivation with the phase shift due to the diusion gradient. It can be written as a sum of two terms, g = g + + g ; (4.70) where the + and subscripts represent the shifts associated with the rst and second pulses of diusion gradients. We have ! g = ~ g + " ~ r + ~ r 0 ~ r 1 t # ; (4.71) T D = + 2; (4.72) g + = Z t 1 + t 1 ! g + (t)dt = ~ g + " ~ r + ( 2 + 2t 1 ) ~ r 0 ~ r 2 # ; (4.73) g = Z t 1 +2+ t 1 ++ ! g (t)dt = + ~ g + " ~ r + (2 2 + 2t 1 + 2) ~ r 0 ~ r 2 # ; (4.74) 118 and, consequently, we obtain g = ~ g + ( ~ r 0 ~ r): (4.75) = +2~ q m ( ~ r 0 ~ r); (4.76) ~ q m = ~ g + : (4.77) With similar procedures, the phase shift due to the application of two G x pulses can be summarized as x (t) = ~ G x " (T s =2t 0 )~ r + (T 2 s =4 +T x T s t 02 2t 0 T x t 0 T s ) ~ r 0 ~ r 2T D # ; (4.78) t 0 = tT E : (4.79) For long diusion time, we have T x T s and T D T X . As a result, Eq. (4.78) can be simplied as x (t) = (T s =2t 0 ) ~ G x ~ r 0 ; = 2k x x 0 : (4.80) Finally, for the G y gradient, we have y = ~ G y " y ~ r + ( 2 y 2 y T y ) ~ r 0 ~ r 2T D # ; (4.81) 119 which can be simplied as y = y ~ G y ~ r 0 = 2k y y 0 : (4.82) given that T y y and T y T D . Following a procedure similar to the MRI case, the k-space signal for the Diusion MRI, which includes the eects of relaxations and can accommodate wide diusion gra- dient pulses, can be summarized as S( ~ k;~ qm) = R ~ r !r ? exp i ~ rec M ? (~ r;t=0) exp i2 ~ k~ r R ~ R P ~ r (~ r+ ~ R;T D ) exp T E =T e 2 (~ r;~ r+ ~ R) exp 2kx GxT y 2 (~ r+ ~ R) exp i2( ~ qm ~ k) ~ R d ~ R d~ r; (4.83) In our derivation, it is assumed that the receiver coil magnetic eld is uniform over the entire sample. It can be safely assumed that ~ q m ~ k ~ q m . The same discussion in Sec. 4.2.1 for interpretation of the reconstructed MR image, applies to Eq. (4.83) as well. It is also easy to validate that, by ignoring the relaxation parameters and having very narrow diusion pulse, Eq. (4.83) reduces to Eq. (4.20). To see how dierent is Eq.(4.83) from that of thekspace signal under ideal imaging conditions in Eq.(4.20), we consider a simple 1D reconstruction scenario. Let us assume that the environment is homogeneous, which implies that the self-diusion propagator is governed by the Gaussian distribution. Also, the important identity in MRI is the end image not the kspace signal. For the ideal image reconstruction case in Eq.(4.20) 120 the image is simply obtained by applying the inverse transform to the kspace signal. Provided that the Fourier transform of the Gaussian function, is the Gaussian function, we have the following relation for the image, E(~ q) = M 0 ~ r! r ? e i e 2 2 2 q 2 ; (4.84) = p 2DT D : For practical image reconstruction situations, however, following the discussions on the interpretation of theMRI image, we arrive at the following derivation for the recon- structed image, E m (~ q) = M 0 ~ r! r ? e i (1e T R =T 1 )e T E =T 2 : 1 2K M C [T 2 (1e CK M T 2 ) T + 2 (1e CK M T + 2 )]e 2 2 2 q 2 m ; (4.85) where C = 2 Gx , and T 2 and T + 2 are the values of relaxation parameter before and after the echo time, respectively. We evaluate the relative error, RE = Em(~ q)E(~ q) Em(~ q) , under two extreme cases. First, assume that the relaxation parameters are very large. Fig.4.13 shows the variation of the relative error as a function of the pulse width duration. To have a fair comparison, the magnetic energy measure in two experiments is kept equal to g = 0:8. It can be seen that, for the narrow pulse width, the error is negligible, but as we increase the gradient duration, the discrepancy between the two derivations becomes 121 Figure 4.13: Relative error: Relative error as a function of diusion the pulse width. signicant. For the case study involves the condition where the diusion gradient is very narrow- = 1 [ms], the variable parameter is the T 2 parameter. Fig.4.14 illustrate the variation of the relative error, as a function of the T 2 parameter. Fig.4.14 shows that, in the range of typical relaxation parameter value of white matter, the error is signicantly high. 4.3 Conclusion and Major Contributions In this chapter, we generalized the MRI and diusion MRI signal generation formulations. The research presented in this chapter has several major contributions. Our rst major contribution is the inclusion of theT 1 and theT 2 relaxation parameters into the MRI signal generation formulation. The denition of the eective spin magnet density is modied under the relaxation inclusion. The new denition of the eective density is a function of the k-space variable. The perfect Fourier relations between the 122 Figure 4.14: Relative error as a function of T 2 . k-space signal and the eective density no longer holds and the application of the inverse transform toS( ~ k) will not result in the density. However, the common practice is to apply the inverse Fourier transform to thek-space signal for image reconstruction. Knowing that the reconstructed image diers from the eective spin density, we devised a formulation based on the properties of the Fourier transform and quantied the relation between them. This will be of special importance for our ecient computation method for S( ~ k) signal simulation in complex biological tissue geometries which is inevitable for the simulated reconstruction as opposed to real data acquisitions. In Section 4.1, it was stated that the formulation for Diusion MRI signal genera- tion was applicable if the diusion gradient was innitesimally narrow. In reality, this condition never holds and a large portion of the imaging time is taken by wide diusion gradients. For the narrow pulse scenario, the spin phase shifts occur instantaneously. On the other hand, for wide pulses, the phase encodings happen dynamically during the time 123 interval of wide pulse application. As a results, it is necessary to model the motions of the spin magnets. Our second major contribution is to present a novel yet simple method for the prediction of spin magnet motions. Based on our motion prediction scheme, an ana- lytical method was provided for the computation of the phase shifts due to the presence of wide diusion gradients. Finally, the modied formulation for generation of diusion MRI signal under wide diusion pulse assumption was provided. 124 Chapter 5 : Conclusion and Future Work 5.1 Conclusion In this work, we provided a methodology for realistic simulation of D-MRI signal in the brain white matter tissues. It accommodates the actual values of the parameters in real world MRI scanners. Since the diusive translations of particles are not in uenced by the presence of the environment magnetic eld, the self-diusion process occurs independently of the MRI parameters and can be studied separately. In developing the method, we took a two-stage approach { 1) quantication of self-diusion in the white matter media, and 2) developing inclusiveD-MRI theory that incorporates real imaging and tissue structural parameters. The rst part of this study focused on the quantication of the self-diusion process. The problem was handled by addressing the solution of the governing PDE. Due to the white matter complicated geometry, analytical solutions cannot be derived and numeri- cal methods were employed. Considering the irregular boundary conditions imposed by the white matter environment, the FEM provided the best tool in problem formulation 125 and solution. Due to the discontinuity of the self-diusion coecient on the material compartments of the white matter, the standard FEM cannot be properly applied in a straightforward manner. To resolve the issue mentioned above, we devised two theoret- ical rules (i.e., constraining the FEM mesh generation procedure and the basis function design) in the development of feasible FEM solutions. Several major contributions were made along this line and summarized below. First, to handle the practical issues of ini- tial solution accurate approximation and ecient allocation of computational resources, we proposed concepts of delayed initialization and the eective radius of propagation, respectively. Second, for the validation of the self-diusion solutions for the cases with no ground truth, a self-validation technique was proposed to ensure the accuracy of the results. Third, the relation between the white matter microscopic parameters and the self-diusion aggregate behavior within the including MRI voxel was analytically derived. Finally, the above mentioned concepts were implemented in a software tool to address the solution of self-diusion in the white matters. One of the characteristics of the developed tool is that it allows the modular study of the constituting material of the white matter and its eects on self-diusion. Consequently, abnormalities such as demyelination can be precisely studied and their eects on the self-diusion process can be evaluated. In the second part of this study, we considered the extension of the proposed MRI and D-MRI formulation to accommodate a larger range of parameter variations. Specically, the T 1 and T 2 relaxation parameters were included in the formulations. We examined the general condition where those parameters were spatially varying. It was also argued that, under the severe relaxation eects, the k-space signal was not a Fourier transform pair with the spin magnet density. Then, we provided a mathematical formulation so as 126 to interpret the reconstructed image based on the properties of the Fourier transform. Finally, we provided a formulation for the D-MRI signal where the data acquisition was performed with a wide diusion gradient. 5.2 Future Work The study conducted in this dissertation can be further extended in a number of dimen- sions in the future as detailed below. The ultimate goal is to leverage lessons learned from disintegrated D-MRI formation processes and simulations to devise better tractography schemes. The immediate research topic to pursue is the development of a D-MRI phantom for evaluation of tractography methods. As previously stated, in vivo validation of D-MRI tractography is a dicult task which remains an open issue. It is desired to be condent whether an identied ber tract by a method truly represents a white tract in the brain and how closely they are registered together. This issue becomes more sensitive since the precision of the brain mapping technique plays a substantial role in brain surgeries and even a few millimeters matter. It was shown that the D-MRI signal is characterized by the self-diusion propagators of voxels. Since our method provides a solution to the average propagator, it can be incorporated in a simulator framework for the generation of the gold standard for evaluation and comparison of tractography methods. This will be the center of our endeavor for future exploration. 127 Another application of our current results is the reverse engineering of neurode- generative diseases/abnormalities and quantitative study of resultant eects on the D-MRI signal. As mentioned before, D-MRI has found many applications in di- agnoses and studies of a number of conditions associated with the CNS. Many of those disorders occur at the microscopic level, which go beyond the scope ofD-MRI voxel-level measurements. For the diagnostic purpose, the trend is to scan patients with the abnormality condition and search for common signatures of the disease. The potential problem with this approach is that it is subject to errors and it lacks locality since the ground truth is not available on the ne resolution scale. On the other hand, it is possible to induce those microscopic abnormalities in our method and study their eects on simulated D-MRI signals. The advantage of this tech- nique is that it provides the ground truth at the voxel level which makes it possible to have more accurate studies. Since our approach represents the self-diusion propagator as a function of spatial and temporal variables, it should give invaluable insights to the design of more ecient parameters of the D-MRI scanner. For example, we showed that the self- diusion prole in the case of demyelinated axons is lateral in the beginning but becomes symmetric gradually as time elapses in Chapter 3. This information can be exploited in the design of pulse sequences of diusion gradients in data acquisition. 128 Bibliography [1] C. Adler, S. Holland, V. Schmithorst, and et al., \Abnormal frontal white matter tracts in bipolar disorder: a diusion tensor imaging study," Bipolar Disorders, pp. 197{203, 2004. [2] M. Aksoy, C. Liu, M. Moseley, and R. Bammer, \Single-step nonlinear diusion tensor estimation in the presence of microscopic and macroscopic motion," Magnetic Resonance in Medicine, pp. 1138{1150, 2008. [3] M. Ashtari, S. Kumra, S. Bhaskar, and et al., \Attention-decit/hyperactivity dis- order: a preliminary diusion tensor imaging study," Biological Psychiatry, pp. 448{ 455, 2005. [4] R. Bammer, M. Augustin, S. Strasser-Fuchs, and et al., \Magnetic resonance diu- sion tensor imaging for characterizing diuse and focal white matter abnormalities in multiple sclerosis," Magnetic Resonance in Medicine, pp. 5813{591, 2000. [5] N. Barnea-Goraly, H. Kwon, V. Menon, S. Eliez, L. Lotspeich, and A. Reiss, \White matter structure in autism: preliminary evidence from diusion tensor imaging," Biological Psychiatry, pp. 323{326, 2004. [6] Basser, Mattielo, Turne, and L. Bihan, \Mr diusion tensor spectroscopy and imag- ing," Biophy Journal, pp. 259{267, 1994. [7] Basser, Mattielo, Turner, and L. Bihan, \Diusion tensor echo-planner imaging of human brain," SMRM, p. 584, 1993. [8] P. Basser, S. Pajevic, C. Peirpaoli, and A. Aldroubi, \Fiber tract following in the human brain using dt-mri data," IEICE Trans. Inf. Syst, pp. 15{21, 2002. [9] P. Basser, S. Pajevic, C. Pierpaoli, J. Duda, and A. Aldroubi, \In vivo ber tractog- raphy using dt-mri data," Magnetic Resonance in Medicine, pp. 625{632, 2000. [10] P. Basser, Mattielo, Turner, and L. Bihan, \Estimation of the eective self-diusion tensor from the nmr spin echo," Journal of Magnetic Resonance, pp. 247{254, 1994. [11] A. Berezhkovskii and G. Weiss, \Diusion in multilayer media: Transient behavior of the lateral diusion coecient," Chemical Physics, p. 154710, 2006. [12] M. Bozzali, A. Falini, M. Franceschi, and et al., \White matter damage in alzheimer's disease assessed in vivo using diusion tensor magnetic resonance imaging," Journal of Neurol Neurosurgery Psychiatry, pp. 742{746, 2002. 129 [13] K. Brownstein and C. Tarr, \Importance of classical diusion in nmr studies of water in biological cells," Physical Review A, pp. 2446{2453, 1979. [14] J. C. Butcher, Numerical methods for Ordinary Dierential Equations. John Wiley and Sons, Inc., New Jersey, 2003. [15] M. Caan, A. de Vries, G. Khedoe, E. Akkerman, L. van Vliet, K. Grimbergen, and F. Vos, \Generating ber crossing phantoms out of experimental dwis," in Medical image computing and computer-assisted intervention, 169-176. [16] P. Callaghan, Principles of nuclear magnetic resonance microscopy. Oxford Uni- versity Press, 1994. [17] A. Castriota-Scanderbeg, F. Fasano, G. Hagberg, U. Nocentini, M. Filippi, and C. Caltagirone, \Coecient d(av) is more sensitive than fractional anisotropy in monitoring progression of irreversible tissue damage in focal nonactive multiple scle- rosis lesions," AJNR American Journal of Neuroradiology, pp. 663{670, 2003. [18] B. Chen and A. Song, \Diusion tensor imaging ber tracking with local tissue property sensitivity: phantom and in vivo validation," Magnetic Resonance Imaging, pp. 103{108, 2008. [19] M. Chen and D. Mogul, \A structurally-detailed nite element human head model for brain-electromagnetic eld simulations," in Proceeding of 3rd International IEEE/EMBS Conference on Neural Engineering, 2007. [20] P. Cheng, V. Magnotta, D. Wu, P. Nopoulos, D. Moser, J. Paulsen, R. Jorge, and N. Andreasen, \Evaluation of the gtract diusion tensor tractography algorithm: a validation and reliability study," Neuroimage, pp. 1075{1085, 2006. [21] M. Chung, \Diusion smoothing on brain surface via nite element method," in In Proceedings of IEEE International Symposium on Biomedical Imaging, 2004. [22] C. Clark, D. Werring, and D. Miller, \Diusion imaging of the spinal cord in vivo: estimation of the principal diusivities and application to multiple sclerosis," Mag. Res. In Medicine, pp. 133{138, 2000. [23] L. Demkowicz, W. Rachowicz, and P. Devloo, \A fully automated hp-adaptivity," in TICAM Report , The University of Texas at Austin, 01-28. [24] H. Edelsbrunner, Geometry and topology for mesh generation. Cambridge Mono- graph on Applied and Computational Mathematics, 2001. [25] S. EE and Y. Song, \Multiple echo diusion tensor acquisition technique," Magnetic Resonance Imaging, pp. 7{18, 2006. [26] M. Filippi, M. Cercignani, M. Inglese, M. Horseld, and G. Comi, \Diusion tensor magnetic resonance imaging in multiple sclerosis," Neurology, pp. 304{311, 2001. 130 [27] J. Ford and D. Hackney, \Numerical model for calculation of apparent diusion coecients (adc) in permeable cylinders- comparison with measured adc in spinal cord white matter," Magn Reson Med., pp. 387{394, 1997. [28] J. Ford, D. Hackney, E. Lavi, M. Phillips, and U. Patel, \Dependence of apparent diusion coecients on axonal spacing, membrane permeability,and diusion time in spinal cord white matter," Magn Reson Med, pp. 778{782, 1998. [29] C. Gao, F. Tay, and W. Nowinski, \A nite element method based deformable brain atlas suited for surgery simulation," in Proc IEEE Eng Med Biol Soc, 2005. [30] F. L. Gossl, C and, B. Putz, and D. Auer, \Fiber tracking from dti using linear sate space models: detectability of thr pyramidal tract," NeuroImage, pp. 378{388, 2002. [31] E. Haake, R. Brown, M. Thompson, and R. Venkatesan, Magnetic Resonance Imag- ing: Physical Principles and Sequence Design. Wiley-Liss, 1999. [32] R. Henry, J. Berman, S. Nagarajan, P. Mukherjee, and M. Berger, \Subcortical path- ways serving cortical language sites: initial experience with diusion tensor imaging ber tracking combined with intraoperative language mapping," Neuroimage, pp. 616{622, 2004. [33] R. Hiptmair, \Canonical construction of nite elements," Math. Comput, pp. 1325{ 1346, 1999. [34] S. Hwang, C. Chin, F. Wehrli, and D. Hachney, \An image-based nite dierence model for simulating restricted diusion," Mag. Res. In Medicine, pp. 373{382, 2003. [35] K. Kamada, T. Todo, A. Morita, and et al., \Functional monitoring for visual path- way using real-time visual evoked potentials and optic-radiation tractography," Neu- rosurgery, pp. 121{127, 2005. [36] M. Kubicki, C. Westin, P. Nestor, and et al., \Cingulate fasciculus integrity dis- ruption in schizophrenia: a magnetic resonance diusion tensor imaging study," Biological Psychiatry, pp. 1171{1180, 2003. [37] D. Le Bihan, E. Berton, D. Lallemand, P. Grenier, E. Cbanis, and M. Laval-Jeantet, \. mr imaging of intravoxel incoherent motions: application to diusion and perfusion in neurologic disorders," Radiology, pp. 401{407, 1986. [38] J. Lee, M. Han, S. Kim, O. Kwon, and J. Kim, \Fiber tracking by diusion tensor imaging in corticospinal tract stroke: topographical correlation with clinical symp- toms," Neuroimage, pp. 771{776, 2005. [39] T. Li, V. Mathews, Y. Wang, D. Dunn, and W. Kronenberger, \Adolescents with dis- ruptive behavior disorder investigated using an optimized mr diusion tensor imaging protocol," Annals NY Academy of Science, pp. 184{192, 2005. [40] K. Lim, M. Hedehus, M. Moseley, A. de Crespigny, E. Sullivan, and A. Pfeerbaum, \Compromised white matter tract integrity in schizophrenia inferred from diusion tensor imaging," Arch Gen Psychiatry, pp. 367{374, 1999. 131 [41] N. Lori, E. Akbudak, J. J Shimony, T. Cull, A. Snyder, R. Guillory, and T. Conturo, \Diusion tensor ber tracking of human brain connectivity: aquisition methods, reliability analysis and biological results," NMR in Biomedicine, pp. 493{515, 2002. [42] S. Lu, D. Ahn, G. Johnson, and S. Cha, \Peritumoral diusion tensor imaging of high-grade gliomas and metastatic brain tumors," AJNR American Journal of Neuroradiology, pp. 937{941, 2003. [43] Y. Lu, A. Aldroubi, J. Gore, A. Anderson, and Z. Ding, \Improved ber tractography with bayesian tensor regularization," Neuroimage, pp. 1061{1074, 2006. [44] D. Merhof, G. Soza, A. Stadlbauer, G. Greiner, and C. Nimsky, \Correction of susceptibility artifacts in diusion tensor data using non-linear registration," Medical image analysis, pp. 588{603, 2007. [45] S. Mori and P. van Zijl, \Fiber tracking principles and strategies," NMR in Biomedicine, pp. 468{480, 2002. [46] K. Mrboldt, W. Hnicke, and J. Frahm, \Self-diusion nmr imaging using stimulated echoes," Journal of Magnetic Resonance, pp. 479{486, 1985. [47] P. Nucifora, R. Vrma, S.-K. Lee, and E. Melhem, \Diusion-tensor mr imaging and tractography: Exploring brain microstructure and connectivity," Radiology, p. 245, 2007. [48] E. Ozarslan, C. Koay, and P. Basser, \Remarks on q-space mr propagator in par- tially restricted, axially symmetric, and isotropic environment," Magnetic Resonance Imaging, pp. 834{844, 2009. [49] E. Ozarslan and T. Mareci, \Generalized diusion tensor imaging and analytical relationships between diusion tensor imaging and high angular resolution diusion imaging," Magnetic Resonance in Medicine, pp. 955{965, 2003. [50] S. Pajevic, A. Aldroubi, and P. Basser, \A continuous tensor eld approximation of discrete dt-mri data for extracting microstructural and architectural features of tissue," Journal of Magnetic Resonance, pp. 85{100, 2002. [51] J. Powel, M. Mallett, G. Rickayzen, and W. Evans, \Exact analytic solutions for diusion impeded by an innite array of partially permeable barriers," in Proceedings of Mathenmatical and Physical Sciences, 1992. [52] G. Price, M. Bagary, M. Cercignani, D. Altmann, and M. Ron, \The corpus callosum in rst episode schizophrenia: a diusion tensor imaging study," Journal of Neurol Neurosurgery Psychiatry, pp. 585{587, 2005. [53] S. Reading, M. Yassa, A. Bakker, and et al., \Regional white matter change in pre-symptomatic huntington's disease: a diusion tensor imaging study," Psychiatry Res, pp. 55{62, 2005. 132 [54] S. Rose, F. Chen, J. Chalk, and et al., \Loss of connectivity in alzheimer's disease: an evaluation of white matter tract integrity with colour coded mr diusion tensor imaging," Journal of Neurol Neurosurgery Psychiatry, pp. 528{530, 2000. [55] L. S, D. Ahn, G. Johnson, M. Law, D. Zagzag, and R. Grossman, \Diusion-tensor mr imaging of intracranial neoplasia and associated peritumoral edema: introduction of the tumor inltration index," Radiology, pp. 221{228, 2004. [56] P. Sen and P. Basser, \A model for diusion in white matter in the brain," Biophy- isical Journal, pp. 2927{2938, 2005. [57] P. Solin, Partial dierential equations and the nite element method. John Wiley and Sons, Inc., Hoboken, New Jersey, 2006. [58] P. Solin, K. Segeth, and I. Dolezel, Higher-order nite element methods. Chapman and Hall/CRC, 2004. [59] G. Stanotsz, A. Szafer, G. Wright, and R. Henkelman, \An analytical model of restricted diusion in bovine optic nerve," Magn Reson Med, pp. 103{111, 97. [60] E. Stejskal and J. Tanner, \Spin diusion measurements: spin echoes in the presence of time-dependent eld gradient," Journal of Chemical Physiscs, pp. 288{292, 1965. [61] A. Sukstanskii and D. Yablonsky, \Eect of restricted diusion on mr signal forma- tion," Magn Reson, pp. 92{105, 2002. [62] A. Szafer, J. Zhong, and J. Gore, \Theoretical model for water diusion in tissues," Magn Reson Med, pp. 697{712, 1995. [63] P. Szeszko, B. Ardekani, M. Ashtari, and et al., \White matter abnormalities in obsessive-compulsive disorder: a diusion tensor imaging study," Arch Gen Psychi- atry, pp. 782{790, 2005. [64] K. Tabelow, J. Polzehl, V. Spokoiny, and H. Voss, \Diusion tensor imaging: struc- tural adaptive smoothing," Neuroimage, pp. 1763{1773, 2008. [65] J. Tanner, \Transient diusion in system partitioned by permeable barriers. ap- plication to nmr measurement with pulsed eld gradient," J. Chem. Physics, pp. 1748{1754, 1978. [66] D. Tayler and M. Bushnell, \The spatial mapping technique," Phys. Med. Biol, pp. 345{349, 1985. [67] C. Tench, P. Morgan, L. Blumhardt, and C. Constantinescu, \Improved white matter ber tracking using stochastic labeling," Magnetic Resonance in Medicine, pp. 677{ 683, 2002. [68] J. Tournier, F. Calamante, M. King, D. Gadian, and A. Connelly, \Limitations and requirements of diusion tensor ber tracking: An assessment using simulations," Magnetic Resonance in Medicine, pp. 701{708, 2002. 133 [69] A. Tropine, G. Vucurevic, P. Delani, and et al., \Contribution of diusion tensor imaging to delineation of gliomas and glioblastomas," Journal of Magnetic Resonance Imaging, pp. 905{912, 2004. [70] D. Tuch, \Diusion mri of complex tissue structure," in Ph.D. dissertation MIT, 2002. [71] D. Tuch, T. Reese, M. Wiegell, and V. Wedeen, \Diusion mri of complex neural architecture," Neuron, pp. 885{895, 2003. [72] S. Wakanna, H. Jiang, L. Nagae-Poetscher, P. van Zijl, and S. Mori, \Fiber tract- based atlas of human white matter anatomy," Radiology, pp. 77{87, 2004. [73] D. Werring, C. Clark, G. Barker, A. Thompson, and D. Miller, \Diusion tensor imaging of lesions and normal-appearing white matter in multiple sclerosis," Neu- rology, pp. 1626{1632, 1999. [74] D. Werring, A. Toosy, C. Clark, and et al., \Diusion tensor imaging can detect and quantify corticospinal tract degeneration after stroke," Journal Neurol Neurosurgery Psychiatry, pp. 268{272, 2000. [75] K. Yamada, H. Ito, H. Nakamura, O. Kizu, W. Akada, T. Kubota, M. Goto, J. Kon- ishi, K. Yoshikawa, K. Shiga, M. Nakagawa, S. Mori, and T. Nishimura, \Stroke patients' evolving symptoms assessed by tractography," Journal of Magnetic Reso- nance Imaging, pp. 923{929, 2004. [76] K. Yoshikawa, Y. Nakata, K. Yamada, and M. Nakagawa, \Early pathological changes in the parkinsonian brain demonstrated by diusion tensor mri," Journal of Neurol Neurosurgery Psychiatry, pp. 481{484, 2004. 134

Abstract (if available)

Abstract Diffusion MRI (D-MRI) has opened a new front for uncovering the convoluted structure of the central nervous system by providing capability for the non-invasive identification of white matter tract geometries in the brain. One of the major open issues in fully extending this technology to the clinical domain is the lack of in vivo validation of the results, which makes it diffcult to have objective comparisons of different D-MRI techniques. To this end, the application of simulated data from known ground truths appears to be the second best choice. It is well understood that, this imaging modality is characterized by the shape of the self-diffusion (SD) profile within the brain fibers. The previous methods for the quantification of the SD process in the white matter environments suffer from the lack of enough generality or solution precision, and often impose excessive computational complexity, which limits their full extent of applicability. The contributions of this research are two folds: 1) the development of a new numerical method to compute the self-diffusion SD profile and 2) the provision of a generic framework for reconstruction of diffusion MR images under imperfect imaging conditions.

Linked assets

University of Southern California Dissertations and Theses

Conceptually similar

PDF

Diffusion MRI white matter tractography: estimation of multiple fibers per voxel using independent component analysis

PDF

Diffusion MRI of the human brain: signal modeling and quantitative analysis

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

PDF

Novel theoretical characterization and optimization of experimental efficiency for diffusion MRI

PDF

Location of lesions on postmortem brain by co-registering corresponding MRI and postmortem slices

PDF

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

PDF

Signal processing methods for interaction analysis of functional brain imaging data

PDF

Representation problems in brain imaging

PDF

Pediatric magnetic resonance image processing: applications to posterior fossa cancer and normal development

PDF

Velocity-encoded magnetic resonance imaging: acquisition, reconstruction and applications

PDF

New theory and methods for accelerated MRI reconstruction

PDF

Visualizing and modeling vocal production dynamics

PDF

Characterization of the brain in early childhood

PDF

Statistical signal processing of magnetoencephalography data

PDF

Susceptibility-weighted MRI for the evaluation of brain oxygenation and brain iron in sickle cell disease

PDF

Applications of graph theory to brain connectivity analysis

PDF

Techniques for efficient cloud modeling, simulation and rendering

PDF

Signal processing methods for brain connectivity

PDF

Time-resolved laser induced fluorescence spectroscopy for detection of brain tumors

Asset Metadata

Creator Karimi-Ashtiani, Shahryar (author)

Core Title Theory and simulation of diffusion magnetic resonance imaging on brain's white matter

School Viterbi School of Engineering

Degree Doctor of Philosophy

Degree Program Electrical Engineering

Publication Date 08/09/2012

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

Tag brain imaging,MRI,OAI-PMH Harvest,white matter

Language English

Contributor Electronically uploaded by the author (provenance)

Advisor Kuo, C.-C. Jay (committee chair), Leahy, Richard M. (committee member), Singh, Manbir (committee member)

Creator Email karimias@usc.edu,shahryar.karimi@gmail.com

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

Unique identifier UC1478183

Identifier etd-KarimiAshtiani-3752 (filename),usctheses-m40 (legacy collection record id),usctheses-c127-379286 (legacy record id),usctheses-m3353 (legacy record id)

Legacy Identifier etd-KarimiAshtiani-3752.pdf

Dmrecord 379286

Document Type Dissertation

Rights Karimi-Ashtiani, Shahryar

Type texts

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

Repository Name Libraries, University of Southern California

Repository Location Los Angeles, California

Repository Email cisadmin@lib.usc.edu