Index of /rozprawy2/11225

Pełen tekst

(1)Faculty of Physics and Applied Computer Science. Doctoral thesis. Krzysztof Kłodowski. Development of the b-matrix spatial distribution diffusion tensor imaging with applications in porous media and soft tissue imaging. Supervisors: prof. dr hab. Henryk Figiel Faculty of Physics and Applied Computer Science, AGH, Cracow. dr Artur T. Krzyżak Faculty of Geology Geophysics and Environmental Protection, AGH, Cracow. Cracow, March 2017.

(2) Declaration of the author of this dissertation: Aware of legal responsibility for making untrue statements I hereby declare that I have written this dissertation myself and all the contents of the dissertation have been obtained by legal means.. data, podpis autora. Declaration of the thesis Supervisor: This dissertation is ready to be reviewed.. data, podpis promotora rozprawy. 2.

(3) I would like to cordially thank the following persons who supported research presented in this thesis: prof. Henryk Figiel for encouraging me to follow the route of scientific research in NMR, Artur T. Krzyżak for sharing with me his reach experience, ideas and thoughts and giving me the opportunity to carry out the research within his group, Paweł Banyś from John Paul II hospital in Cracow, who gratefully helped with data acquisition on Siemens scanner, Iwona Habina, Karol Borkowski and Grzegorz Machowski – colleges from my research group who took over lots of my responsibilities for the time of writing of this thesis. The authorities of the John Paul II hospital in Cracow for providing the clinical scanner for the purpose of diffusion tensor imaging. I would also like to thank for the financial support from: Marian Smoluchowski Scientific Consortium KNOW, National Centre of Research and Development within projects: 1. PBS2/A2/16/2013, 2. STRATEGMED2/265761/10/NCBR/2015.. 3.

(4) Acronyms BSD b-matrix spatial distribution. BSD-DTI b-matrix spatial distribution in diffusion tensor imaging. DAI diffusion anisotropic indices. DTI diffusion tensor imaging. DWI diffusion weighted imaging. EPI echo planar imaging. FA fractional anisotropy. FD fibre density. FLGS Gauss filter. FLMD median filter. FLMN mean filter. HARDI high angular resolution diffusion imaging. MRI magnetic resonance imaging. NMR nuclear magnetic resonance. PGSE pulsed gradient spin echo. PMMA polymethyl methacrylate. PVE partial volume effect. RF radio frequency. ROI region of interest. 4.

(5) RSD relative standard deviation. sBSD simplified b-matrix spatial distribution. SD standard deviation. SE-EPI spin echo echo planar imaging. SNR signal to noise ratio. sRA scaled relative anisotropy. uBSD uniform b-matrix spatial distribution.. 5.

(6) Contents Acronyms. 4. Streszczenie. 9. Abstract. 11. 1 Introduction. 13. 1.1. Aim and motivation. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 13. 1.2. Structure of this thesis . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 13. 1.3. Author’s contribution . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 14. 2 Theory. 15. 2.1. The theoretical description of an experiment . . . . . . . . . . . . . . . . . . . . 15. 2.2. Theoretical model of diffusion . . . . . . . . . . . . . . . . . . . . . . . . . . . . 16. 2.3. Nuclear magnetic resonance . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 17 2.3.1. 2.4. Quantum mechanical description of NMR . . . . . . . . . . . . . . . . . 18. Magnetic resonance imaging . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 23 2.4.1. Selective excitation and detection . . . . . . . . . . . . . . . . . . . . . . 23. 2.5. Diffusion tensor imaging . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 25. 2.6. b-matrix spatial distribution . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 27 2.6.1. Normalization of the b effective . . . . . . . . . . . . . . . . . . . . . . . 27. 2.6.2. Calibration of the b-matrix. 2.6.3. Derivation of the calibrated tensor . . . . . . . . . . . . . . . . . . . . . 29. 2.6.4. Convenient simplifications . . . . . . . . . . . . . . . . . . . . . . . . . . 29. . . . . . . . . . . . . . . . . . . . . . . . . . 29. 3 Experiments 3.1. Phantoms . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 32 3.1.1. 3.2. 32. BSD-DTI phantoms . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 33. BSD-DTI experiments . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 38 3.2.1. Bruker 9.4 T . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 38. 3.2.2. GE 3.0 T . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 39. 3.2.3. Siemens 3.0 T . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 41 6.

(7) 3.3. Data collections, datasets, experiments . . . . . . . . . . . . . . . . . . . . . . . 42. 4 Calculations. 46. 4.1. Software . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 46. 4.2. Standard DTI . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 47. 4.3. BSD-DTI . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 47. 4.4. 4.3.1. Spatial b effective . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 47. 4.3.2. Spatial b-matrix. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 48. Filters . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 48 4.4.1. Gaussian blur . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 49. 4.4.2. Median and mean filters . . . . . . . . . . . . . . . . . . . . . . . . . . . 49. 4.5. Bivariate polynomial fits . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 49. 4.6. b-matrix shimming . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 49. 4.7. Tensor analysis . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 50. 5 Results. 51. 5.1. Preselection of the results . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 51. 5.2. Visualization of the b-matrices . . . . . . . . . . . . . . . . . . . . . . . . . . . . 52. 5.3. Quantitative results . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 52. 5.4. Bruker . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 54. 5.5. 5.6. 5.7. 5.4.1. Isotropic phantom . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 54. 5.4.2. Anisotropic phantoms . . . . . . . . . . . . . . . . . . . . . . . . . . . . 57. GE . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 60 5.5.1. Isotropic phantom . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 60. 5.5.2. Anisotropic phantoms . . . . . . . . . . . . . . . . . . . . . . . . . . . . 62. Siemens . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 68 5.6.1. Isotropic phantom . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 68. 5.6.2. Anisotropic phantom . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 71. Soft tissue imaging . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 74. 6 Discussion and conclusions. 81. 6.1. Performance . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 81. 6.2. Summary . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 83. 6.3. Perspectives . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 83. Bibliography. 84. A Bruker data tables. 92. B GE data tables. 125. 7.

(8) C Siemens data tables. 179. 8.

(9) Streszczenie Zaproponowane niedawno podejście do obrazowania tensora dyfuzji z wykorzystaniem przestrzennego rozkładu macierzy b jest sposobem kalibracji, który umożliwia wyeliminowanie błędów systematycznych pojawiających się w macierzach b wyliczanych z wykorzystaniem powszechnego założenia o ich przestrzennej jednorodności. Kalibracja taka wymaga zastosowania dobrej jakości fantomów anizotropowych oraz ich precyzyjnego pozycjonowania w czasie akwizycji obrazów kalibrujących. Kompletna kalibracja wymaga wykonania obrazowania dyfuzyjnego fantomu izotropowego oraz sześciu położeń fantomu anizotropowego, co wymaga siedem razy więcej czasu niż zajmuje akwizycja obrazu badanej próbki. Ponadto, taka kalibracja pomimo, że bardzo dokładna, jest ograniczona tylko do przestrzeni obejmowanej przez fantom. Niniejsza rozprawa ma na celu pokonanie powyższych niedogodności i ograniczeń.. W. tym celu wykorzystane zostały dwa uproszczenia; jedno wynikające z założenia jednorodnego przestrzennie ważenia dyfuzyjnego oraz drugie zakładające jednorodność fantomu anizotropowego. Uproszczenia te, zredukowały liczbę niezbędnych do kalibracji położeń fantomu anizotropowego do trzech, jak również uprościły obliczenia. Ponadto, dokonana została próba zastosowania omawianej kalibracji po raz pierwszy do danych uzyskanych poprzez obrazowanie tensora dyfuzji in vivo dla ludzkiego mózgu. Otrzymane w wyniku kalibracji przestrzenne rozkłady macierzy b zostały poddane dalszym modyfikacjom z wykorzystaniem filtrów wygładzających, dopasowania dwuwymiarowych wielomianów oraz numerycznej procedury iteracyjnego dopasowania (shimmowania) macierzy b. Wprowadzone modyfikacje miały na celu wyznaczenie możliwie najbliższego rzeczywistemu rozkładu macierzy b oraz ekstrapolację wyników poza obszar fantomu wykorzystywanego do kalibracji. Wymienione modyfikacje zostały również zastosowane do uproszczonej kalibracji wykonanej jedynie w oparciu o fantom izotropowy. Dane wykorzystane do analizy zaproponowanych rozwiązań zostały uzyskane przy pomocy systemu badawczego (9.4 T) oraz dwóch skanerów klinicznych (3.0 T). Poprawa precyzji wyznaczenia tensora dyfuzji, mierzona jako obniżenie odchylenia standardowego średniej wartości własnych tensora dla fantomu izotropowego, zawiera się pomiędzy 20 % a 70 %, w zależności od analizowanego zestawu danych i wybranego obszaru. Poprawa precyzji sięgająca 62 % została zaobserwowana dla obszarów większych od fantomów wykorzystywanych do kalibracji, potwierdzając możliwość ekstrapolacji.. 9.

(10) Błędy systematyczne zostały zwizualizowane przy pomocy map macierzy b. Wykazały one największe zniekształcenia w przypadku użycia silnych gradientów dyfuzyjnych skierowanych wzdłuż osi głównych układu laboratoryjnego. Kalibracja dla danych pochodzących z obrazowania tkanek miękkich dała w efekcie obrazy wolne od artefaktów oraz wykazała rozsądne wartości anizotropii frakcyjnej zarówno wewnątrz obszaru fantomu, jak i poza nim. Zaproponowane modyfikacje mogą zostać wykorzystane również razem z kompletną kalibracją przestrzennego rozkładu macierzy b dając najwyższą możliwą poprawę precyzji, bądź też zostać użyte razem z uproszczonym podejściem, które nadal daje znaczną poprawę, ale przy ograniczeniu czasu kalibracji i uproszczeniu procedury.. 10.

(11) Abstract The recently proposed b-matrix spatial distribution approach to diffusion tensor imaging is a comprehensive calibration technique able to reduce systematic errors appearing in the bmatrices derived under common assumption of their spatial uniformity. The calibration requires good quality anisotropic phantoms and its precise positioning during the acquisition of calibration images. Complete calibration requires obtaining data of an isotropic phantom and six various positions of an anisotropic phantom, what gives seven times longer acquisition time than takes imaging of the analysed sample. Moreover the calibration, despite being very precise, is limited only to the space covered by the phantoms. This thesis aimed to overcome these inconveniences and limitations. In order to achieve that goal two simplifications have been introduced; the assumptions of spatially uniform diffusion weighting and the homogeneity of the anisotropic phantom. The simplifications reduced number of required anisotropic phantom positions to three and made the calculations easier. Moreover, the first attempt of calibration applied to the human brain data obtained by means of in vivo diffusion tensor imaging has been made. Spatially distributed b-matrices derived from the calibration were further modified with usage of smoothing filters, bivariate polynomial fitting and numerical b-matrix shimming procedure. The modifications aimed to derive the most accurate b-matrix distribution and to extrapolate the results beyond the region covered by the calibration phantom. The same procedures were also applied to the simplistic calibration with the isotropic phantom. Data used for the analysis were obtained on a research 9.4 T and two clinical 3.0 T scanners. The improvement of accuracy of the diffusion tensor imaging measured in terms of reduction of standard deviation of mean eigenvalues for the isotropic phantom ranged from 20 % to 70 % depending on the dataset and selected region of interest. Up to 62 % better accuracy was observed for the regions bigger than the calibration phantom, validating the extrapolation approach as well. The systematic errors were characterized with b-matrix maps and showed the most severe distortions in case of strong diffusion gradients in directions along the main axes of the laboratory reference frame. The calibration of the soft tissue images resulted with images devoid of artefacts and yielding reasonable quantitative values of fractional anisotropy in regions within and outside the. 11.

(12) calibration phantom. The developed modifications can be also applied to the complete b-matrix spatial distribution calibration giving the best possible accuracy, or can be used together with the simplified approach resulting with still significant improvement, but being easier to carry out and less time consuming.. 12.

(13) Chapter 1 Introduction Ten, komu wiadomo, iż nie należy jeść ryby nożem, może jeść rybę nożem.1 W. Gombrowicz. 1.1. Aim and motivation. The most general aim of this work was to develop b-matrix spatial distribution in diffusion tensor imaging (BSD-DTI) technique. Treated as a starting point this general aim branched into several more specific goals. These include experimental verification of the novel phantoms, experiments on clinical scanners, software development, further improvement of the calibration through application of b-matrix filters, b-matrix shimming and b-matrix polynomial fits. One of the serious limitations of the BSD-DTI is restriction of the calibrated region to the area covered by the calibration phantom. The proposed modifications aimed also to make possible reasonable extrapolation of the calibration. The final goal was application of the developed approaches to porous media and soft tissue imaging. Capillary anisotropic phantoms were treated as model porous media in the first case, whilst human brain in vivo diffusion tensor imaging (DTI) served as an example of the latter. The BSD-DTI is a novel technique with great potential for improvement of all kinds of nuclear magnetic resonance (NMR) diffusion measurements. I believe the simplified form of the BSD-DTI presented in this thesis together with the proposed modifications can become a robust and reliable tool for precise quantitative diffusion measurements in variety of applications.. 1.2. Structure of this thesis. The thesis starts with a theoretical description of diffusion models, NMR phenomenon, magnetic resonance imaging (MRI) and DTI followed by the b-matrix spatial distribution (BSD) 1. The one who knows he shouldn’t eat fish with a knife, can eat fish with a knife.. 13.

(14) approach to DTI and two simplifications of this technique making the measurements and tensor calculations less complicated, however, at the cost of accuracy. Chapter 3 starts with a description of three generations of the phantoms used during measurements carried out on a 9.4 T research scanner and two 3.0 T clinical tomographs. Hierarchy of the analysed data, and the detailed acquisition parameters are also presented. The following chapter 4 describes the algorithms of b-matrix spatial distribution derivation and introduces several modifications of the BSD proposed to achieve better improvement of the results through BSD-DTI calibration. Selected results and their analysis are presented in chapters 5 and 6. For the complete sets of results the reader is referred to the appendices.. 1.3. Author’s contribution. Unless one is dealing with pure mathematics or some aspects of theoretical physics, the scientific work is not lying in the domain of individualism any more. It is barely possible to be a sole author of a good experimental research. I do not treat this fact as a drawback, because actually it gives much more possibilities and makes the scientific research more sound in general. But in accordance to that I believe a clean distinction of what was exactly done by the author of the thesis is utmost appropriate. Therefore I greatly acknowledge the following contributions to the presented thesis: 1. The general concept of BSD-DTI calibration technique as well as design and idea of the anisotropic phantoms described in this thesis are by doctor Artur T. Krzyżak. 2. The phantoms were manufactured by professor Zbigniew Raszewski’s group at the Faculty of Advanced Technologies and Chemistry at Military University of Technology in Warsaw. 3. The data obtained on Bruker and GE scanners were provided by doctor Artur T. Krzyżak. 4. Karol Borkowski participated in software development and greatly contributed to the part responsible for visualization of the results.. 14.

(15) Chapter 2 Theory Was sich überhaupt sagen lässt, lässt sich klar sagen; und wovon man nicht reden kann, darüber muss man schweigen.1 L. Wittgenstein. 2.1. The theoretical description of an experiment. In general, description of the theoretical background of an experimental work is balancing between two extreme approaches. The former I would call a laconic experimentalist approach and it can be recognized from the two with the extremely concise form. The description is in general true, but nothing more than necessary is written. Usually such a description is readable only for those involved in the particular narrow field of research, and even they have to presume most of the content. The latter approach, which I call a garrulous experimentalist approach, is also true in general, but such a description flounders in a barren struggle of explaining way too much details than necessary. Proper balance between the two is the key to good theoretical description, sadly it is not easy to find. Good practise is to start somewhere in the middle and add a tendency towards one of the limiting approaches judging by the target reader. Living in hope this thesis may be useful for anyone interested in DTI I assumed the reader has a non-negligible background in mathematics, physics and some familiarity with NMR. I wish the following theoretical introduction avoided hidden (and often false) pre-assumptions commonly repeated in number of NMR textbooks. The common practise of presenting approximate description at first, and subsequent addition of consecutive corrections bringing closer to the truth was purposely abandoned in this thesis in favour of the completely inverse approach. The reason for such description is that starting with the possibly most precise picture constituents a core which is always valid. This 1. What can be said at all can be said clearly, and what we cannot talk about we must pass over in silence.. 15.

(16) core can be looked at through particular filters (approximations) giving valid picture for particular aspects. I believe such strategy roots the core as a natural reference when thinking of a problem, and all the simplistic views have to be derived from that core, thus it requires some effort thanks to which the reader should be much more aware of the simplifying assumptions he is taking.. 2.2. Theoretical model of diffusion. Driven by the lack of general law governing diffusion in a unit volume, and inspired by Fourier’s mathematical description of heat transport, Adolph Fick endeavoured to fill this gap. As an analogy to Fourier’s law he assumed the transfer of salt in water must be directly proportional to the gradient of concentration and inversely proportional to the distance [1]. From those assumptions bolstered by simple experiments he derived his famous law binding the diffusive flux with the concentration under the assumption of the steady state: J~ = −D∇ϕ,. (2.1). where: • J~ is the diffusion flux, • D is the diffusion coefficient, • ϕ is the concentration. For a diffusion process in which particles neither can be created nor destroyed the continuity equation is valid: ∂ϕ + ∇ · J~ = 0, (2.2) ∂t which combined with the Fick’s first law and under assumption the diffusion coefficient is constant in space yields the diffusion equation: ∂ϕ = D∇2 ϕ. (2.3) ∂t The here above equation even though still useful today is somewhat problematic in case of diffusion due to the Brownian motions in homogeneous fluid. If there is no concentration gradient there should be no diffusion in accordance to the Fick’s law. The problem can be worked around by addition of contrast agents of given concentration to the medium under question. Such approach is, however, undesirable or often not even possible in most of the biomedical experiments. The Brownian motions inspired Albert Einstein to take a look at the problem from microscopic perspective. Basing on the kinetic theory of heat Einstein considered probability of 16.

(17) finding a particle at a given position in space [2]. Interestingly such approach lead him to the classical diffusion equation, with the difference that instead of concentration the probability density appeared in the equation. Although mathematically the final results are equivalent, conceptually the microscopic approach gives better insight into the physical process of diffusion. One of the interesting solutions of the diffusion equation is a free diffusion case, which will be referred to throughout this work. Assuming the solution is separable, i.e. can be expressed as a product of separate functions for each of the space dimensions and time: ϕ(~r, t) = X(x) · Y (y) · Z(z) · T (t),. (2.4). and under initial condition the probability density at time zero is described with Dirac’s delta, after somewhat laborious math we obtain the solution: −(~ r −~ r 0 )2 1 4Dt · e , (2.5) (4πDt)n/2 where n is the number of spatial dimensions. Due to the Gaussian form of the above solution,. ϕ(~r − ~r 0 , t) =. the free diffusion is often referred to as Gaussian diffusion. From the above equation it follows that the root mean square motion is expressed with: q. h¯ r 2 i − h¯ r i2 =. 2.3. √. 2nDt.. (2.6). Nuclear magnetic resonance. Since the very first nuclear magnetic resonance curve was recorded by Rabi et al. in 1938, announcing a novel method of measuring nuclear magnetic moment [3], the phenomenon found a panoply of applications including spectroscopy, relaxometry, diffusometry and imaging. Even though the history of NMR is long and rich its theoretical description still meets with some inconsistencies. The problem stems from the fact that nuclear spin, being a prerequisite to observe the magnetic resonance, is purely quantum mechanical property, but at the same time the resonance itself can be successfully described by means of classical mechanics. The two descriptions are often mixed together resulting in a model being rather hard to comprehend. The source of this diversity can be found in two pioneering works on NMR observed in bulk matter by groups from Harvard (under the leadership of Purcell) and Stanford (directed by Bloch) [4, 5]. The papers were published almost at the same time and described the same phenomenon, however, seen from different perspectives. Referring to the Rabi’s pioneering paper the Harvard group wrote about transitions between energy levels split due to the Zeeman effect, whereas Bloch himself has seen it rather as a change in magnetic moment orientation on a resonance condition [6]. The two visions were so divert that Purcell and Bloch had to talk over half an hour before they realized they were talking about the same phenomenon [7]. The former approach is the hallmark of quantum mechanical description of NMR which is more 17.

(18) accurate, but sometimes counter intuitive, whilst the latter relates to classical description which is often sufficient - especially in case of imaging.. 2.3.1. Quantum mechanical description of NMR. Nuclei are characterized, among the others, with a quantum number called spin. Spin can take half integer or integer values including zero and is directly proportional to the magnetic moment of nuclei: ˆ µ ˆ = γ¯hI,. (2.7). where: µ ˆ is magnetic moment, γ is gyromagnetic ratio, h ¯ is reduced Planck constant, Iˆ is nuclear spin, and hats denote operators. Thus, nuclei with spin zero do not posses magnetic moments what makes them invisible for NMR experiments. In order to observe NMR phenomenon the sample has to be immersed in external homogeneous magnetic field B0 . In such conditions the Hamiltonian for a non zero spin nucleus can be written in a general form: ˆ = Hˆm + Hê , H. (2.8). where: Hê is an electric spin Hamiltonian, i.e. part representing interactions of the nucleus with electric field, and Hˆm is a magnetic spin Hamiltonian representing interactions with magnetic field. The former can be expanded to terms describing subsequent orders of electric multipoles. The number of terms in the expansion depends on the value of nuclear spin, and in the case of spin. 1 2. only the first term representing spherical charge distribution remains. It means. that electric energy is not affected by the rotations of the nucleus. Since hydrogen nucleus, i.e. protons being the only NMR active nuclei considered throughout this work have spin. 1 2. I. will limit further description to this particular case. Thus the remaining part of the general Hamiltonian is just the magnetic spin Hamiltonian: ~ Hˆm = −ˆ µ · B.. (2.9). Through substitution from equation 2.7 and the Larmore resonance condition: ωL = −γB,. (2.10). ˆ Hˆm = h ¯ ωL I,. (2.11). we obtain: where Iˆ is an adequate angular momentum operator:. 18.

(19) . . 1 0 1 Iˆx =  2 1 0 . 0. 1. (2.12a). . 1  Iˆy =  2i −1 0 . (2.12b). . 1 1 0 Iˆz =  2 0 −1. (2.12c). For a static magnetic field along the z axis the single spin Hamiltonian2 is: ˆ = ωL Iˆz . H. (2.13). In general a finite basis for the spin operator is formed by 2I + 1 Zeeman eigenstates. Thus in our case we have only two eigenstates, let denote them |αi and |βi. The eigen equations for those states are the following:. ˆ |αi = + 1 ωL |αi H 2 (2.14) ˆ |βi = − 1 ωL |βi H 2 with two eigenvalues: ± 21 ωL . The eigenstates form a convenient basis, however, are of little. importance in NMR since are stationary, i.e. do not change over time. Moreover, even though some of the spins3 can be found in eigenstates, much more probable for a spin is to occupy a superposition state: |ψi = cα |αi + cβ |βi ,. (2.15). where: cα and cβ are complex superposition coefficients and they are normalized: |cα |2 + |cβ |2 = 1.. (2.16). Time evolution In order to characterize time evolution of a single spin one has to solve time dependent Schrödinger equation:. d ˆ |ψt i . |ψt i = −iH dt After substitution from equation 2.13 and integration we obtain the solution: 0 |ψt i = exp −iωL τ Iˆz |ψt i ,. . . (2.17). (2.18). which can be shorten with usage of rotation operator Rˆz (φ) : 2 3. For simplicity the Hamiltonian was written in natural units, i.e. divided by ¯h. I take the liberty to use word spin interchangeably with phrase ’nuclei with spin. accepted jargon.. 19. 1 2’. as it is commonly.

(20) . . 1 exp −i 12 φ Rˆz (φ) = exp −iφIˆz = √  2 0 . . . 0. . (2.19).  ,. . exp +i 12 φ. to: 0 |ψt i = Rˆz (ωL τ ) |ψt i .. (2.20). Now, consider time evolution of two separate spins; one in the Zeeman eigenstate |ψ1 i = |αi, and another in a superposition state |ψ2 i = |+xi = √12 |αi + √12 |βi. The latter is denoted +x because it is an eigenstate of Iˆx operator. Assume time interval is equal to 2ωπL , so ωL τ = π2 . In the first case we get:. π |ψ1 i = Rˆz |αi , 2 what after writing in explicit matrix representation transforms to: . 0. . . 1 exp −i π4 |ψ1 i = √ 2 0. .  . 0. 0. 1. . exp +i π4. (2.21). π 1 = √ exp −i |αi . 4 2 0 . .   . (2.22). The latter case differs a little: . . 1 exp −i π4 |ψ2 i = 2 0 0. . 0.  . .  . . . 0 1 1 1 − i 1 1−i   =     =  = |+yi , π 2 2 1+i 0 1+i 1 exp +i 4 1 1. (2.23). where |+yi is an eigenstate of Iˆy operator. The first result proves what was stated above: the Zeeman eigenstates are stationary. After given period of time the spin is still in |αi state multiplied by a phase factor, which is not important in this case. However, the superposition state |+xi evolved into state |+yi. The time evolution of the superposition state resembles classical precession movement. The spin polarization changes its orientation over time as if was precessing about the magnetic field direction with Larmore frequency ωL . Spin ensemble The above description considers behaviour of a single spin. 1 2. in conditions of NMR experiment.. However in such a case more interesting is behaviour of the whole spin ensemble, which actually generates a measurable signal. The expectation value of an operator acting on a general superposition state: . cα. . |ψi =   , cβ. (2.24). ˆ = hψ|Q|ψi ˆ hQi .. (2.25). can be written as:. Using the operator defined as below: 20.

(21) . cα c∗α cα c∗β. |ψi hψ| = . cβ c∗α cβ c∗β. . (2.26). ,. the expectation value can be also expressed as: ˆ = T r{|ψi hψ| Q}. ˆ hQi. (2.27). For a spin ensemble consisting of N members the spin density operator can be defined as: . . . . N ρ ραβ cα c∗ cα c∗β 1 X ,  =  αα ρˆ = |ψn i hψn | =  α N n=1 ρβα ρββ cβ c∗α cβ c∗β. (2.28). and thus the expectation value for the whole ensemble becomes: ˆ hQi = T r{ˆ ρQ}.. (2.29). The diagonal elements of the density matrix represent the average populations of energy levels. The difference between populations of the two levels describes spin polarization along the external magnetic field. The off-diagonal elements represent quantum coherences, i.e. spin polarization perpendicular to the external field. Thermal equilibrium In thermal equilibrium the populations are governed by Boltzman distribution which incorporates the ratio of magnetic energy to thermal energy: exp (−¯ hωi /kB T ) ρii = P , hωj /kB T ) j exp (−¯. (2.30). where kB is Boltzman constant and T is temperature. Defining for simplicity the factor Bf : h ¯ γB0 , kB T and substituting ωα and ωβ from equation 2.14, we obtain: Bf =. ραα =. exp . . . . exp − 12 Bf + exp . ρββ =. 1 B 2 f. exp − 12 Bf . . (2.31). . 1 B 2 f. . . exp − 21 Bf + exp. . 1 B 2 f. (2.32). .. In usual conditions of imaging experiment there is huge excess of thermal energy over magnetic energy. For example in 3 T field in a room temperature the magnetic energy of proton ensemble is equal to 8.5 × 10−26 J while thermal energy is equal to 4.1 × 10−21 J. Thus the Bf factor is very small and the exponentials in equation 2.32 can be expanded to just first two terms of Taylor series. Then: 21.

(22) ραα =. 1 + 21 Bf 1 1 = + Bf , 1 1 2 4 1 − 2 Bf + 1 + 2 Bf. (2.33). and analogically:. 1 1 − Bf . (2.34) 2 4 On the other hand the coherences in thermal equilibrium are equal to zero. It means the ρββ =. spin polarizations are distributed isotropically in plane perpendicular to the magnetic field. Thus the explicit form of density matrix in thermal equilibrium is the following: . ρêq =. 1 2. + 14 Bf 0. 0 1 2. − 14 Bf.  .. (2.35). Sample excitation The density matrix in thermal equilibrium can be visualized as a net magnetization vector pointing along the external magnetic field. If the sample is placed in a radio frequency (RF) coil perpendicular to the external magnetic field and an alternate current of resonance frequency flows through the coil the difference of populations becomes smaller. The loss of the energy stored in this difference is transferred to coherences which represent spin polarization in the transverse plane. If the RF pulse of given amplitude is long enough the populations will equalize and the initial polarization is fully transferred to the transverse plane. Such RF pulse is called π 2. pulse because it can be visualized as rotation of the net magnetization vector about 90◦ . If. the RF pulse is about two times longer the inversion of the populations occurs, i.e. the diagonal elements of the density matrix replace each other, what can be visualized as a rotation of the net magnetization vector about 180◦ and thus the name π pulse. Such geometrical representation turns out to be very useful tool to describe the behaviour of an ensemble, not only in NMR, but in general case of two quantum level systems [8]. Mathematically the shape of density matrix after particular RF pulse can be calculated through application of a rotation operator to the density matrix in thermal equilibrium4 in order to obtain the results discussed above. The time evolution of the ensemble is similar to the evolution of a single spin. The net magnetization vector when flipped from the thermal equilibrium state is precessing about the direction of external magnetic field. The equation of motion derived from Schrödinger equation5 is equivalent to the classical equation of precession [10]: ~ dM ~ × B, ~ = γM dt 4. (2.36). Precisely it is much more convenient to transform the density matrix to rotating frame of reference and. then apply a rotation operator. 5 More precisely from Liouville equation which is a Schrödinger equation equivalent for density matrix [9].. 22.

(23) ~ is magnetization vector. The above result is a marvellous example of the correwhere M spondence principle. The quantum mechanical description of behaviour of macroscopic net magnetization is in perfect agreement with the classical description. From now on all further description will be mostly presented from this classical perspective. Relaxation The transverse component of precessing magnetization vector induces current flow in the receiver coil. Accordingly to the equation 2.36 the detected signal should be a simple sinusoidal wave. However, flipped magnetization can be expected to return to the thermal equilibrium state. In case of fast molecular motion (like in the case of liquids) the relaxation process can be described by phenomenological Bloch equations [10]: dMx ~ × B) ~ x − Mx = γ(M dt T2 dMy ~ × B) ~ y − My (2.37) = γ(M dt T2 dMz ~ × B) ~ z − Mz − M0 , = γ(M dt T1 where T2 and T1 are transverse and longitudinal relaxation time constants, respectively. The former defines how fast the magnetization disappears in the transverse plane, whilst the latter describes how fast the magnetization along the external field is restored. The transverse relaxation depends on direct spin-spin interactions, i.e. process of energy exchange between spins until they reach thermal equilibrium. The longitudinal relaxation depends on the spin-lattice relaxation, i.e. energy exchange between spins and the surrounding thermal reservoir. Usually T2 is visibly shorter than T1 .. 2.4. Magnetic resonance imaging. The above theoretical description is sufficient to design simple pulse experiments, predict time evolution of single spins as well as spin ensemble including relaxation process. However, in such perspective both excitation and detection of the signal pertain to the whole sample. In order to divide the analysed sample into small voxels and ascribe adequate portions of the measured signal to particular position in space one needs both selectively excite the sample and detect the signal.. 2.4.1. Selective excitation and detection. The resonance condition (equation 2.10) binds the resonance frequency with the magnetic field experienced by the local magnetization. Thus application of a magnetic field gradient in a given direction change the resonance condition along this direction: 23.

(24) ~ · ~r), ω(~r) = −γ(B0 + G. (2.38). ~ is magnetic field gradient. Combined with the selective RF pulse (soft pulse) exiting where G only narrow range of frequencies one can flip the magnetization in a selected thin slice of the sample. The stronger the applied gradient the thinner the excited slice can be. Set of three gradient coils generating magnetic field gradient in three orthogonal directions suffice to select slice in any direction. After excitation of a chosen slice the evolution of the magnetization in transverse plane can be derived by solving Bloch equations (2.37): Mxy (t) = (M0 cos(ωt)ˆ x + M0 sin(ωt)ˆ y ) exp (−t/T2 ) ,. (2.39). what can be shorten using complex number notation: Mxy (t) = M0 exp(iωt)exp (−t/T2 ) .. (2.40). Consider a signal measured from volume element dV . The maximum of the signal would be defined by the local spin density ρ(~r): ~ t) = ρ(~r)dV exp (iω(~r)t) exp (−t/T2 ) . dS(G,. (2.41). Inserting resonance frequency for the field modified by the gradient and neglecting the last term we get: . . ~ t) = ρ(~r)dV exp −iγ G ~ · ~rt . dS(G,. (2.42). The omission of the relaxation term is not harmful since the active gradient in usual conditions destroys signal much faster. In the above equation the γB0 term was also omitted since it can be thought of as a representation of the carrier frequency and thus the actually measured signal ~ · ~r. The signal integrated over entire volume would be: oscillates just at γ G S(t) =. ZZZ. . . ~ · ~rt d~r. ρ(~r)exp −iγ G. (2.43). For simplification we define the reciprocal space vector ~k : ~k = 1 γ Gt, ~ 2π. (2.44). and rewrite the signal formula: S(~k) =. ZZZ. . . ρ(~r)exp −i2π~k · ~r d~r.. (2.45). As it can be noticed the measured signal is conjugated with spin density through the Fourier transform. The inverse is also valid: ρ(~r) =. ZZZ. . . S(~k)exp i2π~k · ~r d~k. 24. (2.46).

(25) The crucial part is that the reciprocal space in which the signal is measured can be traversed either by change in time or change of gradient. Thus through precise manipulation of magnetic field gradients and time dependent signal detection one can sample whole k space in order to measure the desired spatially encoded signal.. 2.5. Diffusion tensor imaging. The phenomenological Bloch equations (2.37) neglect the influence of diffusion on the evolution of the net magnetization vector. A correction including diffusion terms was proposed by Torrey [11]:. ∂Mx ~ × B) ~ x − Mx + ∇ · D∇(Mx − Mx0 ) = γ(M ∂t T2 ∂My M ~ × B) ~ y − y + ∇ · D∇(My − My0 ) (2.47) = γ(M ∂t T2 ∂Mz ~ × B) ~ z − Mz − M0 + ∇ · D∇(Mz − Mz0 ), = γ(M ∂t T1 where: D is diffusion coefficient and Mi0 are components of equilibrium magnetization in given. direction. In order to make nuclear magnetic resonance experiment sensitive to diffusion one has to introduce diffusion sensitizing gradients. Consider simple spin-echo sequence [12]. In this sequence some time after exciting the sample with. π 2. an additional π pulse is applied. The evo-. lution of partial net magnetizations from particular voxels start with tipping them to transverse plane, they precess with Larmore frequency, after some time they are all flipped 180◦ about e.g. x axis, and precess further. Assuming the external field is perfectly homogeneous evolution of each of the partial magnetizations is exactly the same. However, if the π pulse is sandwiched between two strong gradient pulses along given direction, the evolution of the magnetizations differs along that direction. Those in the higher field region precess faster, whilst those in the lower field precess slower. If there is no movement along the direction of the diffusion gradient, the π pulse followed by second gradient pulse refocuses the net signal, i.e. all magnetizations get back in phase. However, if some nuclei change their positions between application of the diffusion gradient pulses, the magnetizations do not get back in phase perfectly. The stronger the diffusion, the less signal is refocussed with the second gradient. This difference in signal measured with and without diffusion gradient gives information about the diffusion coefficient. Such modified spin-echo sequence is called pulsed gradient spin echo (PGSE) [13]. The quantitative relation between the ratio of measured signals and the diffusion coefficient was given by Stejskal and Tanner [14]: Sx ln = −bD, (2.48) S0 where: Sx and S0 are signals measured with and without the diffusion gradient, respectively, b . . 25.

(26) is so called attenuation factor, and D is diffusion coefficient. In order to quantify diffusion in any direction, the simplest but quite robust is diffusion tensor model [15, 16]. A 3 × 3 second rank tensor diagonal elements describe magnitude of diffusion along three orthogonal axes of the laboratory reference frame. The off-diagonal elements indicate the rotation of diffusion ellipsoid, being visualization of the diffusion tensor, in respect to particular axes of laboratory frame. The explicit expression of the diffusion tensor is the following6 : . .  .  . Dxx Dxy Dxz. D = Dyx Dyy Dyz  .   Dzx Dzy Dzz. (2.49). In most cases, i.e. as long as uncharged molecules such as water are under investigation, the tensor is symmetrical: Dij = Dji ,. (2.50). thus consists of only six independent elements. The Stejskal-Tanner equation adapted to the diffusion tensor model becomes then: SB ln SB0. !. =−. 3 X 3 X. (Bij − B0ij )Dij ,. (2.51). i=1 j=1. where Bij is an element of the so called b-matrix which contains information about the diffusion gradients used while acquiring diffusion weighted signals and B0ij is an element of the b-matrix representing gradients active during acquisition of the reference image. Precise expression describing b-matrix elements depends on the shape of diffusion gradients. For example, in case of rectangular gradient pulses the ij element of the b-matrix is [17]: !. Bij = γ Gi Gj δ 2. 2. δ ∆− , 3. (2.52). where: Gi and Gj are the maximum field gradients along the ith and jth directions, respectively, δ is gradient duration and ∆ is time between the onset of the first and second gradient pulses. Therefore if the reference image is acquired without diffusion gradients the B0 matrix is equal to zero, and the number of B-matrices is determined by the number of diffusion gradient directions, i.e. at least six b-matrices describing six independent gradient directions are required to calculate the diffusion tensor. 6. Tensors and matrices are noted in bold throughout this thesis.. 26.

(27) 2.6. b-matrix spatial distribution. The principles of DTI described here above define how quantitatively measure diffusion tensor. However, the reality is a little more complicated. DTI inherent proneness to artefacts make the accuracy of the results questionable [18]. The problem stems from both various approaches to determination of the b-matrix and eddy currents induced by strong magnetic fields, disturbing the b-matrix. Several strategies were applied to work against the eddy current induced artefacts. Introducing bipolar diffusion sensitizing gradients was one of the predictive approaches [19]. It weakens the eddy current induced artefacts, but at the same time reduces the effective diffusion weighting [20]. Other predictive approaches focused on adjusting gradient sampling scheme which lessen the artefacts as well, but also vastly increase acquisition time [21–23]. On the other hand, the post processing strategies usually amount to either cross-correlation of the diffusion and baseline images [24], or calibration of the diffusion weighted data to the previously acquired images of a phantom [25]. An alternative approach claims that the commonly used approximate form of the b-matrix introduces systematic errors into the DTI experiment [26]. Thus in order to ameliorate the problem it was proposed to derive experimentally the spatial distribution of the b-matrix in order to calibrate the imaging sequence. The original BSD calibration procedure consists of the following steps: 1. Determination of the DTI sequence parameters and acquisition of the analysed specimen. 2. DTI acquisition of the isotropic phantom images. 3. DTI acquisition of the anisotropic phantom images in six independent positions. 4. Derivation of the b-matrix spatial distribution. 5. Derivation of the diffusion tensors.. 2.6.1. Normalization of the b effective. The first step is exactly the same as in the standard DTI procedure. Data acquisition is done in ordinary way, the only differences are the additional acquisitions of the phantoms and tensor derivation. The isotropic water phantom is used for normalization of the measured signal. Assuming the experiment is carried out in known temperature the self diffusion coefficient of water can be derived theoretically from the Stokes-Einstein equation for low Reynold’s number regime: Dt =. kB T , 6πηr. 27. (2.53).

(28) where: Dt is theoretically derived diffusion coefficient, η is water viscosity at given temperature, r is water molecule radius, and the remaining symbols has the same meaning as in previous equations. Thus in 21 ◦C, assuming water viscosity equal to 0.8 mPa s and water molecule radius to 1.375 Å, the self diffusion coefficient is roughly equal to 2 × 10−9 m2 /s. Since for isotropic phantom the diffusion ellipsoid should be symmetrical, the off-diagonal elements of the tensor are unimportant. The spatial orientation of the sphere does not matter due to its symmetry. From that it follows the off-diagonal elements of the b-matrix are also unimportant in this case, so only the b effective (Bef f ) is normalized to the theoretical value of the diffusion coefficient. Bef f = Bxx + Byy + Bzz .. (2.54). The normalization is done through the Stejskal-Tanner equation in a simplified form: Ni Bef f. 1 S = ln , Dt S0 . . (2.55). where superscript N i denotes we are dealing with the b effective derived through normalization to isotropic phantom. The normalization is done for each voxel in each slice and for each gradient direction, therefore the precise formula incorporating adequate indices7 is the following: Sg,s,m,n 1 ln . = Dt S0g,s,m,n !. Ni Bef f g,s,m,n. (2.56). The obtained map of spatial distribution of Bef f makes sense only within the region occupied by the isotropic phantom. Thus it is important to position the phantom inside the scanner in such a way to encompass the space occupied by region of interest of the given specimen with the phantom volume. Basing on such a simple normalization each element of the b-matrix can be spatially normalized in order to obtain spatial distribution of the b-matrix. This simplistic approach completely lacks anisotropic information, but even though can improve the results. The normalized diagonal elements of the b-matrix for each gradient direction, slice and voxel are given by the following equation: BiiN i. Ni Bef f = Bii · Bef f. (2.57). and the off-diagonal ones can be also normalized: BijN i = Bij. v u Ni uB t ii. Bii. ·. Ni Bjj , Bjj. (2.58). where: i 6= j. 7. g - gradient direction, s - slice, m - horizontal voxel coordinate for a given slice, n - vertical voxel coordinate. for a given slice.. 28.

(29) 2.6.2. Calibration of the b-matrix. Having determined Bef f from isotropic phantom one can derive the complete spatial distribution of each of the b-matrix elements. In order to do that six acquisitions of the anisotropic phantom placed in six various positions have to be acquired. Change of the position means rotations of the phantom in such a way the anisotropic core of the phantom is constantly occupying given volume of interest. The b-matrix elements for a given voxel can be then derived from the following system of six equations:. . ln. S0 Sl. . l l l l l l l l l l l l = Bxx Dxx + Byy Dyy + Bzz Dzz + 2Bxy Dxy + 2Bxz Dxz + 2Byz Dyz ,. (2.59). where l takes values from 1 to 6 and denotes either particular gradient direction in case of signal and b-matrix elements, or anisotropic phantom position in case of diffusion tensor elements. The above equation assumes there was no active diffusion gradient while acquisition of the reference signal S0 . The left hand side of the equation is known from the experiment, on the right hand side, the diffusion tensor elements are known for an anisotropic phantom with well defined structure, and the elements of the b-matrix can be derived. Let denote such element with superscript a for they are derived from the anisotropic phantom. Each b-matrix element for a given gradient direction and voxel is subsequently normalized to b effective from isotropic phantom: BijN a. =. Bija. Ni Bef · a f, Bef f. (2.60). where: i can be equal to j.. 2.6.3. Derivation of the calibrated tensor. Having specified all the b-matrix elements for each gradient direction, slice and voxel one can derive the diffusion tensor of the examined specimen. It can be done by substituting to the equation 2.59 the measured signal of the specimen and derived spatial distribution of the b-matrix, and solve the equation with respect to diffusion tensor elements. The procedure of tensor calculation is therefore almost the same as in the standard case, with this little difference, that for each voxel a dedicated b-matrix is incorporated into equations. On the contrary to the standard DTI it enables one to take into account spatial variability of the b-matrix and reduce present systematic errors [27].. 2.6.4. Convenient simplifications. Two significant simplifications can be introduced to the BSD calibration procedure. The former reduces the number of anisotropic phantom acquisitions by half under assumption of spatially 29.

(30) uniform diffusion weighting [28], and the latter simplifies the derivation of the spatial distribution of the b-matrix in case of high homogeneity of anisotropic phantom [29]. Spatially uniform diffusion weighting assumption A commonly accepted approach to DTI assumes the same diffusion weighting in each direction. The b-matrix can be expressed then as a product of a b effective and a gradient matrix constructed as a product of a unit gradient vector defining the direction with its transposition [30]: . gx2. gx gy gx gz. . ~G ~ T = Bef f ·  gy gx B = Bef f · G . gy2. gz gx gz gy.   . . gy gz  . (2.61). gz2. Such approach neglects the influence of the background and imaging gradients and omits their interaction with the diffusion sensitizing gradients. Nevertheless, the simplification reduces the number of unknown variables in the b-matrix from 6 to 3. In the forward calculation of the b-matrix the omission of the cross-terms and all the gradients besides diffusion ones, is one of the factors leading to introduction of systematic errors into the b-matrix. However, in a BSD approach, where the actual distribution of the b-matrix is measured, it only changes the meaning of the gradient matrix elements. They do not represent only the diffusion gradients any more, but rather an apparent diffusion gradient active during the sequence. This apparent gradient is the result of superposition of all the gradients and cross-terms active during the sequence. Separation of this apparent gradient into particular components would be challenging, but fortunately is not necessary for the calibration purposes. Phantom homogeneity assumption The crucial part of BSD-DTI is precise positioning of the anisotropic phantom inside the scanner. The spatial distribution of the b-matrix is calibrated in stationary volume inside the magnet, thus it is important to precisely ascribe phantom volume element to position in space. The operation has to be repeated after each rotation of the phantom. Nevertheless, if the phantom is homogeneous, i.e. its diffusion properties are spatially uniform, one can take just the average value of each of the diffusion tensor elements from the region occupied by the phantom and ascribe it to each voxel in a given phantom position. In reality phantom cannot be perfectly uniform, however, it was shown that if the relative standard deviation (RSD) of the distortion in phantom’s homogeneity does not exceed 0.2 %, the performance of the BSD under assumption of phantom homogeneity is almost as good as in the case of full BSD-DTI [29]. For bigger distortions the performance drops, however, as long as RSD is below 1 % there is still improvement in comparison to standard DTI [29].. 30.

(31) Further reading It is not possible to include precise references throughout theoretical description of the phenomenon discovered almost 100 years ago, successfully applied in number of fields of interest and described in number of books. Thus I restricted myself only to citation of the original papers where need be, and the textbooks I found most helpful in preparing this chapter are listed below in the alphabetical order. 1. Abragam J., Principles of Nuclear Magnetism [31] 2. Callaghan P. T., Principles of Nuclear Magnetic Resonance Microscopy [32] 3. Hecke W., Diffusion Tensor Imaging [33] 4. Keeler J., Understanding NMR Spectroscopy [34] 5. Levitt M. H., Spin Dynamics [35] 6. Mori S., Introduction to Diffusion Tensor Imaging [36] 7. Slichter C. P., Principles of Magnetic Resonance [37] 8. Szantay C., Anthropic Awareness [38]. 31.

(32) Chapter 3 Experiments La physique théorique est l’alliance de la physique sans l’expérience, et des mathématiques sans la rigueur.1 J. M. Souriau. 3.1. Phantoms. In order to experimentally validate BSD-DTI calibration one needs anisotropic phantoms. Number of various approaches to the construction of diffusion phantoms have been noted in the literature. Since in the clinical domain diffusion MRI studies focus on various neurological disorders in most cases the goal is to obtain a phantom with a well defined anisotropic structure and achieve as high fractional anisotropy (FA) as possible, in order to mimic neuronal tissue, which can achieve FA in order of 0.8 [39]. One of the most straight forward approaches to manufacturing phantoms of anisotropic biological structure of the neural tissue is to mimic it with another anisotropic biological structure. Asparagus turned out to be a useful model with well defined anisotropy [40]. However, the structure of such biological phantom is neither perfectly known nor reproducible from one specimen to another [41]. It is also likely to change over time, and the fractional anisotropy, even though uniform across the phantom, does not exceed 0.3 [42]. The artificial phantoms usually consist of hydrophobic fibres immersed in liquid. The diffusion properties of the phantom depends on the geometry and density of fibres packing. The higher the fibre density (FD) the higher the FA. A simple example of such a phantom may be just a bunch of fibres held tightly with a shrinking tube [43]. A more realistic phantom should distinguish with more precisely defined geometry and also can contain areas with crossing fibres [44]. There were also attempts to manufacture a diffusion phantom with thermal control in order to minimize the influence of the temperature on the measured diffusion coefficient [45]. 1. The theoretical physics is a combination of physics bereft of experience and mathematics bereft of rigour.. 32.

(33) Some of the phantoms are dedicated to calibrate a particular acquisition scheme or test specific reconstruction technique. It is mostly used in high angular resolution diffusion imaging (HARDI) together with the algorithms able to resolve crossing fibres within a single voxel, like Q-ball [46–48]. However, a caution should be taken in case of the phantoms which structure is precisely defined in scale smaller then the spatial resolution of the scanner due to possible artefacts caused by partial volume effect (PVE) [49]. An alternative approach to construction of diffusion phantoms is usage of capillaries [50] or hollow fibres obtained through electrospinning technique [51]. The latter if made from a hydrophobic material, requires saturation with liquid other than water. Such hollow fibre phantoms can be also adapted to organ specific design, like e.g. cardiac imaging [52]. On the other hand there is a need for a phantom which is not specific to a particular acquisition or reconstruction protocol, nor to the subject under investigation. A phantom which enables one a multi-centre comparison of the diffusion measurements. Example of a universal (high resolution) anisotropic DTI phantom may be phantom consisting of an array of channelled silicon plates [53]. The channels depths ranging from 50 µm to 300 µm obtained with high precision resulted in desired very small diffusivity across the channels. However, the diffusivity obtained along the channels turned out to be almost 50 % higher than free diffusion in the surrounding, what the authors called as an unexpected result needing explanation. This unpredictable effect can probably stem from the inhomogeneity of the b-matrix, but susceptibility artefacts can be also an issue.. 3.1.1. BSD-DTI phantoms. The BSD-DTI method is yet another, more profound, way to alleviate the inherent problems of DTI. Instead of a simple correction of the artefacts it focuses on the source of the problem i.e. systematic errors occurring in the spatial distribution of the b-matrix. The anisotropic diffusion phantoms constituent a core of this technique [54]. A complete set required for BSD calibration consist of an isotropic and an anisotropic phantom. Additionally another anisotropic phantom can be used for verification purposes. The isotropic phantom (F1) is simply a reservoir filled with liquid. The volume of the container has to be adjusted to the dimensions of the anisotropic phantom. The development of the BSD-DTI anisotropic phantoms was a long process with many trials and resulting with number of prototypes [55]. Various materials and geometries of anisotropic structures have been analysed [56]. Eventually, two geometries of the anisotropic phantoms emerged as the most useful. The former is a cubic plate (F2) phantom consisting of an array of thin glass plates separated with thin layers of liquid, whilst the latter is a cylinder capillary (F3) phantom in which held tightly together bundles of capillaries constituent a cylinder. The development of the phantoms can be divided into three stages, representing three generations of the phantoms. 33.

(34) Figure 3.1: MRI image of the F2 phantom. Copyright 2015 IEEE [56].. First generation The first generation isotropic phantom was a simple bottle of water of 3 cm diameter, i.e. large enough to encompass the anisotropic phantoms. The anisotropic phantom F2 distinguished with very high precision. The glass plates were very thin (100 µm), and the water layers between the plates were 20 µm thick. The glass-water array formed a cube of 2.5 cm side length (Fig. 3.1). The F3 phantom consisted of bundles of capillaries held tightly together and forming circular shape. Two types of the F3 phantom were manufactured, which differed in shape and size of capillaries, as well as, material of the capillary bundle. Polymethyl methacrylate (PMMA) was used in one case and silica glass in the other (Fig. 3.2). Each bundle constituted a matrix for 360 capillaries. The diameter of a single capillary was 36 µm and 49 µm for silica glass and PMMA, respectively. The diameter of a whole bundle was 920 µm in case of silica glass and 760 µm for PMMA bundle. The PMMA bundles were found more likely to deform at the edges and to have bigger spread in diameter of capillaries within a bundle. The silica glass bundles have, however, thicker walls and therefore the active volume of the phantom is reduced. The outer diameter of the F3 phantom was 1.5 cm and height was 2 cm (Fig. 3.3). Second generation The second generation consisted of similar set of phantoms, however, adapted to clinical scanner with significantly bigger bore. The F1 isotropic phantom was a bottle of water treated with copper sulphate. The bottle had 12 cm in diameter and was 18 cm high. The anisotropic phantoms were based on the first generation, however, small modifications have been done. The outer dimensions were enlarged. The F2 phantom side length was increased to 5 cm and the F3 phantoms diameters changed to 2.5 cm and their height to 3 cm.. 34.

(35) Figure 3.2: Bundles of capillaries made of different materials (left: PMMA, right: silica glass). Copyright 2015 IEEE [56].. Figure 3.3: Cross section of the F3 phantom. Bundles of capillaries are visible. Copyright 2015 IEEE [56].. 35.

(36) The inner diameter of single PMMA capillary was reduced to 40 µm. One of the sources of severe artefacts in high magnetic fields are differences in magnetic susceptibility between neighbouring regions. Local strong magnetic field gradients appear on such interfaces and cause distortions of the images. In order to reduce susceptibility artefacts the anisotropic phantoms were covered with 3 mm thick glass. Thanks to that, first trials of diffusion imaging with usage of echo planar imaging (EPI) sequence became possible. Due to very high cost of the first generation F2 phantom, the glass plates were thicken to 150 µm in order to reduce the cost, however 20 µm liquid layer spacing was kept. Third generation The phantoms of third generation met with three major requirements; further reduction of the susceptibility artefacts, easier and more precise phantom positioning, compatibility with a clinical scanner head coil. The thick glass cover used in the second generation anisotropic phantoms gave not fully satisfying results in case of EPI. The susceptibility artefacts were still present. In order to reduce the unwanted artefacts, the phantom cores were placed inside significantly bigger glass balls filled with liquid (Fig. 3.4). Such solution keeps the problematic air-glass-water interface far enough from the phantom core. The remaining issue is the liquid-glass interface all around the phantom core. Since volume susceptibility of water is slightly lower than volume susceptibility of glass (−9.021 × 10−6 and −8.470 × 10−6 for water and glass, respectively) addition of small amount of substance with positive susceptibility to the water, can level the difference. Luckily commonly used in MRI phantoms copper sulphate has positive volume susceptibility (3.773 × 10−4 ). Theoretically derived mass fraction of CuSO4 levelling the difference in susceptibilities between water and glass is equal to 0.51 % [55]. However, CuSO4 also visibly shortens the T2 of the treated liquid, what is very undesired effect in diffusion measurements, which are based on measuring the signal loss due to the diffusion gradients. A compromise was achieved by carrying out DTI measurements for the dummy phantom with various amounts of CuSO4 additions. Analysis of the signal profiles going through liquid-glass, liquid-PMMA and air-glass interfaces established the optimum mass fraction of CuSO4 between 0.1 % and 0.2 % [55, 56]. New external geometry of the phantoms not only reduced the susceptibility artefacts, but also made precise positioning easier. Circular basis of the phantom fits the goniometer compatible with the clinical scanner head coil. The goniometer has two degrees of freedom: rotation about the axis along the scanner (Z axis) and about the axis going from top to the bottom of the phantom. It suffice to position the phantom in any desired orientation. During the rotations phantom’s centre remains in the isocentre of the head coil (Fig. 3.5). Usage of PMMA capillaries was abandoned due to better quality of the silica glass capillaries. The diameter of the silica glass capillary can be smaller and imperfections in structure are less likely to occur. Detailed phantoms’ specification was presented in table 3.1. 36.

(37) (a) Isotropic. water. reference. phantom (F1).. (b) Anisotropic cubic plate phan-. (c) Anisotropic cylinder capilary. tom (F2).. phantom (F3).. Figure 3.4: Complete set of third generation phantoms enabling calibration normalization and validation of BSD-DTI measurements.. Table 3.1: Third generation phantoms specification.. Parameter. Value. Glass ball radius. 5.3 cm. F2 core dimensions. 5.0 cm × 5.0 cm × 5.0 cm. F2 plate dimensions. 5.0 cm × 5.0 cm × 0.18 cm. F2 water layer thickness. 20 µm. Number of plates in F2 core. 236. F2 phantom weight. 368.15 g. F3 core height. 3.0 cm. F3 core diameter. 3.0 cm. Capillary bundle diameter in F3 core. 830 µm. Single capillary diameter in F3 core. 30 µm. Number of bundles in F3 core. 921. Number of capillaries in a bundle. 300. F3 phantom weight. 351.43 g. 37.

(38) (a) Initial position.. (b) Position after 90◦ rotation about Z axis.. Figure 3.5: The F2 phantom attached to the goniometer placed inside clinical tomograph’s head coil.. 3.2. BSD-DTI experiments. Three experimental sessions have been carried out on three different MRI scanners. During each session one of the phantom generations was used to obtain DTI data. The sessions differed in details, but every time collected data were sufficient to perform BSD calibration. In chronological order the experimental sessions were performed on 9.4 T Bruker research scanner (first generation phantoms), 3.0 T GE clinical scanner (second generation phantoms) and 3.0 T Siemens clinical scanner (third generation phantoms).. 3.2.1. Bruker 9.4 T. The initial experiments were carried out in the Henryk Niewodniczański Institute of Nuclear Physics Polish Academy of Science in Cracow. The imaging was performed on 9.4 T Bruker Biospec 94/20 scanner with 210 mm bore diameter2 . The sequence used was basic PGSE. Sequence parameters were kept the same in each experiment. List of the most important imaging sequence parameters is collected in table 3.2. The reference data were collected without diffusion weighting. Gradient encoding scheme was set to default for this particular system (Tab. 3.3). The data put under analysis in this thesis have been acquired in axial direction. There were two acquisitions of the F1 phantom, single measurement of the F2 phantom in six different positions and one acquisition of the F3 phantom. For further convenience the data will be referred to through the following abbreviation: TOFX, where T is scanner make (i.e. B for Bruker, G for GE and S for Siemens), O is the orientation (A - axial, S - sagittal, C - coronal), 2. The actively shielded gradient coils installed inside the magnet narrow the effectively available bore diameter. to 120 mm.. 38.

(39) Table 3.2: Bruker sequence parameters.. Parameter. Value. Imaging sequence. PGSE. FOV. 40 mm × 40 mm. Matrix size. 64 × 64. Repetition time. 2500 ms. Echo time. 40 ms. Slice thickness. 2 mm. b value. 600 s/mm2. Number of slices. 5. Number of gradient directions. 6. Number of averages. 4. Table 3.3: Bruker diffusion gradient vectors.. Gradient direction. vx. vy. vz. 0. 0.000. 0.000. 0.000. 1. 0.667. 0.333. 0.667. 2. 0.667. -0.333. 0.667. 3. 0.333. 0.667. 0.667. 4. -0.333. 0.667. 0.667. 5. 0.667. 0.667. 0.333. 6. 0.667. 0.667. -0.333. Presented values are rounded. In calculations 6 digit precision was kept.. FX denotes phantom type (F1, F2, F3 according to previous description of the phantoms). Additionally if the acquisition was repeated the acquisition number may appear after underscore (pertains to the F1 phantoms), and for anisotropic phantoms additional number following the letter ’P’ denotes position of the phantom. Thus for example data obtained on Bruker scanner in axial orientation with usage of isotropic phantom may be denoted: BAF1_1. Complete list of the data acquired on Bruker scanner was presented in table 3.10.. 3.2.2. GE 3.0 T. Second experimental session was carried out in Voxel Medical Diagnostic Centre in Cracow on 3.0 T GE Discovery MR 750 clinical scanner. The anisotropic second generation phantoms were immersed in a bucket filled with water in order to further reduce the susceptibility mismatch. The data were collected with spin echo echo planar imaging (SE-EPI) sequence (see table 3.4 for sequence details). Gradient encoding scheme was set default for this scanner (Tab. 3.5). Complete list of the acquired data is presented in table 3.11. 39.

(40) Table 3.4: GE sequence parameters.. Parameter. Value. Imaging sequence. SE-EPI. FOV. 240 mm × 240 mm. Matrix size. 256 × 256. Repetition time. 6000 ms. Echo time. 83 ms. Slice thickness. 3 mm. b value. 1000 s/mm2. Number of slices. 25. Number of gradient directions. 6. Number of averages. 4. Table 3.5: GE diffusion gradient vectors.. Gradient direction. vx. vy. vz. 0. 0.000. 0.000. 0.000. 1. 1.000. 0.000. 0.000. 2. 0.446. -0.895. 0.000. 3. 0.447. -0.275. 0.851. 4. 0.448. 0.723. -0.525. 5. 0.447. 0.724. 0.525. 6. -0.449. 0.277. 0.850. Presented values are rounded. In calculations 6 digit precision was kept.. 40.

(41) Table 3.6: Siemens sequence parameters.. Parameter. Value. Imaging sequence. SE-EPI. FOV. 244 mm × 244 mm. Matrix size. 122 × 122. Repetition time. 5300 ms. Echo time. 95 ms. Slice thickness. 2 mm. b value. 1000 s/mm2. Number of slices. 32. Number of gradient directions. 6. Number of averages. 32 (F1 phantom) 12 (F2 and F3 phantoms). Table 3.7: Siemens diffusion gradient vectors.. Gradient direction. vx. vy. vz. 0. 0.000. 0.000. 0.000. 1. 0.709. -0.005. -0.705. 2. -0.709. 0.005. -0.705. 3. 0.001. 0.711. -0.703. 4. 0.001. 0.711. 0.703. 5. 0.706. 0.708. -0.001. 6. -0.704. 0.710. -0.001. Presented values are rounded. In calculations 6 digit precision was kept.. 3.2.3. Siemens 3.0 T. Third experimental session was carried out on 3.0 T Siemens Skyra clinical scanner in John Paul the Second hospital in Cracow. The third generation phantoms were imaged with usage of SE-EPI sequence (Tab. 3.6) and gradient encoding scheme was set to default for this system (Tab. 3.7). The complete list of the acquired data is presented in table 3.12. The same Siemens 3.0 T scanner was also used for the purpose of soft tissue imaging. The experimental session took place some time after imaging of the phantoms. In the meantime the software was upgraded and the SE-EPI Resolve sequence became available. A young adult healthy male volunteer served as a patient and gave written informed consent. The imaging sequence parameters were collected in table 3.8. Since the sequence parameters and the diffusion gradient encoding scheme differed from the one used in previous study of the phantoms (table 3.9), another acquisition of the isotropic calibration phantom images was carried out. Due to limited access to the scanner being constantly in use for clinical studies and the conclusions 41.

(42) Table 3.8: Soft tissue imaging sequence parameters.. Parameter. Value. Imaging sequence. SE-EPI Resolve. FOV. 300 mm × 300 mm. Matrix size. 160 × 160. Repetition time. 6100 ms. Echo time. 65 ms. Slice thickness. 4 mm. b value. 1000 s/mm2. Number of slices. 27. Number of gradient directions. 6. Number of averages. 4. Table 3.9: Soft tissue imaging diffusion gradient vectors.. Gradient direction. vx. vy. vz. 0. 0.000. 0.000. 0.000. 1. -0.713. 0.187. 0.676. 2. 0.697. 0.091. 0.711. 3. 0.058. 0.553. -0.831. 4. 0.043. 0.830. 0.556. 5. 0.756. 0.644. -0.120. 6. -0.655. 0.740. -0.155. Presented values are rounded. In calculations 6 digit precision was kept.. drawn from the previous phantom imaging session, the acquisition of the anisotropic phantom images was abandoned in this case.. 3.3. Data collections, datasets, experiments. The data analysed in this thesis comprise of four smaller subsets named data collections. Each data collection represent all measurements acquired on a particular scanner using particular generation of the phantoms or measurements of the soft tissue and the isotropic phantom. Each data collection consists of experiments being the most atomic (from further analysis perspective) part of the data structure. Single experiment represent a DTI measurement of a given specimen, i.e. one reference acquisition without diffusion gradients and six acquisitions with various diffusion gradient directions. Each acquisition within the experiment may be repeated several times in order to improve the signal to noise ratio (SNR). By choosing particular experiments from given data collection one can create the dataset. 42.

(43) Table 3.10: Bruker data collection. Description. Phantom. Orientation. Position. Acquisition number. BAF1_1. F1. Axial. -. 1. BAF1_2. F1. Axial. -. 2. BAF2P1. F2. Axial. 1. -. BAF2P2. F2. Axial. 1. -. BAF2P3. F2. Axial. 1. -. BAF2P4. F2. Axial. 1. -. BAF2P5. F2. Axial. 1. -. BAF2P6. F2. Axial. 1. -. BAF3P1. F3. Axial. 1. -. The dataset is a collection of experiments required for a BSD-DTI calibration. In the most simple case of isotropic normalization of the b effective the dataset consists of only two experiments: acquisition of the analysed sample in a given orientation and acquisition of the isotropic phantom in the same orientation. Another option is a dataset consisting of five experiments; beside the two from the previous case, acquisitions of the anisotropic phantom in three orthogonal positions are also included. In the case of full BSD-DTI (without simplifications) six experiments with anisotropic phantom oriented in six different positions should be considered. However, such a case is not analysed in this thesis, since it was shown, the simplified approach is accurate enough [29].. 43.

(44) Table 3.11: GE data collection. Description. Phantom. Orientation. Position. Acquisition number. GAF1_1. F1. Axial. -. 1. GAF1_2. F1. Axial. -. 2. GSF1_1. F1. Sagittal. -. 1. GSF1_2. F1. Sagittal. -. 2. GCF1_1. F1. Coronal. -. 1. GCF1_2. F1. Coronal. -. 2. GAF2P1. F2. Axial. 1. -. GAF2P2. F2. Axial. 2. -. GAF2P3. F2. Axial. 3. -. GSF2P1. F2. Sagittal. 1. -. GSF2P2. F2. Sagittal. 2. -. GSF2P3. F2. Sagittal. 3. -. GCF2P1. F2. Coronal. 1. -. GCF2P2. F2. Coronal. 2. -. GCF2P3. F2. Coronal. 3. -. GAF3P1. F3. Axial. 1. -. GAF3P2. F3. Axial. 2. -. GAF3P3. F3. Axial. 3. -. 44.

(45) Table 3.12: Siemens data collection. Description. Phantom. Orientation. Position. Acquisition number. SAF1_1. F1. Axial. -. 1. SAF1_2. F1. Axial. -. 2. SSF1_1. F1. Sagittal. -. 1. SSF1_2. F1. Sagittal. -. 2. SCF1_1. F1. Coronal. -. 1. SCF1_2. F1. Coronal. -. 2. SAF2P1. F2. Axial. 1. -. SAF2P2. F2. Axial. 2. -. SAF2P3. F2. Axial. 3. -. SSF2P1. F2. Sagittal. 1. -. SSF2P2. F2. Sagittal. 2. -. SSF2P3. F2. Sagittal. 3. -. SCF2P1. F2. Coronal. 1. -. SCF2P2. F2. Coronal. 2. -. SCF2P3. F2. Coronal. 3. -. SAF3P1. F3. Axial. 1. -. SAF3P2. F3. Axial. 2. -. SAF3P3. F3. Axial. 3. -. SSF3P1. F3. Sagittal. 1. -. SSF3P2. F3. Sagittal. 2. -. SSF3P3. F3. Sagittal. 3. -. SCF3P1. F3. Coronal. 1. -. SCF3P2. F3. Coronal. 2. -. SCF3P3. F3. Coronal. 3. -. 45.