Abib

Kostia, Isa Samuel

Federica

Amandine Mauro Patrick

Gabriel

Daniela Sara

James Mclauren

Marco, Konstantinos, Maureen James Inglis

Anna Vendri

Tia Ketty

Pippo, Cesco Bebe

Diffusion and Perfusion Magnetic Resonance Imaging (dMRI & pMRI) represent two modalities that allow sensing important and different but complementary aspects of brain imaging.

This thesis presents a theoretical and methodological investigation on the MRI modalities based on diffusion-weighted (DW) and dynamic susceptibility contrast (DSC) images. For both modalities, the contributions of the thesis are related to the development of new methods to improve quality, processing, and exploitation of the obtained signals.

With respect to contributions in diffusion MRI, the nature of the complex DW signal is investigated to explore a new potential contrast related to tissue microstructure. In addition, the complex signal is exploited to correct a bias induced by acquisition noise of DW images, thus improving the estimation of structural scalar metrics.

With respect to contributions in perfusion MRI, the DSC signal processing is revisited in order to account for the bias due to bolus dispersion. This phenomenon prevents the correct estimation of perfusion metrics but, at the same time, can give important insights about the pathological condition of the brain tissue.

The contributions of the thesis are presented within a theoretical and methodological framework, validated on both synthetic and real images.

Résumé (en français)

Le travail effectué dans cette thèse touche aux modalités d'IRM basées sur des images pondérées en diffusion (IRMd) et sur des images de perfusion (IRMp) par injection de produit de contraste. Pour chacune des deux modalités, des contributions computionnelles, aussi bien théoriques que méthodologiques, sont développées pour améliorer la qualité, le traitement et l'exploitation des signaux acquis.

Ces deux modalités sont basées sur des mécanismes de contraste complémentaires de l'imagerie cérébrale: sur la base du mouvement brownien des molécules d'eau dans le cerveau, l'IRMd permet de remonter à une information sur la connectivité structurale cérébrale, alors que l'IRMp permet d'inférer le flux sanguin qui s'y trouve, ce qui peut être lié au degré de certains métabolismes tissulaires et d'états pathologiques.

Cette thèse est structurée en 3 parties principales.

Partie I: "La nature du signal"

Cette partie est fondamentalement constituée par le chapitre 3 qui pose les bases de la résonance magnétique nucléaire (RMN), avec un accent particulier sur le contraste pondéré en diffusion et sur celui basé sur la susceptibilité dynamique, c'est-à-dire les principes de l'IRMd et l'IRMp. Le chapitre fait un relevé des concepts fondamentaux qui jouent un rôle dans la génération du signal: les acquisitions, comme la Pulsed Gradient Spin Echo, et les phénomènes physiques, tels que le mouvement brownien de particules. Les chapitres d'introduction à l'IRM de diffusion (chapitre 4) et de perfusion (chapitre 9) font référence aux concepts présentés dans ce chapitre.

Part II: "Diffusion"

Cette partie porte sur l'IRM pondérée en diffusion (IRMd). Après un premier chapitre d'introduction, les trois suivants étudient la nature complexe du signal de diffusion, et son lien avec le phénomène physique de diffusion. Un dernier chapitre est ensuite consacré à l'exploitation du signal complexe dans le cadre d'un processus de débruitage des images DWI par correction de phase.

Chapitre 4 illustre les bases du signal de diffusion, tels que le concept d'atténuation, sa dépendance sur les paramètres contrôlables d'une séquence PGSE, et son in-terprétation en termes de diffusion dans le tissu micro-structurel. Le chapitre présente ensuite la dérivation de la relation de Fourier entre le signal et la densité de probabilité moyenne des déplacements des particules, connue comme l'EAP (Ensemble Average Propagator), un concept fondamental au coeur des trois chapîtres suivants. Enfin, le chapitre décrit deux modèles bien connus pour représenter les signaux de diffusion et qui seront utilisés dans le chapitre 8: le tenseur de diffusion (DTI), et le propagateur apparent moyen (MAP-MRI).

Chapitre 5 explore la nature complexe de l'atténuation due à la diffusion. En particulier, le chapitre présente des exemples de la façon dont l'atténuation complexe peut être liée au processus de diffusion sous-jacent dans certaines conditions et hypothèses spécifiques. Au moyen de ces exemples, le concept d'asymétrie de l'EAP sera introduit et clarifié, et étudié plus en détail dans les deux chapitres suivants. Ce chapitre contribue particulièrement à la formalisation de l'idée d'asymétrie de l'EAP et à la clarification des hypothèses derrière l'existence de l'atténuation complexe.

Chapitre 6 propose une nouvelle forme d'exploitation de la relation de Fourier basée sur les principes et hypothèses décrits dans le chapitre précédent. Pour cela, on propose un exemple de géométrie axonale qui montre une asymétrie intrinsèque. L'EAP asymétrique est dérivée analytiquement pour des axones ondulés et comprimés. Ensuite, son asymétrie -liée à l'atténuation complexe -est quantifiée pour caractériser la quantité de compression axonale dans laquelle les particules d'eau diffusent. La contribution de ce chapitre réside dans la proposition d'un paradigme géométrique où l'asymétrie de l'EAP est étroitement liée aux caractéristiques de la microstructure du tissu sous-jacent.

Chapitre 7 propose un paradigme différent où l'asymétrie de l'EAP est induite par les conditions expérimentales, précisément en imposant un gradient de diffusivité auquel les particules d'eau sont soumis. En particulier, un fantôme théorique est développé qui rend compte des paramètres expérimentaux avec différents types de tissus cérébraux. Dans cette étude théorique, le gradient de diffusivité -qui causerait la diffusion asymétrique -est imposé par un gradient de température. L'EAP est dérivé analytiquement pour un tel scénario. In-fine, ce chapitre contribue à clarifier la dépendance de l'asymétrie de l'EAP avec les conditions expérimentales.

Chapitre 8 mène une exploitation du signal complexe dans le cadre d'un processus de débruitage des images DWI par correction de phase. Avec cette procédure, le problème de la gestion du bruit Rician qui affecte l'amplitude des données DWIs peut être évité tout en menant une caractérisation et une quantification rigoureuse des effets de la correction de phase sur le signal de diffusion, et sur les métriques scalaires calculées avec les modèles DTI et MAP-MRI (espace-q). Une contribution complémentaire a été la définition d'une stratégie pour trouver la quantité optimale de correction de phase nécessaire pour chaque donnée DWI. Cela permet d'éviter le réglage empirique des paramètres qui pourraient se traduire par une corruption plutôt qu'une correction des données DWIs. Le résultat est un cadre de correction de phase plus général dans lequel aucun réglage n'est nécessaire et où des caractéristiques spécifiques du bruit affectant les images, telles que la variabilité spatiale de sa variance, peuvent être prises en compte. La correction de phase proposée est ensuite validée sur des images synthétiques et réelles. Au final de cette seconde partie dédiée à l'IRMd, les contributions computationnelles proposées et validées sur plusieurs jeux de données synthétiques et réelles, constituent un cadre bien adapté pour le prétraitement des images DWI tout en y apportant des solutions pratiques dans le domaine.

Partie III: "Perfusion"

Cette partie est axée sur l'IRM de perfusion, et en particulier sur sa modalité de contraste de susceptibilité dynamique (DSC). Après un premier chapitre introductif, où les bases de la théorie de la perfusion sont révisées, l'accent est mis sur les concepts liés à la caractérisation d'un phénomène connu sous le nom de dispersion du bolus, qui polarise les cartes scalaires du flux sanguin cérébral (CBF) et du temps moyen de transit (MTT) qui constituent des résultats importants en IRMp. Cependant, cette information de dispersion est potentiellement liée à la présence de maladies steno-occlusives et peut donc fournir des informations importantes. Pour cette raison, dans les chapitres suivants il est proposé un cadre de traitement de données de perfusion qui se base sur des bases fonctionnelles (DCB) afin d'intégrer cette présence éventuelle de dispersion du bolus et arriver à la caractériser. La méthode est validée sur des images synthétiques et réelles. Le cadre présenté, qui évite la modélisation explicite de la dispersion, se révèle utile en tant qu'étape de prétraitement du signal afin d'améliorer les résultats de techniques basées sur des modèles explicites. Le dernier chapitre utilise le cadre proposé pour inférer in vivo la nature de la dispersion, une procédure originale et efficace qui pourrait ouvrir de nouveaux scénarios potentiels de caractérisation et de quantification de la dispersion à des fins cliniques.

Chapitre 9 décrit les bases de la théorie sur laquelle repose la modélisation de la perfusion. En particulier, on introduit les paramètres de perfusion d'intérêt, tels que le flux sanguin cérébral (CBF) et le temps moyen de transit (MTT). Ensuite, on introduit le concept fondamental de la fonction de transfert du tissu, R(t), une fonction inconnue, dépendant du temps, qui décrit entièrement la perfusion dans un voxel, et qui est généralement estimée par déconvolution. Le chapitre présente également les défis liés à l'estimation de la fonction de réponse, abordés dans les chapitres suivants, tels que la question du retard, et la dispersion de bolus.

Chapitre 10 présente la modélisation mathématique de la dispersion du bolus, et son interprétation phénoménologique. On y rapporte aussi les effets de la dispersion sur la fonction de transfert et sur les paramètres de perfusion estimées. A cette fin, le chapitre introduit la notion de temps de dispersion et clarifie aussi le concept selon lequel, en présence de dispersion, la déconvolution conduit à l'estimation de paramètres de perfusion et de fonction de transfert qui sont de perfusion effective par opposition à une estimation sans dispersion. Cette phase d'estimation sera effectuée avec la méthode DCB qui est proposée dans le chapitre suivant.

Chapitre 11 constitue le coeur méthodologique de la partie de perfusion. Il propose une description de la fonction de transfert effective au moyen de bases fonctionnelles "Dispersion-Compliant Bases (DCB)": celles-ci sont proposées dans des versions non linéaires et linéarisées. Le chapitre présente également des détails sur les plus récentes techniques qui sont utilisées dans cette thèse pour effectuer des comparaisons. Après avoir proposé la déconvolution DCB, il est également présenté une méthode pour récupérer les paramètres de perfusion sans dispersion. Le cadre sera validé dans le chapitre qui suit.

Chapitre 12 présente des expériences pour valider in silico et in vivo le cadre de déconvolution DCB introduit dans le chapitre précédent. On montre que la méthode DCB a une plus grande précision, par rapport aux techniques de l'état de l'art dans l'estimation de la fonction de transfert effective et des paramètres de perfusion. De plus, il est montré que l'utilisation de la déconvolution DCB en tant qu'étape de prétraitement du signal améliore l'estimation de la fonction de transfert et des paramètres de perfusion sans-dispersion, quand ceux-ci sont estimés via des méthodes utilisant un modèle de dispersion.

Chapitre 13 propose une méthode permettant d'inférer in vivo le modèle correct à adopter pour représenter la dispersion, afin d'obtenir une estimation correcte de la fonction de réponse et des paramètres de perfusion sans dispersion. En même temps, le chapitre illustre l'importance d'effectuer une déconvolution de perfusion avec peu d'hypothèses concernant le processus de dispersion éventuellement en cours, par ex. sans adopter un modèle de dispersion au cours de la déconvolution, tel que proposée dans l'approche DCB.

Partie IV: "Conclusion"

Cette partie, constituée par le chapitre 15, présente le résumé des contributions pour la partie II (diffusion) et pour la partie III (perfusion) et conclut par quelques perspectives liées aux contributions proposées dans ce travail.

Life is a wonderful journey. And I had the luck of making a three years long stop in Athena. Actually, this stop was a journey itself that brought me around the world, where I met incredibly interesting people, friends, and had unique experiences.

For this great opportunity, but not only, I deeply thank Rachid. I remember when he showed me the view on the forest from the balcony of the 5 th floor... He knows how to attract people here!! He believed in me, and when I finally arrived he said "Marco, here you are in family", and it really felt like so. Even now, three years from then, I see day by day how he puts my happiness before his own interests. Rachid, I thank you for being such a great human person who believes in people.

My time here wouldn't have been as good without all the friends I made in the team. Starting with Ruggy-Ruggy, the person I have been sharing most of my days with.. and not just because we had to share the office.. We traveled around the world, spent countless hours talking about life, relationships, emotional crisis.. and even science sometimes!! Rutger, I can't thank you enough for being the supportive and understanding shoulder I could rely on during these years.. Plus, let me tell you how much I could appreciate your genuine and contagious passion for enjoying life.. As you said, let's keep on drinking beers together after this PhD experience!! Guillermo Galladro was not just my friend and colleague in these years.. but also a flatmate despite he had his own apartment.. Guishi, after our countless chatsagain women and life -I always had the impression that I have to learn how to be mature from you, despite I like to proudly think of you as my younger brother. The more I write the more I realize that I have been just sharing a lot with many people. Stefano BigHouse, you are a great friend, and you managed to make me feel like nella mia bella Italia, sharing thoughts, moods, and values that are precious to me. My dear friend Auro, who pushed me through my first steps in the PhD, basically setting the path for me. Never been more honored of being the "workslave" of somebody. Nathanael and Sebastian were my great philosophers.. I loved spending time talking with you.. Brahim, "Abibi".. I don't know how I will do without your daily "Hola amigo".. and your wisdom.. Christos, a genuine nice true person with too less time due to kitesurf.. but always full of emotions in our chats.. Nathalie, thank you for being such a (half) Greek force of nature and kind at the same time.. Kai, a sweet person, big heart, a pleasure to meet during the days.. unforgettable mushrooms on the bbq.. Demian, I should really take inspiration from your charisma.. thank you for pushing, and pushing, and pushing me to believe in crazy ideas.. and for the positivity you bring around yourself.. I owe a lot to the friendship with Rodrigo such as great fun moments, a bbq on my third day in the region, half of my social life, and much more.. I have a nice thought for everybody Introduction 21 Magnetic resonance imaging (MRI) is a versatile technique with the unique capability of producing images having diverse contrast mechanisms. This thesis is focused on theoretical and methodological aspects related to the processing of diffusion and perfusion MRI (dMRI & pMRI) data in the brain. In fact, although their contrast mechanisms are directly or indirectly related to the Brownian motion of water particles, these two modalities represent two different and complementary sides of brain imaging: dMRI gives a structural information about the brain tissue, whereas pMRI allows inferring on the blood flow therein, which can be related to the grade of tissue metabolism and to pathological conditions. This thesis aims at investigating in both directions. The contributions of this thesis are organized and structured in three main parts.

Part I: "Common Grounds: The Nature of the Signal"

This part is fundamentally constituted by chapter 3 which reports the basics of nuclear magnetic resonance (NMR), with a special focus on diffusion-weighted and dynamic susceptibility contrast NMR, i.e. the principles of the corresponding MRI modalities. The chapter makes a survey of the fundamental concepts that play a role in the signal generation: the acquisitions, such as Spin-Echo and the Pulsed Gradient Spin Echo, and the physical phenomena, such as the Brownian motion of spin bearing particles. Both the next introductory chapters about diffusion (chapter 4) and perfusion (chapter 9) MRI will heavily refer to the concepts presented in this chapter.

Part II: "Diffusion"

This part focuses on the signal contrast based on diffusion-weighted NMR. After a first introductory chapter, the three following ones will investigate on the complex nature of the diffusion signal, and on its link with the underlying diffusive physical phenomenon. The last chapter, instead, will present a technique that exploits the complex nature of the diffusion-weighted images (DWIs) to eliminate a severe bias in the conventionally used magnitude DWI, which is related to the non-Gaussian nature of the noise.

Chapter 4 illustrates the basics of the diffusion signal, such as the concept of attenuation, its dependency on the controllable parameters of a PGSE sequence, and its interpretation in terms of the underlying diffusion within tissue microstructure. The chapter then reports the derivation of the Fourier relationship between the signal and the average probability density function of the spins' displacements, namely the ensemble average propagator (EAP), which is a fundamental concept at the heart of the next three chapters. Finally, the chapter describes two signal representations that will be used in chapter 8. These are the diffusion tensor imaging (DTI), and the mean apparent propagator (MAP-MRI), with the corresponding scalar metrics.

Chapter 5 explores the nature of the complex attenuation due to diffusion, trying to shed light on the subject which for many aspects is controversial. Particularly, the chapter presents examples of how the complex attenuation can be linked to the underlying diffusion process under some specific conditions and assumptions. By means of these examples, the concept of EAP asymmetry will be introduced and clarified, to be further investigated in the next two chapters. This chapter particularly contributes in the formalization of the idea of EAP asymmetry and in clarifying assumptions behind the existence of the complex attenuation.

Chapter 6 proposes a new form of exploitation of the Fourier relationship based on the principles and assumptions described in the previous chapter. In fact, it proposes an example of axonal geometry that shows intrinsic asymmetry. Particularly, the asymmetric EAP for undulated and compressed axons is derived analytically, and then its asymmetry -linked to the complex attenuation -is quantified to characterize the amount of axonal compression within which water spins undergo diffusion. The contribution of this chapter lies in the proposition of a geometrical paradigm where the asymmetry of the EAP is closely linked with features of the underlying microstructure.

Chapter 7 proposes a different paradigm where EAP asymmetry is induced by the experimental conditions, precisely by imposing a diffusivity gradient to which spins are subject to. Particularly, a theoretical phantom is developed which accounts for different experimental parameters and brain tissue types. In this theoretical investigation, the diffusivity gradient -that would cause asymmetric diffusion -is imposed by a temperature gradient. Also in this case the EAP is derived analytically for such a scenario. The chapter contributes in further clarifying the dependence of the potential EAP asymmetry on the experimental conditions.

Chapter 8 makes a change of perspective and exploits the complex nature of the DWIs in order to study how the Rician bias of magnitude DWIs -which are commonly used in the current practice -can be avoided by means of phase correction. Particularly, the chapter contributes by characterizing the debiasing effects of phase correction on the signal, and on DTI and q-space metrics, e.g. those calculated with MAP-MRI. The proposed analysis reveals that phase correction is potentially useful despite it is affected by some procedural limitations. In particular, one current limitation is represented by the need of having to manually tune the "amount" of desired phase correction, which is inefficient and might lead to data alterations. The chapter then contributes by proposing a strategy to improve over the current technique, thus overcoming this limitation. The result is a more general framework where no tuning is required and where specific characteristics of the noise affecting DWIs, such as the spatial variability of its variance, can be taken into account. The proposed phase correction is then validated in silico and in vivo.

Part III: "Perfusion"

This part focuses on theory and processing related to perfusion MRI, and particularly to the dynamic susceptibility contrast (DSC) modality. After a first introductory chapter, where the basics of perfusion theory are illustrated, the focus is on the concepts related to the characterization of a phenomenon known as bolus dispersion. Dispersion constitutes a bias for the estimation of scalar maps of the cerebral blood flow (CBF) and the mean transit time (MTT), which are the main outcomes of perfusion data processing. However, this bias is also potentially related to the presence of steno-occlusive diseases and can therefore provide important information. For this reason, the next chapter 11 proposes a perfusion data processing framework, based on proposed dispersion-compliant bases (DCB), that embraces the eventual presence of bolus dispersion and characterizes it. The method is validated with synthetic experiments and on real data. Moreover, the presented framework -that avoids explicit modeling -reveals useful as a pre-processing step for model-based techniques that can quantify dispersion, improving their results. The final chapter uses the proposed framework to infer in vivo on the nature of dispersion, an innovative procedure that could open new potential scenarios in dispersion characterization and quantification for clinical purposes.

Chapter 9 describes the basics of the indicator dilution theory on which all the perfusion modeling is based. In particular, it introduces the perfusion parameters of interest, such as the cerebral blood flow (CBF) and the mean transit time (MTT). Then, it introduces the fundamental concept of response function, R(t), an unknown time-dependent function that describes entirely the perfusion within a voxel, and that is generally estimated via deconvolution. The chapter also presents the challenges related to the estimation of the response function, tackled in the next chapters, such as the time-delay issue, and bolus dispersion itself.

Chapter 10 presents the mathematical description of bolus dispersion, and the proposed phenomenological interpretation of it. It also reports the effects of dispersion on the estimated response function and perfusion parameters, and introduces the proposed dispersion time to characterize them. It finally clarifies the concept that, in the presence of dispersion, deconvolution leads to the estimation of effective, i.e. apparent, response function and perfusion parameters, as opposed to a dispersion-free estimation. These, will be estimated with the framework proposed in the next chapter.

Chapter 11 constitutes the main methodological chapter of the perfusion part. It proposes a description of the effective response function by means of Dispersion-Compliant Bases (DCB): these are proposed in a non-linear and linearized versions. The chapter also reports details about state-of-the-art techniques implemented in this thesis to perform comparisons. After proposing DCB deconvolution, it is also presented a method to recover dispersion-free perfusion parameters. The framework will be validated in the next chapter.

Chapter 12 presents experiments to validate in silico and in vivo the performance of the proposed DCB deconvolution framework, introduced in the previous chapter. It shows that the DCB method has higher precision, compared to stateof-the-art techniques, in the estimation of the effective response function and perfusion parameters. Moreover, it is shown that the use of DCB deconvolution as a pre-processing step improves the accuracy of dispersion-free response function and perfusion parameters via methods that employ a model to describe dispersion.

Chapter 13 proposes a method to infer in vivo on the correct model to adopt to represent dispersion, in order to obtain a correct estimation of dispersion-free response function and perfusion parameters. At the same time, the chapter provides evidence of the importance of performing perfusion deconvolution with the least amount of assumptions about the possible undergoing dispersion process, e.g. without adopting a model of dispersion during the deconvolution-based estimation procedure, which is the case of the proposed DCB approach.

Part IV: "Conclusion"

This part is fundamentally constituted by chapter 14 which reports the summary of the contributions for both part II (diffusion) and part III (perfusion). In addition, some perspectives are described.

Chapter 2

Introduction (en français)

CHAPTER 2. INTRODUCTION (EN FRANC ¸AIS)

Les techniques d'imagerie par résonance magnétique de Diffusion (IRMd) et de Perfusion (IRMp) permettent la détection de divers aspects importants et complémentaires en imagerie cérébrale. Le travail effectué dans cette thèse présente des contributions théoriques et méthodologiques sur les modalités d'IRM basées sur des images pondérées en diffusion, et sur des images de perfusion par injection de produit de contraste. Pour chacune des deux modalités, les contributions de la thèse sont liées au développement de nouvelles méthodes pour améliorer la qualité, le traitement et l'exploitation des signaux acquis.

En IRM de diffusion, la nature complexe du signal est étudiée avec un accent sur l'information de phase. Le signal complexe est ensuite exploité pour corriger le biais induit par le bruit d'acquisition des images, améliorant ainsi l'estimation de certaines métriques structurelles. Après un premier chapitre d'introduction, les trois suivants étudient la nature complexe du signal de diffusion, et son lien avec le phénomène physique de diffusion. Un dernier chapitre est ensuite consacré à l'exploitation du signal complexe dans le cadre d'un processus de débruitage des images DWI par correction de phase.

En IRM de perfusion, le traitement du signal est revisité afin de tenir compte du biais dû à la dispersion du produit de contraste. On montre comment ce phénomène, qui peut empêcher l'estimation correcte des métriques de perfusion, peut aussi donner des informations importantes sur l'état pathologique du tissu cérébral. Après un premier chapitre introductif, où les bases de la théorie de la perfusion sont révisées, l'accent est mis sur les concepts liés à la caractérisation d'un phénomène connu sous le nom de dispersion du bolus, qui polarise les cartes scalaires du flux sanguin cérébral (CBF) et du temps moyen de transit (MTT) qui constituent des résultats importants en IRMp. Cependant, cette information de dispersion est potentiellement liée à la présence de maladies steno-occlusives et peut donc fournir des informations importantes. Pour cette raison, dans les chapitres suivants il est proposé un cadre de traitement de données de perfusion qui se base sur des bases fonctionnelles (DCB) afin d'intégrer cette présence éventuelle de dispersion du bolus et arriver à la caractériser. La méthode est validée sur des images synthétiques et réelles. Le cadre présenté, qui évite la modélisation explicite de la dispersion, se révèle utile en tant qu'étape de prétraitement du signal afin d'améliorer les résultats de techniques basées sur des modèles explicites. Le dernier chapitre utilise le cadre proposé pour inférer in vivo la nature de la dispersion, une procédure originale et efficace qui pourrait ouvrir de nouveaux scénarios potentiels de caractérisation et de quantification de la dispersion à des fins cliniques.

Les contributions apportées dans cette thèse sont présentées dans un cadre théorique et méthodologique validé sur de nombreuses données synthétiques et réelles.

Overview

What are diffusion and perfusion? How do they relate to the nature of the Nuclear Magnetic Resonance (NMR) signal? How can the measured signal be related to physical quantities of interest, such as the diffusion coefficient or the concentration of a tracer? This chapter starts by summarizing the concepts of NMR, necessary to interpret the related images obtained in Magnetic Resonance Imaging (MRI), and introduces the fundamental equations used to describe the diffusion signal acquired with Diffusion-Weighted MRI, and the perfusion signal acquired with Dynamic Susceptibility Contrast MRI.

Keywords NMR; spin-echo; diffusion; perfusion.

Concepts of Nuclear Magnetic Resonance

The physical phenomenon of NMR allows observing the quantum mechanical magnetic properties of the different isotopes of the atomic nucleus. These properties are relevant in medical imaging since, for instance, they allow to differentiate between different tissue types. The fundamental property of the atom nucleus we are concerned with is the spin. Indeed, the isotopes' nuclei that have an odd mass number (the sum of the number of protons and neutrons) possess an intrinsic magnetic angular moment µ or, in other words, a spin. A spin induces a magnetic field. However, within an ensemble of spins the resulting moment is zero because they are randomly oriented in the space, as shown in fig. 3.1a.

Equilibrium

When the ensemble is placed in a powerful external magnetic field, B 0 , this imposes a mechanic moment on each spin according to

d µ dt = γ µ × B 0 (3.1)
where γ is the gyromagnetic ratio that depends on the type of nucleus. For instance, the hydrogen's isotope 1 H, very relevant for studying biological tissues, has γ = 267.513 10 6 rad/sT . The mechanic moment in combination with gravity induces a precession, i.e. a rotation within a cone about the field's direction, with an angle θ ≈ 54 • , and an angular speed ω 0 defined by the Larmor's law

ω 0 = 2π B 0 = γ B 0 (3.2)
where = h/2π with h the Planck's constant. Under the effect of the static field, B 0 , the energy state of a spin is described by two eigen-states, i.e. the directions of angular momentum arising from nuclear spin are quantized: one associated with the alignment of an individual spin towards the direction imposed by the field, i.e. the parallel direction associated to a low energy level E up ; one associated with the opposite alignment, i.e. the anti-parallel direction, associated to a high energy level E down . These eigen-states form a basis for describing any spin-state, the spin being capable of assuming any orientation in the 3D space [START_REF] Hanson | Is quantum mechanics necessary for understanding magnetic resonance?[END_REF]. The difference between the energy-states is

∆E = E down -E up = γ B 0 . (3.3)
At thermal equilibrium, in the presence of B 0 , an observable ensemble of spins will have a global orientations distribution that can be described as if the parallel spins slightly outnumbered the anti-parallel ones according to

N up N down = exp - ∆E k B T (3.4)
where k B = 1.38×10 -23 J/K is the Boltzmann constant and T the absolute temperature in Kelvin. This relationships shows that NMR is a low sensitivity technique. Indeed, at B 0 = 1T and T = 300 • K, and considering that ∆E

k B T N up -N down N up + N down ≈ γ B 0 2k B T ≈ 3 × 10 -6 (3.5)
which means that about 3 spins per million contribute to obtain a positive net magnetization

M = N n=1 µ n (3.6)
where N = N up + N down . At equilibrium the magnetization vector M is aligned with B 0 , that is along the z axis of the spatial frame of reference xyz. The vector can thus be decomposed into a parallel longitudinal component, M z = 0, and a perpendicular transverse one, M xy = 0. This situation is illustrated in fig. 3.1b. Spins, however, have random phase. Since spins are precessing at ω 0 about B 0 , it is useful to visualize the system in a frame of reference, xyz ω 0 , that rotates at the same speed about the z axis, as shown in the second line of fig. 3.1a.

Excitation

The system at equilibrium can be excited by applying a second magnetic field, B 1 , of lower intensity and pulsating at radio frequency ω RF = ω 0 to obtain nuclear magnetic resonance. The field is perpendicular to B 0 and in xyz ω 0 can be visualized as oriented towards the negative x axis (fig. 3.1c). Depending on the duration of the application of the RF the net magnetization M is flipped with a nutation angle α. For instance, to perform a flip α = π/2, the 90 • RF pulse has duration

τ 90 • = π/2 γ B 1 .
(3.7)

B(t) = (B x (t), B y (t), B 0 + ∆B z (t)
) is described by the phenomenological equations derived by [START_REF] Bloch | Nuclear induction[END_REF]

dM x (t) dt = γ M (t) × B(t) x - M x (t) T 2 dM y (t) dt = γ M (t) × B(t) y - M y (t) T 2 dM z (t) dt = γ M (t) × B(t) z - M z (t) -M eq z T 1 (3.8)
where M eq z is the value of the z component at thermal equilibrium.

Relaxation

In NMR the focus is on studying the restoration of the system's equilibrium after it has been perturbed, particularly the reappearance of a non-zero component M z , The two types of relaxation happen, in a counterintuitive manner, according on two different physical phenomena. For this reason, they are described by exponential laws with different time rates, T 1 and T 2 [START_REF] Bloch | Nuclear induction[END_REF]. In the case of a 90 • RF pulse

M z = M eq z 1 -exp - t T 1 M xy = M eq z exp - t T 2 (3.9)
where T 1 corresponds to the time in which longitudinal component M z regain a value which is the 37% (1/e) lower of its equilibrium value M eq z , and T 2 corresponds to the time in which the transverse component is reduced to 37% of it. The relaxations are shown in fig. 3.3. The longitudinal relaxation occurs according to the energy exchange between spins with the surrounding molecules, the "lattice", therefore T 1 is either called the longitudinal or spin-lattice relaxation time. Instead, the transverse relaxation is mainly due local interactions between spins, without necessarily involving energy exchange, therefore T 2 is either called the transverse or spin-spin relaxation time.

The measurement of the transverse relaxation time T 2 is not trivial. Indeed, within a real ensemble of spins there are minor differences in the chemical environment, such that the pulsations of the individual spins distribute around the ideal ω 0 . Over time, this leads to a loss of phase coherence among the spins, with a consequent signal loss measured by the Free Induction Decay (FID). The situation is depicted in fig. 3.2 The measured decay of M xy is then always faster than the one depending on T 2 . In fact, the FID decays with a time rate 1/T * 2 that is always larger than 1/T 2 according to

1 T * 2 = 1 T 2 + γ∆ B 0 (3.10)
where the second term includes the contribution to the local variation of the static magnetic field intensity.

Spin-Echo

The variation of the filed intensity in eq. (3.10) largely depends on the location of the spins, and can be compensated unless spins are subject to motion, such as diffusion. In fact, T 2 can be measured with the Spin-Echo (SE) sequence described by [START_REF] Erwin L Hahn | Spin echoes[END_REF] and illustrated in fig. 3.4. After the 90 • RF pulse the magnetization vector M lies in the xy plane, i.e. M xy = 0 and M z = 0. Bearing in mind that M is the sum of the single spins' moments, according to eq. (3.6), it is possible to see that if spins rotate at the same speed ω 0 about the z axis, then the magnitude of the magnetization is preserved during the rotation, neglecting T 2 effects for the sake of clarity. However, the field inhomogeneity causes the ensemble of spins to follow a distribution of speeds around the ideal ω 0 . Therefore, some spins rotate faster and some others slower as shown in fig. 3.2. This loss of phase coherence causes the spins' moments to add destructively, thus decreasing the magnitude of the net magnetization over time, much beyond the desired T 2 contribution. Hahn had the idea of adding an additional RF pulse of 180 4). This has the effect of flipping the spins of 180 • in the xy plane such that now, the moments of the spins that were faster than ω 0 are still faster but are located behind (in the rotating frame) the net magnetization vector, and conversely the moments of the spins that were slower are ahead. The idea is that there is a time instant in which the faster spins located behind will catch up with the slower ones, refocusing the spins and reestablishing the phase coherence. This means that the magnitude of the net magnetization vector, measured during time, first decreases due to destructive summation, i.e. the measured FID, and after will increase due to refocusing, thus reaching a maximum. The time of the maximum is the echo-time TE, that is the time at which the echo of the signal occurs. By repeating the sequence with increasing distances between the RF pulses, that is half the TE, the peaks of the echoes will reveal the successive samples of the actual T 2 decay.

Diffusion in NMR: a Dephasing Force

Water diffusion plays a central role in the generation of the NMR signals on which are based the MRI modalities discussed in this thesis: Diffusion-Weighted (DW-MRI) and Dynamic Susceptibility Contrast (DSC-MRI). In the first case, diffusion of water spins is the object of interest, the measured contrast, in the second it is a major phenomenon determining the contrast. In both cases diffusion is a natural process that in NMR is an established cause of spin dephasing, that for both of the techniques means signal loss. The situation is depicted in fig. 3.5. In the [START_REF] Edward | Spin diffusion measurements: spin echoes in the presence of a time-dependent field gradient[END_REF].

left, a spin displaces -due to diffusion -within a constant magnetic filed along the z axis with intensity B 0 that imposes a precession speed ω 0 according to eq. (3.2). While moving, the spin does not experience changes in the precession speed, hence it does not accumulate phase. On the right, the magnetic field varies along z according to a gradient. The spin moving within the gradient experiences different velocities ω (increasing from left to right in the figure), which causes the spin to accumulate phase based on the amount of displacement. Different spins within an ensemble moves differently and accumulate different phase contributions.

After some time the spins will be out of phase, thus reducing the total amount of transverse magnetization that can be measured. This simplified concept is useful to understand both the principles of DW-MRI and DSC-MRI.

Measuring Diffusion with DW-MRI

As previously described in section 3.3, diffusion leads to some peculiar effects in NMR. After [START_REF] Erwin L Hahn | Spin echoes[END_REF] noticing these effects, Carr and Purcell [1954] formalized them and Torrey [1956] modifies the phenomenological eq. (3.8) of Bloch by adding the Fick's diffusion equation to it

dM x (t) dt = γ M (t) × B(t) x - M x (t) T 2 + ∇ • D∇(M x -M eq x ) dM y (t) dt = γ M (t) × B(t) y - M y (t) T 2 + ∇ • D∇(M y -M eq y ) dM z (t) dt = γ M (t) × B(t) z - M z (t) -M eq z T 1 + ∇ • D∇(M z -M eq z ) (3.11)

MEASURING DIFFUSION WITH DW-MRI

where D is the diffusion coefficient, and where isotropic diffusion is considered.

Based on these works [START_REF] Edward | Spin diffusion measurements: spin echoes in the presence of a time-dependent field gradient[END_REF] proposed to modify the SE sequence by adding two identical magnetic gradients around the 180 • RF pulse, in order to encode the diffusion of the spins. This sequence, that established the field of diffusion NMR, is the pulsed gradient spin echo (PGSE) which is shown in fig. 3.6. The idea is basically to exploit the phenomenon illustrated in fig. 3.5.

In the PGSE sequence, each gradient pulse G has amplitude g = | G|, direction g, duration δ, and they are separated by a time ∆. However, the two pulses effectively generate gradients of opposite sings due to the effect of the 180 • RF pulse in between them. Considering the total magnetic field experienced by a spin in position r, B( r, t) = (0, 0, B 0 + r • G(t)), Torrey shows how to apply the new eq. (3.11) to find the spin-echo amplitudes M x (t) and M y (t) in the presence of an inhomogeneous field. He considers the complex xy magnetization at the spin location r, with pulsation ω 0 and accounting for T 2 relaxation

M x ( r, t) + jM y ( r, t) = ψ( r, t) exp -jω 0 t - t T 2 (3.12)
and derives that the amplitude ψ(t) follows

∂ψ( r, t) ∂t = -jγ( r • G)ψ + D∇ 2 ψ (3.13)
from which it can be seen that the diffusion contribution, the real part of the equation, is responsible for the amplitude attenuation (the imaginary part gives a phase contribution instead). Following the formalism of Torrey, the effects of diffusion can be separated into a function A(t), and the solution of eq. (3.13) in case of the PGSE, developed by [START_REF] Edward | Spin diffusion measurements: spin echoes in the presence of a time-dependent field gradient[END_REF], is simplified as

ψ( r, t) = A(t) exp -jγ r • F (t) (3.14)
where, rewriting eq. 3 found in the same reference,

F (t) = t 0 G(t )dt -2Θ t - T E 2 T E 2 0 G(t )dt (3.15)
with Θ(t) the Heaviside step function with convention Θ(t) = 1 ∀t > 0. Assuming a solution to eq. (3.13) identical to that in eq. (3.14), substitution leads to a first order differential equation describing the evolution over time of the diffusion contribution A(t) to the amplitude attenuation of ψ according to

dA(t) dt = -γ 2 D| F (t)| 2 A(t) (3.16)
that when integrated between the time of the 90 • RF pulse, t = 0, and the echo, T E, leads to [START_REF] Edward | Spin diffusion measurements: spin echoes in the presence of a time-dependent field gradient[END_REF] ln

A(T E) A(0) = -γ 2 δ 2 g 2 ∆ - δ 3 D (3.17)
which describes the attenuation of M xy only due to diffusion: this is the reduction of the echo amplitude A(T E) compared to that right after the 90 • RF pulse, A(0). Transposing this concept when taking into account T 2 relaxation, defining S the signal at T E acquired in the presence of diffusion gradients, and S 0 the signal acquired without diffusion gradients, leads to

S(b) = S 0 e -bD (3.18) where b = γ 2 δ 2 g 2 ∆ - δ 3 (3.19)
is the b-value describing the global diffusion weighting. Equation 3.18 reveals that, in the case of isotropic free diffusion, the signal decays exponentially as a function of the diffusion weighting. [START_REF] Stejskal | Use of spin echoes in a pulsed magnetic-field gradient to study anisotropic, restricted diffusion and flow[END_REF] proposed an extension of the theory to anisotropic media by introducing the concept of diffusion tensor D. This theory was formalized later on by Basser et al. [1994a,b] who provided a an efficient way of estimating the tensor allowing one to represent the signal as

S(b • g) = S 0 e -b g T D g (3.20)
which forms the basis for the technique known as Diffusion Tensor Imaging (DTI).

Perfusion: a functional indicator for the brain

The term perfusion indicates the physiological phenomenon of the passage of blood in the parenchyma, the latter being the functional part of a tissue. The amount of perfusion is related to the metabolism of the parenchyma and to its grade of activity. This thesis focuses on the brain tissue parenchyma, which is the relevant domain of application of DSC-MRI. Perfusion can be schematized as in fig. 3.7, where the blood coming from an artery goes through the capillary bed, from where the red cells can exchange oxygen for carbon dioxide. Then, the red cells are transported until a venous blood vessel until they reach the lungs, where the respiration process takes place. One of the main indicators related to perfusion is the blood flow. Different tissue types, for instance within the brain, have different normal levels of blood flow but these might change in the presence of pathology.

The measure of perfusion is typically performed by means of tracers. A tracer is a substance that is not metabolized by the tissue and that can be detected (and measured) separately from the blood, but at the same time, that is transported together with it. In this way it is possible to measure the behavior of blood perfusion by just measuring the one of the tracer. In MRI the perfusion can be measured via different techniques. The main division among these regards the nature of the tracer which can be endogenous or exogenous. In the first category falls the Arterial Spin Labeling (ASL), which involves the magnetic labeling of the water particles in the arterial blood that will be used as tracer. In the second category there are the Dynamic Susceptibility Contrast MRI (DSC-MRI) and the Dynamic Contrast Enhanced MRI (DCE-MRI). These last two techniques mainly differ in the physiological assumption regarding the integrity of the blood-brain barrier (BBB), i.e. the epithelial barrier occurring along capillaries in the brain that separates the circulating blood from the brain extracellular fluid. Indeed, the tracer is confined within the capillaries by the BBB. However, if this is damaged then the tracer can leak out. Whether DCE explicitly accounts for leakage, allowing for permeability modeling (fig. 3.7), DSC assumes BBB integrity which is an important notion for the understanding of the origin of the signal.

Measuring Perfusion with DSC-MRI

Dynamic Susceptibility Contrast MRI (DSC-MRI) exploits the decrease of signal intensity, i.e. the transverse magnetization M xy , in the presence of a paramagnetic (contrast) agent (PA). To understand the nature of the measured signal it is therefore important to understand the effects of paramagnetism on NMR relaxation. One of the first complete dissertations on the subject is that written by [START_REF] Bruce R Rosen | Perfusion imaging with nmr contrast agents[END_REF].

Effect of Paramagnetic Agents on NMR Relaxation

Paramagnetism is the phenomenon for which some materials are attracted by external magnetic fields, in contrast with diamagnetic materials which are repelled. The most studied PA for in vivo applications is the gadolinium, Gd, coupled with the diethylenetriamine penta-acetic acid in the form of Gd-(DTPA) -2 . Microscopically, the atoms of PAs have a permanent magnetic moment due to the presence of at least one unpaired electron in the external orbital. This unpaired paramagnetic electron has a net spin which is not canceled out by any other electron in the orbital, thus causing the entire atom to have a net magnetic moment.

The contrast agents used in NMR can affect the signal mainly according to two different mechanisms: relaxation and susceptibility. In the first case, the PA enhances the relaxation due to dipole-dipole interactions between the free electron spin and the surrounding proton spins. This causes the relaxation enhancement of both T 1 and T 2 but, since the rate R 1 = 1/T 1 is typically smaller than R 2 = 1/T 2 the enhancement due to dipole-dipole interactions is more evident, in percentage, with T 1 contrast [START_REF] Bruce R Rosen | Perfusion imaging with nmr contrast agents[END_REF]. Hence, the variation in relaxivity, ∆R 1 , depends on the PA concentration in the tissue, [T issue], according to a linear relationship

∆R 1 = k 1 [T issue] (3.21)
where k 1 is a tissue and field strength specific constant. However, in the brain tissue, the assumptions of the expected mono-exponential longitudinal relaxation do not comply with the nature of the T 1 relaxivity enhancement due to the presence of an intact BBB [START_REF] Bruce R Rosen | Perfusion imaging with nmr contrast agents[END_REF]. Indeed, relaxivity effects dominate within aqueous materials when the PA is distributed homogeneously. However, when the BBB is intact the PA is compartmentalized within the capillaries thus rendering an inhomogeneous distribution [START_REF] Villringer | Dynamic imaging with lanthanide chelates in normal brain: contrast due to magnetic susceptibility effects[END_REF]. For this reason, relaxivity is thus exploited outside the brain where there is no BBB or, within the brain, by means of DCE-MRI which explicitly models a non-intact BBB.

The second mechanism of action of PAs is susceptibility, which mainly influences T * 2 and T 2 relaxations causing a signal drop, as opposed to the enhancement due to the relaxation mechanism. Also in this case [START_REF] Bruce R Rosen | Perfusion imaging with nmr contrast agents[END_REF]]

∆R 2 = k 2 [T issue] (3.22)
where k 2 is a different tissue and field strength specific constant. Similarly to what described for nuclear spins in section 3.1, the unpaired electrons of elements such as Gd have moments that respond to an externally applied magnetic field. When the unpaired spins are present in sufficient concentration, their orientation distribution induces a net magnetization that is aligned with the applied field, i.e. paramagnetic materials display a positive magnetic susceptibility that is the proportionality constant between an applied magnetic field and the resultant magnetization established within the material. In the presence of an intact BBB, as previously mentioned, the PA is compartmentalized within the capillaries leading to local differences of susceptibility within the tissue. These differences imply local differences of the magnetic field, affecting the local resonance pulsation ω of the water proton spins according to T * 2 effects, similarly to what pictured in the right image of fig. 3.2. However, these static effects are compensated when using a SE sequence. On the other hand, there are also dynamic effects caused by the diffusion of water protons within the magnetic field gradient exerted by the heterogeneous distribution of PA, or by the variation of the gradient itself (similarly to what illustrated in the right image of fig. 3.2). On top of these mesoscopic effects there also are microscopic and macroscopic effects influencing the relaxivities. Whether the macroscopic effects are due to MRI issues, such as imperfect shimming, the microscopic ones, due to dipole-dipole interactions, only account fr about the 16% of the total signal drop due to transverse relaxation in SE 1 sequences [START_REF] Valerij | Transverse relaxation effect of MRI contrast agents: A crucial issue for quantitative measurements of cerebral perfusion[END_REF]. The transverse relaxation is then mainly due to diffusion.

Diffusion: a dephasing force in perfusion

When the PA is compartmentalized, a magnetic field shift is present in the intravascular tissue whereas field gradients are generated around the vessel's lumen [START_REF] Villringer | Dynamic imaging with lanthanide chelates in normal brain: contrast due to magnetic susceptibility effects[END_REF]. Water protons diffusing within the gradients generated locally around the blood vessels (containing the PA) accumulate phase, as shown in fig. 3.5 which globally causes a loss of proton spins coherence within the voxel, thus resulting in a reduced transverse magnetization and signal loss. This scenario is depicted in fig. 3.8. Simulations performed by [START_REF] Villringer | Dynamic imaging with lanthanide chelates in normal brain: contrast due to magnetic susceptibility effects[END_REF] show that considering a tissue volume with capillaries, almost the 60% of the protons within the intercapillary space are exposed to significant gradients. Furthermore, the amount of diffusion influencing the signal is directly related to the time the protons can move before the signal acquisition, therefore the signal drop, as will be more clear later, directly depends on the echo-time TE [START_REF] Bruce R Rosen | Perfusion imaging with nmr contrast agents[END_REF]. Indeed, the contrast of MRI images acquired via SE increases with TE. The contrast also depends on the size of the vessels. Models of the role of diffusion in the signal generation are also reported by [START_REF] Dmitriy | Theory of nmr signal behavior in magnetically inhomogeneous tissues: the static dephasing regime[END_REF], [START_REF] Jerrold L Boxerman | Mr contrast due to intravascular magnetic susceptibility perturbations[END_REF], and are extended by [START_REF] Kiselev | Analytical model of susceptibility-induced MR signal dephasing: effect of diffusion in a microvascular network[END_REF] who take into account also the smaller contribution from the intravascular compartment. where κ = k 2 . The linear relationship in eq. (3.22) does not hold in many cases. Indeed, it can be superlinear or sublinear depending on the type of PA used [START_REF] Weisskoff | Microscopic susceptibility variation and transverse relaxation: theory and experiment[END_REF], therefore correction is often performed. Following an injection of tracer into the system, its concentration measured within a voxel changes during time. Consequently, the susceptibility in the tissue and the gradients of the magnetic field to which water proton are exposed change dynamically, during time. Dynamic susceptibility implies a varying relaxivity during time. By repeating the acquisition for consecutive time instants, while keeping TE fixed, it is possible to measure the dynamic of the signal that is described by

S(t) = S 0 e -κC(t)T E (3.25)
where C(t) is the tracer's concentration during time. From eq. (3.25) it is straightforward to determine the concentration from the measured signal as

C(t) = - 1 κT E ln S(t) S 0 . (3.26)
Figure 3.9 shows, in principles, the signal behavior as the concentration of tracer within a voxels varies.

Conclusion

The signal generation of Diffusion and Perfusion NMR is linked by the action of the diffusion of water particles. From a mathematical point of view it is straightforward to notice the parallelism between eq. (3.18) and eq. (3.23), here reported

   S(b) = S 0 e -bD S(T E) = S 0 e -T E∆R 2
the first being exploited in DW-MRI to infer the tissue structure, and the second to retrieve functional information in DSC-MRI. Both of the techniques involve processing at least four-dimensional (4D) data, that means one additional dimension to the three spatial ones. In diffusion it is typical to acquire images -the diffusion-weighted images (DWIs) -along multiple directions in the space, and/or with multiple b-values. In perfusion, images are acquired in successive time points, i.e. perfusion-weighted images (PWIs), to allow characterizing the tissue concentration of the PA during time. However, since the aims of the techniques are orthogonal, the underlying processing are quite different. The rest of the thesis will be divided in two parts, one for diffusion and one for perfusion, where contributions will range from the physical interpretation of the signal and signal denoising in the case of diffusion, and modeling in the case of perfusion.

Overview

In the most general terms, diffusion of water spins is described via the probability density (pdf) function of the displacements that spins undergo during the measuring time interval. This chapter describes the link between the measured diffusion signal, described in chapter 3, and the description of the diffusion process as a pdf. The thesis builds on top of this link to explore different aspects of its nature, specifically involving the complex nature of the diffusion signal. On one side it provides a theoretical exploration on the diffusion process that generates a complex signal. On the other side, the complex signal is exploited for the sake of denoising purposes.

The chapter also describes the diffusion tensor (DT) and mean apparent propagator (MAP) signal representations that will be used in chapter 8.

Keywords microstructure; ensemble average propagator; asymmetry; EAP.

Diffusion in tissue microstructure

The characteristics of the diffusion process in the brain depend on how water spins are restricted by the surrounding tissue environment, that is the action of those phenomena and physical boundaries that prevent the free movement of spins themselves. To understand this relationship, it is convenient to focus on the movement of an individual spin and then expand this concept to the movement of an ensemble of spins. The movement of the ensemble is better understood in terms of the probability density function of the spins' displacements from a given position, due to diffusion. In fact, water spins undergo random motion caused by collisions with other neighboring ones. While every diffusing spin follows its own deterministic path, when considering an ensemble only the probability that a specific displacement occurs can be reasonably characterized. This pdf describes the probability density P (r, τ ) that a particle will travel a distance r ∈ R 3 during diffusion time τ ∈ R + , and is often referred to as the ensemble average propagator (EAP) [START_REF] Kärger | The propagator representation of molecular transport in microporous crystallites[END_REF]. [START_REF] Einstein | Investigations on the Theory of the Brownian Movement[END_REF] showed that when diffusion occurs in an unrestricted environment, i.e. a free diffusion, the EAP is described by a Gaussian distribution as

P (r, τ ) = 1 √ 4πDτ e -r 2 4Dτ (4.1)
where D is the diffusion coefficient. Equation (4.1) shows that the probability of larger displacements increases with D or τ . isocontour for τ 3 . The shape of the isocontour depends on the structure of the surrounding tissue. For instance, in the left image free water particles are unobstructed and can travel further with isotropic probability. Gray matter tissue, in the second image, is generally unorganized and hinders diffusion in all directions, also producing isotropic, but smaller contours compared to free water. White matter, instead, is mainly organized in bundles of axons. This is illustrated in the third image, where axons, illustrated as gray lines, are mostly aligned in one main direction.

Particle movement is obstructed perpendicularly to this direction and is relatively free along it, causing anisotropic isocontours. Finally, in areas where two bundles cross, there is a mix between the isocontours of each bundle. Note that, in the right two white matter tissues, isocontours are intentionally drawn more isotropic for τ 1 than for τ 2 and τ 3 . Indeed, for shorter τ particles have not had much time to interact with the surrounding tissue, resulting in similar probability that a particle travels in any direction. The isocontours of lim τ →0 P (r, τ ) will therefore always be isotropic. intuitively, for longer τ particles potentially render more anisotropic profiles.

Images in fig. 4.1 illustrate very simplified scenarios. In fact, gray matter is distinguished from white matter in that the former contains numerous cell bodies and relatively few myelinated axons, whereas white matter is composed mainly of long-range myelinated axons and contains relatively few cell bodies. The brain also contains glial cells of various kinds having the main function of supporting the activity of neurons. Diffusion occur within and between all the various types of structures present in the tissue.

In the following section the mathematical object EAP is defined, arriving to the establishment of the main relationship between the acquired signal and the EAP itself.

Ensemble average propagator

The EAP is a mathematical object that describes the probability density function of spins' displacements due to diffusion in the three-dimensional space, and it is related to the measured DW-MRI signal -under some simplifications and assumptionsvia a Fourier relationship. Although grounded on the mathematics of the PGSE described in section 3.4, the derivation of this relationship is based on a probabilistic "random walk" approach rather than a phenomenological one.

The diffusion encoding can be interpreted as the result of a differential mechanism. Indeed, in the PGSE sequence shown in fig. 3.6 spins are "labeled", before and after the 180 • RF pulse, by the two gradient pulses of G(t) = G(t). In the hypothesis of gradient pulses applied along one constant direction let us consider a spin that displaces during the time between the two pulses. For convenience, let us consider that the spin is located in a position r 0 during the first gradient, and, after moving during the time ∆ between the pulses, it is found in position r 1 during the second. The influence of each pulse can be quantified in the dephasings

φ 0 = γ δ 0 G(t)r 0 dt = γgδr 0 φ 1 = γ ∆+δ ∆ G(t)r 1 dt = γgδr 1 (4.2)
which holds if G(t) = g and the spins are static during the application of the gradient pulses. This second condition is assumed by considering δ ∆ which is known as the Narrow Pulse Approximation (NPA). Therefore, the displacement of spins during the time ∆ causes a net dephasing

ϕ = φ 1 -φ 0 = γδG • (r 1 -r 0 ) (4.3)
where the minus sign is due to the effect of the 180 • RF pulse. Equation (eq. ( 4.3)) expresses an important concept: if, between the two pulses, the spin has moved w.r.t. its initial position, r 0 , with net displacement, then it is subject to a nonnegative phase accumulation ϕ = 0; if, on the other side, the final position corresponds to the initial one, r 1 = r 0 , then φ = 0. The complex signal contribution of a dephasing, T 2 effects apart, generated for an individual spin with phase ϕ is [START_REF] Paul | Principles of nuclear magnetic resonance microscopy[END_REF] e jϕ = e jγδG • (r 1 -r 0 ) (4.4) which clearly shows the absence of contribution in the case there is no phase accumulation. The signal attenuation E(G, ∆) = S(G, ∆)/S 0 at the echo time acquired within a voxel is the average of the contributions of the ensemble of spins therein.

If ρ(r 0 , 0) is the probability of finding a spin initially at r 0 , and P (r 1 |r 0 , ∆) the probability of a spin to displace to r 1 given the initial position r 0 in a time ∆, then the signal attenuation is E(G, ∆) = ρ(r 0 , 0) P (r 1 |r 0 , ∆)e jγδG • (r 1 -r 0 ) dr 1 dr 0 (4.5) which was initially shown by [START_REF] Stejskal | Use of spin echoes in a pulsed magnetic-field gradient to study anisotropic, restricted diffusion and flow[END_REF]. Since all the terms but the exponential are probability density functions, in the case either there is no diffusion gradient G or no displacement, r 1 = r 0 , then it is straightforward to see that the attenuation is equal to the unity. Particularly, in the first case it is

E(0, t) = 1 (4.6)
where t is used instead of ∆ because the latter has no meaning when there are no gradients applied. Henceforth, there will follow some assumptions and simplifications that are useful to arrive at the Fourier relationship mentioned above. Indeed, Callaghan [1993] noticed that through some assumptions it is possible to link directly the signal to the average displacement probability via an inverse Fourier transform. To do so, he groups some quantities present in the exponent of the exponential term in a factor q expressed as q = γδG 2π (4.7)

such that now eq. ( 4.5) can be rewritten as E(q, ∆) = ρ(r 0 , 0) P (r 1 |r 0 , ∆)e j2πq • (r 1 -r 0 ) dr 1 dr 0 (4.8) where the Fourier exponential term is now visible. The signal acquired with the PGSE depends only on the displacement r = r 1 -r 0 , therefore the probability can be written as P (r, t) which is the propagator for an individual spin: the probability that a spin displaces of a distance r in a time t. Then, by manipulating the integrals it is possible to make the ensemble average propagator (EAP), P (r, t), to emerge, which is expressed as P (r, t) = ρ(r 0 , 0)P (r, t) dr 0 (4.9) that combined with eq. ( 4.8) gives E(q, t) = P (r, t)e j2πq • r dr (4.10) that establishes de facto an inverse Fourier relationship between the EAP, P (r, • ), and the signal attenuation, E(q, • ). Now, r denotes the displacement vector, normally expressed in micrometers (µm), hence q is a spatial frequency vector. Equation (4.10) describes the main result of what is known as the q-space formalism. This relationship is at the heart of the more theoretical contribution of this thesis. However, when assuming antipodal symmetry the signal attenuation, E(q) = E(-q), this translates to E(q, t) = P (r, t)e -j2πq • r dr (4.11) that is the relationship generally used in diffusion MRI. Generally, higher q-value, q , and diffusion time, τ , lead to a larger signal attenuation. However, they influence differently the diffusion-weighting, i.e. the b-value expressed in eq. (3.19) which can be rewritten as b = 4π 2 q 2 τ (4.12)

where q = q and b = b . In fact, q is a spatial frequency, and by increasing its value it is possible to achieve a higher spatial resolution, of P (r, τ ) in the displacement space described by r = r . On the other hand, the value of the diffusion time expresses the state of the diffusion that one wants to measure, similarly to the exposure time of a camera: the longer one waits the higher the amount of light the film is exposed to, and the lighter the resulting picture becomes. This concept, of taking a picture of the tissue microstructure at a specific time, is developed in the next section.

A spatial picture at a specific time

The diffusion time τ expresses the time interval during which spins are allowed to diffuse before measurement. A longer diffusion time allows the spins to move a longer distance causing, in the absence of restrictions to the diffusion process, a larger net phase shift, i.e. a stronger attenuation of the signal. Therefore, expressing the diffusion-weighting in terms of q and τ can provide useful insights on the signal nature. In the absence of restrictions to the diffusion process, eq. ( 4.11) has a closed form. This is obtained by substituting eq. (4.12) in eq. (3.18) as

E(q, τ ) = S(q, τ ) S 0 = e -4π 2 q 2 τ D (4.13)
which expresses a Gaussian attenuation profile as function of q, and corresponds, via the Fourier relationship, to the EAP expressed before in eq. ( 4.1). However, eqs. (3.18) and (4.13) are valid when the diffusion process can be considered unrestricted, e.g. when the movement of spins is not obstructed by the presence of a barrier. In the case of restricted diffusion, for instance when the signal is measured along a direction perpendicular to a barrier, these equations are no longer valid. directions parallel and perpendicular to the restriction. The figure also illustrates the signal attenuation in the case of parallel, unrestricted diffusion (a,b), and in the case of perpendicular, restricted diffusion (c,d ). The unrestricted attenuations are obtained with eqs. (3.18) and (4.13), whereas the restricted ones are simulated as the diffusion signal attenuation generated within an ensemble of cylinders along the direction perpendicular to the cylinders' axes [START_REF] Paul T Callaghan | Pulsed-gradient spin-echo nmr for planar, cylindrical, and spherical pores under conditions of wall relaxation[END_REF]. Moreover, the curves in a and c are reported as function of the b-value with diffusion time τ 2 , whereas the curves in b and d are functions of the q-value and are reported for increasing diffusion times τ 1 < τ 2 < τ 3 . In the graphs, points of each curve corresponding to b = 1000 s/mm 2 are highlighted with a dot.

Along the unrestricted direction (a,b) the attenuation values, indicated by dots, are lower than the corresponding ones along the restricted direction (c,d ). Indeed, when diffusion is restricted by the presence of the barrier, the spins are subject to a smaller net displacement and the signal attenuates less.

The choice of q and τ to obtain a certain diffusion-weighting, i.e. a specific b-value, assumes different relevance in terms of signal attenuation depending on whether diffusion is restricted or not. In the absence of restrictions (a,b), an increase of q-value or diffusion time always attenuates the signal, and points with different q and τ , but with same b-value, render the same amount of attenuation (b). However, when diffusion is restricted (c,d ), an increase of the diffusion time τ implies letting the spins diffuse a longer distance with the consequence of experiencing more restriction. In this case, the Gaussian attenuation expressed by eq. (4.13) is not longer valid. Indeed, different combinations of q and τ render different non-Gaussian profiles of signal attenuation, and points with same b-value -the dots of fig. 4.2d -correspond to different attenuations.

Till now, only the properties of the diffusion signal and EAP have been described. In the next section signal-driven representations used in this thesis will be covered, as long with the scalar parameters that characterize the diffusion process.

Signal-driven representations and scalar metrics

This section describes three signal representations that will be used in the course of the thesis for different reasons. In the final part it also reports scalar metrics, based on these representations, that are normally computed as signal attenuation/EAP features. Although many representation, or better said, models of the diffusion signal are available in literature, for the scope of this thesis only non-parametric representations are adopted. In fact, one of the purposes of the thesis is diffusion signal denoising and a non-parametric representation, such as the MAP-MRI one presented below, seems to be a suitable choice.

Diffusion Tensor Imaging

The idea behind diffusion tensor imaging (DTI) is that the EAP can be represented as a Gaussian function, like eq. (4.1), in any spatial direction. This leads to the the previously shown eq. (3.20), here rewritten with a more appropriate notation S(bg) = S 0 e -bg T Dg (4.14) where g = (g x , g y , g z ) is the gradient orientation, and D is the diffusion tensor

D =     D xx D xy D xz D yx D yy D yz D zx D zy D zz     = λ 1 v 1 v T 1 + λ 2 v 2 v T 2 + λ 3 v 3 v T 3 (4.15)
where D ij is the apparent diffusivity along the ij direction in a Cartesian reference frame, λ i are the eigenvalues with λ i ≥ λ 2 ≥ λ 3 and v i the corresponding eigenvectors. The estimation procedure proposed by Basser et al. [1994a,b] considers diffusion as antipodically symmetric. Then only the following vector has to be determined

d = lnS(0), D xx , D xy , D xz , D yy , D yz , D zz T (4.16)
where the first entry is there to include the signal measured without diffusionweighting. This vector can be estimated via a collection of at least 6 non-collinear samples, i.e. DWIs, plus S(0) with a total amount of N s ≥ 7 samples with the following objective function

f LLS (d) = 1 2 Ns i=1   lnS(b i ) - 7 j=1 B ij d j   2 (4.17) where B = B 1 , B 2 , • • • , B Ns T (4.18)
and

B i = 1, -b i g x i g x i , -b i 2g x i g y i , -b i 2g x i g z i , -b i g y i g y i , -b i 2g y i g z i , -b i g z i g z i (4.19)
the elements of the estimated d to be used to create the symmetric tensor D. However, in this thesis, for reasons that will be explained after, DTI will be estimated in a non-linear fashion constraining the signals samples to be positive and according to the objective function

f N LS (d) = 1 2 Ns i=1   S(b i ) -exp   7 j=1 B ij d j     2 (4.20)
where also in this case the elements of d are used to create the symmetric tensor.

SHORE

In this thesis it will also be used a representation that relaxes the Gaussian assumption along one gradient direction, making the estimation of a non-Gaussian profile possible. This method is called Simple Harmonic Oscillator based Reconstruction and Estimation (SHORE) [ Özarslan et al., 2008a], which is suitable for the q-space formalism. Within this framework, the signal is represented as the linear combination of orthogonal basis functions, result of the multiplication between an exponential and an Hermite polynomial

φ n (u, q) = i -n u √ π 2 n-1 n! e -2π 2 q 2 u 2 H n (2πqu) (4.21)
where n is the order of the basis, H n (x) is the nth-order Hermite polynomial and u is the characteristic data dependent length or scaling factor to be determined, for instance, by fitting the signal to a Gaussian decay according to E(q) = exp(-2π 2 q 2 u 2 ). The formulation of the basis functions in eq. ( 4.21) differs from that given by Özarslan et al. [2008a] with the introduction of the normalizing factor 2u √ π, which renders the bases orthonormal. The bases are well suited for representing the signal in the complex domain: the even order basis functions are real valued and evenly symmetric whereas the odd order basis functions are imaginary and show odd symmetry, which is precisely the case of the real and the imaginary parts of the diffusion signal. The complex signal attenuation, E c (q), can thus be represented with SHORE as

E c (q) = N -1 n=0 c n φ n (u, q) (4.22)
where, after an estimation of the scaling factor u, the coefficients {c n } can be estimated via linear least squares. N is the maximum allowed order and its choice directly affects the signal reconstruction: in the case of noisy data, a high order will potentially cause the reconstruction to follow the noise, whereas a low order will inherently enforce a smoothing effect. For a given order N the signal reconstruction in the complex domain is performed by considering the even and odd coefficients for the real and imaginary parts respectively. After knowing the coefficients, it is possible to reconstruct the magnitude attenuation from the complex one as

M c = |E c (q)| = {E c (q)} 2 + {E c (q)} 2 . (4.23)
However, when only the magnitude samples are available, the magnitude signal reconstruction is

M = E(q) = |E c (q)| = N -1 n∈even c n φ n (u, q) (4.24)
where only the even coefficients are taken into account. However, in the course of this thesis, it is proposed to do a maximum likelihood estimation (MLE) of the coefficients in the case of a magnitude signal reconstruction. This helps taking into account the noise distribution of the signal samples that is known not to be Gaussian. The estimation is presented in the next subsection, where the considered noise distribution is Rician.

Maximum likelihood estimation

To better take into account the Rician distribution of the noisy magnitude signal samples, the estimation of the coefficients can be performed with the Maximum Likelihood Estimator (MLE). In the case of Rician noise, the MLE has been introduced in [START_REF] Sijbers | Maximum-likelihood estimation of rician distribution parameters[END_REF]. Normally several noisy realizations of each signal sample are required to properly estimate via the MLE. However, based on the work developed in [START_REF] Pizzolato | Magnitude and complex based diffusion signal reconstruction[END_REF], it is considered that in the case of one noisy realization per sample, the signal samples M i = S(b i ) = S(q i ) can be interpreted as noisy realizations of a function, thus the likelihood is given by

L = Ns i=1 M i σ 2 e -M 2 i +[Ac] 2 i 2σ 2 I 0 [Ac] i M i σ 2 (4.25)
where n is the number of samples of the magnitude signal, M i is the i-th sample, σ 2 is the noise variance, A is the real-valued SHORE design matrix, c is the vector of the even coefficients, and I 0 is the modified zeroth order Bessel function of the first kind. The MLE is then defined as the estimator maximizing L or equivalently log

L ĉML = argmax c (log L) (4.26)
where ĉML is the vector of coefficients estimated with MLE. The maximum likelihood reconstruction is

M M L = N -1 n∈even c M L n φ n (u, q) (4.27)
where {c M L n } are the coefficients. The performance of the estimator is expected to increase with the increasing number of samples of the signal. It should then be noticed that the MLE requires σ 2 to be given. Thus a prior estimation of the variance of the signal is necessary.

MAP-MRI

MAP-MRI is a functional basis developed by Özarslan et al. [2009a, 2013a,b] that extends the concept behind 1D-SHORE to the 3D space, i.e. the possibility of representing a non-Gaussian 3D signal profile. In this sense it generalizes the DTI tensor representation of the diffusion signal and EAP. It has the convenient property that its basis functions are eigenvectors of the Fourier transform [START_REF] Gilbert | Properties of hermite series estimation of probability density[END_REF]. These bases are given as products of three orthogonal one-dimensional basis functions, which are the previously shown SHORE bases as originally formulated by Özarslan et al. [2008a]. The signal attenuation bases are

Φ n 1 n 2 n 3 (A, q) = φ n 1 (u x , q x )φ n 2 (u y , q y )φ n 3 (u z , q z ) with φ n (u, q) = i -n √ 2 n n! e -2π 2 q 2 u 2 H n (2πuq).
(4.28)

and those representing the EAP are

Ψ n 1 n 2 n 3 (A, r) = ψ n 1 (u x , r x )ψ n 2 (u y , r y )ψ n 3 (u z , r z ) with ψ n (u, r) = 1 √ 2 n+1 πn!u e -r 2 /(2u 2 ) H n (r/u) (4.29)
where ψ = FT -1 (φ), and A is the scaling factor matrix to be estimated. In practice, the signal attenuation and the EAP can be described with the same coefficients as

E(q) = Nmax N =0 {n 1 ,n 2 ,n 3 } c n 1 n 2 n 3 Φ n 1 n 2 n 3 (A, q) P (r) = Nmax N =0 {n 1 ,n 2 ,n 3 } c n 1 n 2 n 3 Ψ n 1 n 2 n 3 (A, r) (4.30)
where N max is the maximum basis order and where the second sum in each equation is taken over all the possible combinations {n 1 , n 2 , n 3 } such that n 1 + n 2 + n 3 = N . As indicated by Özarslan et al. [2013b], to fit MAP-MRI bases for the signal attenuation it is convenient to move to an anatomical frame of reference, xyz a , that is the space where the eigenvalues of the DTI tensor are aligned with the coordinate axes. The eigenvalues of the DTI tensor are then used to set the scaling factor in each direction (i.e. {u x , u y , u z }), thus performing an anisotropic scaling of the bases.

The anatomical space DTI tensor, D a , is obtained from the "image space" one, D i ,

D a = VD i V T (4.31)
where V is an orthonormal rotation matrix. The scaling matrix is then estimated as the covariance matrix of displacements [START_REF] Peter | Relationships between diffusion tensor and q-space mri †[END_REF] as

A = 2D a τ = Diag(u 2 x , u 2 y , u 2 z ) (4.32)
where, in case of u x = u y = u z MAP-MRI corresponds to a 3D-SHORE representation as given by Özarslan et al. [2009a]. The bases are fitted up to a maximum radial order N max . In the case of magnitude signal attenuation, a symmetric propagator is estimated by using, as for the case of SHORE in eq. ( 4.24), only the even basis orders and the total number of basis coefficients is N coef = 1 6 (F + 1)(F + 2)(4F + 3) with F = N max /2. A DTI representation is obtained with N max = 0 whereas higher N max allows representing non-Gaussian signal attenuation profiles and EAPs.

Although an estimation of the bases' coefficients can be done via linear least squares, in this thesis it is used a regularized estimation based on an analytic Laplacian as developed by Fick et al. [2016] combined with positivity constraints on the recovered EAP as shown by Özarslan et al. [2013b]. Indeed, this combination allows complying with the nature of the EAP, i.e. a real positive-valued probability density function, and obtaining at the same time a smooth signal profile.

The next section presents the scalar indexes that can be recovered from DTI and MAP-MRI, which are used later on in this thesis.

Scalar metrics

This section describes the scalar metrics computed from the previously presented DTI and MAP-MRI signal/EAP representations. These will be used later on in this thesis for experimental purposes.

DTI scalar metrics are based on the eigenvalues and eigenvectors obtained from the diffusion tensor D as reported in eq. (4.15). The principal eigenvector, v 1 , is the one associated to the maximum eigenvalue, λ 1 , which corresponds to the principal apparent diffusivity, henceforth denoted PD in this thesis, such as PD = λ 1 .

(4.33)

Similarly, the tensor encodes the value of the mean diffusivity, MD, that corresponds to the average of the eigenvalues as [Kärger, 1988, Basser and[START_REF] Basser | Microstructural features measured using diffusion tensor imaging[END_REF]]

MD = λ 1 + λ 2 + λ 3 3 . (4.34)
From the knowledge of the eigenvalues it is then possible to calculate, within the limitations of a Gaussian approximation of the diffusion process, the "magnitude" of the tensor that can be ascribed to anisotropic diffusion, i.e. how much the tensor differs from a spherical shape, via the fractional anisotropy, FA, computed as [START_REF] Basser | Microstructural features measured using diffusion tensor imaging[END_REF]

] FA = 3 2 (λ 1 -MD) 2 + (λ 2 -MD) 2 + (λ 3 -MD) 2 λ 2 1 + λ 2 2 + λ 2 3 . (4.35)
The scalar metrics based on MAP-MRI representation, measure properties of the recovered signal attenuation, or equivalently of the EAP, and particularly of its so-called boundary cases as given by Özarslan et al. [2011, 2013b]. Indeed, these boundary limits can be related to properties of pores in the sense of the pore theory presented by [START_REF] Paul T Callaghan | Pulsed-gradient spin-echo nmr for planar, cylindrical, and spherical pores under conditions of wall relaxation[END_REF], such as the mean pore volume V , cross-sectional area A , and length L . The signal/EAP boundary limits are the return-to-origin (RTOP), return-to-axis (RTAP), and return-to-plane (RTPP) probabilities

RTOP = P (0) = V -1 (4.36) RTAP = R P (r ) dr = A -1 (4.37) RTPP = R 2 P (r ⊥ ) dr ⊥ = L -1 (4.38)
where RTOP is the only one not assuming the white matter tissue to be modeled by parallel cylinders, with the displacements r and r ⊥ being parallel and perpendicular to the cylinder' axis respectively. Other computed features of the signal attenuation/EAP are the mean squared displacement (MSD) [START_REF] Cheng | Estimation and Processing of Ensemble Average Propagator and Its Features in Diffusion MRI[END_REF],

MSD = R 3 P (r) r 2 d 3 r (4.39)
and the q-space inverse variance (QIV) calculated as [START_REF] Pasha Hosseinbor | Bessel fourier orientation reconstruction (bfor): An analytical diffusion propagator reconstruction for hybrid diffusion imaging and computation of q-space indices[END_REF] QIV

-1 = R 3 E(q) q 2 d 3 q. (4.40)
For all of these metrics, closed formulas are derived by Özarslan et al. [2013b] and Fick et al. [2016].

Conclusion

In this chapter, the Fourier relationship between the signal attenuation and the ensemble average propagator (EAP) has been clarified. The exploitation of this relationship will be at the heart of the next three chapters, where the focus will be on the theoretical aspects that relate the real and imaginary parts of the complex signal to the underlying diffusion process. This proposed theoretical investigation will lead to unconventional shapes of the recovered EAP.

Later on, the theoretical investigation will leave place to the exploitment of the acquired complex signal, i.e. complex DWIs, with the goal of performing signal denoising. More precisely, the goal will be the restoration of Gaussianity assumptions about the noise affecting the DWIs, such that any subsequent processing -DTI or MAP-MRI -can be done considering Gaussian noise.

The next chapters, for the diffusion part of this thesis, are then organized as follows.

Chapter 5 exploits further the Fourier relationship between the complex diffusion signal attenuation and the EAP, describing the principles of an asymmetric EAP.

Chapter 6 proposes a new form of exploitment of the Fourier relationship based on the principles described in the previous chapter.

Chapter 7 proposes a theoretical study about the feasibility of acquiring a complex diffusion signal, by studying the effects on a proposed numerical phantom.

Chapter 8 makes a comparison between the benefits of performing diffusion signal reconstruction in the complex domain rather than on the magnitude, investigating on eventual benefits of the former compared to the latter. Then, the chapter proposes a characterization of existing techniques to exploit the complex diffusion signal, showing the benefits and the downsides of it, which are overcome via a proposed strategy.

A summary of the contributions is given in chapter 14.

Chapter 5

The interpretation of complex diffusion

Overview

In diffusion MRI the magnet acquires the signal in each voxel in quadrature, that is a complex signal. However, when looking at real and imaginary DWIs they are of little intelligibility because their phase is dominated by artifacts. But in the hypothesis of having a signal free from these artifacts, some particular tissue microstructural configurations, under certain assumption, can lead to an asymmetric EAP and might, by the exploitment of the Fourier relationship, generate a complex signal attenuation. The chapter explores this concept by reporting some examples found in the literature, and provides interpretations and clarifications. The concepts are then leveraged in the next chapter to present another type of microstructure configuration that exhibits the same complex characteristics.

Keywords complex signal; EAP asymmetry; geometry; Fourier.

Motivations and limitations

When acquiring magnitude signal, the Fourier relationship according to eq. ( 4.10) implies a symmetric EAP. However, if one imagines water spins moving in a microstructural environment, then it is plausible to expect that there will be some cases in which the distribution of displacements will not be symmetric along all of the spatial directions. This generally holds -as it will be clearer in the course of this chapter -as long as the diffusion process sensed within a voxel accounts for diffusing spins that have experienced interactions, during the diffusion time τ , with the microstructural environment also outside the voxel itself. However, in the next chapter it will be proposed a derivation of an asymmetric EAP for a specific microstructure geometry, that would prompt the intuition to find an exception to the previous rule. Adding to this, within the diffusion MRI community it is not completely clear whether a complex diffusion signal could even be measured with a standard PGSE sequence. Section 5.2 will try to bring clarity on this.

Nevertheless, in the present thesis a property of the Fourier relationship is exploited in its full extent. This is the property specifying that the Fourier transform of an asymmetric function, i.e. the EAP, leads to a complex function, i.e. the signal attenuation. The chapter provides examples of complex signal attenuation based on asymmetric EAPs, featuring in details this relationship and providing intuitive interpretations that support the intelligibility of the more rigorous mathematical derivations.

On the complex nature of the diffusion signal

In the vast majority of the literature, the diffusion signal attenuation E(q) is always treated in terms of magnitude attenuation. Nevertheless, the signal acquired in MRI has a complex nature and the relationship between diffusion of spin bearing particles and imaginary contributions to the signal attenuation is not completely understood. Strictly speaking, in the largest amount of cases in clinical diffusion MRI acquisitions, the signal attenuation due to diffusion can be considered to be real-valued. For instance, chapter 8 discusses a denoising procedure, namely phase correction, that is based on this assumption. However, as it will be clearer in the course of this section, the diffusion signal has, in general, both real and imaginary components. The practical interpretation of this complex nature is the object of the rest of this chapter, and at the heart of the contributions in the next two ones. However, now an explanation of the origin of the real and imaginary components will be given following the line of Grebenkov [2007], which is based on pore theory. To do so, it is convenient to repropose the complex signal contribution of the dephasing of an individual spin, previously reported in eq. ( 4.4), e jϕ where j is the imaginary unit, and where ϕ is the total phase accumulated during a time T . The macroscopic signal attenuation E at time T can then be obtained by averaging within the domain as

E = E{e jϕ } (5.1)
where E indicates the expectation. For a general temporal profile f (t) of an acquisition sequence, the total phase ϕ accumulated during a time T is

ϕ = γβT 1 0 f (t)B(X t ) dt (5.2)
where t ∈ [0, 1] denotes a dimensionless time w.r.t. T , β is the maximum intensity of the magnetic field, and B( • ) is the dimensionless spatial profile of the magnetic field, which is normalized to 1 and accounts for the variation of the experienced magnetic field at the time-varying location X t . Note that in eq. ( 5.2) the effect of the 180 • RF pulse is implicitly accounted by opportunely changing the sign of f (t).

One can see that ϕ in eq. ( 5.2) is the product of two terms

p = γβT and φ = 1 0 f (t)B(X t ) dt (5.3)
where the first, p, finds similarities with the q-value presented in eq. ( 4.7), and the second, φ, is a random variable due to the randomness of X t . Equation (5.1) can then be rewritten as

E = E{e jpφ } (5.4)
and since φ is bounded, its characteristic function is analytic and in particular the attenuation can be written as a series expansion of its moments E{φ n } as

E = ∞ n=0 (jp) n n! E{φ n } (5.5)
where it is clear that, in case the odd moments are non-zero and do not cancel out, the signal attenuation has an imaginary component.

In the literature the presence of the odd moments, i.e. the imaginary part of the signal attenuation, is often skipped or not explicitly accounted for. Indeed, demonstrating the cancellation of these moments, as an universal rule, is not a trivial task. In fact, as it will be clearer later, there are counterexamples that show the existence of these moments despite differently conjectured by Grebenkov [2007]. A complete investigation on this subject is beyond the scope of the thesis. However, some considerations about it are here provided revealing that, while there are specific cases in which it is possible to demonstrate the cancellation of the odd moments, in general this is typically not the case. Again Grebenkov [2007] discusses this issue, providing some general guidelines and the demonstration that in one specific case the odd moments cancel out. However, in all the other cases, odd moments might not cancel out, although a clearer light has still to be shed on. In his work, Grebenkov gives some guidelines regarding the first moment, E{φ}. The first moment is argued to be zero, for any sequence temporal profile, f (t), if the magnetic field is such that it integrates to zero over the considered domain, that is if E{B(X t )} = 0. This condition is not met when gradients are non-linear or when the domain is asymmetric. Alternatively, the first moment is zero when the temporal profile, f (t), is such that it satisfies the refocusing condition

1 0 f (t) dt = 0.
(5.6)

However, these conditions are derived w.r.t. the first moment, thus the other odd moments can still be present. A general rule is instead formulated for antisymmetric temporal profiles. In practice, Grebenkov demonstrates that, in a closed pore, all the odd moments cancel out if the temporal profile, f (t), is antisymmetric w.r.t. the central point, that is when the following condition is met

f (1 -t) = -f (t) ∀t ∈ [0, 1/2].
(5.7)

From the condition expressed in eq. (5.7) it is clear that sequences such as the antisymmetric PGSE, shown in fig. 3.6, do not lead to a complex signal attenuation even in the case the microstructural domain shows asymmetry. To clarify what actually is the signal measured in a antisymmetric PGSE, within the context of pore theory, it is possible to look at the signal attenuation in terms of eq. (4.8) in the previous chapter, where it appears clear the relationship with the diffusion propagator P (r 1 |r 0 , ∆), where r has been rewritten as r 1 to remind that here reference is made to a position and not to a displacement. In the case the diffusion time, i.e. ∆ under NPA, is very long (∆ → ∞) and if the pore is closed, then spins can sense the limits of the pore itself and diffusion becomes independent from the initial position. In practice, this reduces to the equality

P (r 1 |r 0 , ∆ → ∞) = ρ(r 1 ) (5.8)
which corresponds to the final spin density function. Thanks to the equality in eq. ( 5.8), eq. ( 4.8) can be written as the product of separete terms

E aPGSE (q, ∆ → ∞) = ρ(r 0 )e -j2πq • r 0 dr 0 ρ(r 1 )e j2πq • r 1 dr 1 (5.9)
where "a" in front of PGSE symbolizes that it is antisymmetric. The product in eq. (5.9) corresponds to the multiplication of the Fourier transforms, ρ(q), of the spin density functions that, since they are equal to each other due to the long diffusion time regime, is

E aPGSE (q, ∆ → ∞) = ρ(q)ρ * (q) = |ρ(q)| 2 (5.10)
where the presence of the second power renders the function symmetric, leading to a real-valued signal attenuation. However, differently form what conjectured by [START_REF] Denis S Grebenkov | Nmr survey of reflected brownian motion[END_REF], if the antisymmetry of the temporal profile expressed in eq. (5.7) is broken then a complex signal is obtained. Particularly, [START_REF] Bernd Laun | Determination of the defining boundary in nuclear magnetic resonance diffusion experiments[END_REF] showed that by using a long and small first gradient pulse, and a short, tall second one after the 180 • RF pulse -while respecting the refocusing condition of eq. (5.6)the signal attenuation can be written as

E naPGSE (q, ∆ → ∞) = e j2πq • rcm ρ(q) (5.11)
where "na" stands for non-antisymmetric, and where r cm is the pore's center of mass. The advantage of the formulation in eq. (5.11) is that the spin density function can be recovered, giving a complex signal in the case the pore is point asymmetric, that is for instance the case of a triangular pore. For reference, other equivalent methods are proposed by [START_REF] Shemesh | Magnetic resonance imaging by synergistic diffusion-diffraction patterns[END_REF] and further extended by Kuder and [START_REF] Kuder | Nmr-based diffusion pore imaging by double wave vector measurements[END_REF].

Although pore theory predicts cancellation of the odd moments or imaginary signal components when the aPGSE sequence is adopted, these have been observed by [START_REF] Liu | Characterizing non-gaussian diffusion by using generalized diffusion tensors[END_REF] and Özarslan et al. [2008b, 2009b] in simulation scenarios. Particularly, Liu and colleagues performed a Monte Carlo simulation showing that spins diffusing within a Y-shaped domain, i.e. an asymmetric domain, generate odd moments or, in other terms, a complex signal with non-zero phase. The particularity of the case lies in the fact that the simulation was carried out by considering a standard antisymmetric PGSE sequence, thus apparently in open contrast with what demonstrated by Grebenkov. The missing keystone in this controversial scenario might consist in the fact that Liu and colleagues performed a simulation where spins, each one with its own phase contribution ϕ, are allowed to enter and exit the voxel during the acquisition time. For this reason, all the spins whose trajectory ends within the considered voxel at echo-time, TE, contribute to the signal. Henceforth, this situation is denoted as leakage. Imaginary signal components have been shown also by Özarslan et al. [2008b, 2009b] who presented simplified aPGSE-based cases where spins contributing to signal have experienced interaction with an asymmetric microstructural environment outside the voxel itself, deriving analytic formulas for signal attenuation and EAP in such cases. Indeed, previous considerations based on the work of [START_REF] Denis S Grebenkov | Nmr survey of reflected brownian motion[END_REF] relate to pore theory, where the domain is closed and does not allow leakage. In the rest of the chapter, as well as in the contributions of chapters 6 and 7, the focus will be on examples that follow the line of Liu, Özarslan, and colleagues.

Diffusion in the proximity of an infinite plate

When considering a voxel located in the proximity of a restricting boundary such as an infinite plate, with voxel's dimensions significantly smaller than the separation distance between the plate and any other boundary, the magnetization in the voxel is influenced only by the boundary in the vicinity [ Özarslan et al., 2008b]. A graphical representation of this scenario is shown in fig. 5.1a, where an infinite plate is represented with its normal aligned with the z-axis, and a voxel with height z 2 -z 1 is located at a distance z 1 from the plate itself. Using the notation in Özarslan et al. [2008b] it is convenient to define the dimensionless position variable

ζ = z u (5.12)
where u is the characteristic diffusion length given by u = 4D 0 ∆ (5.13) with D 0 being the free diffusion coefficient and ∆ the diffusion time. Figure 5.1a also shows the gradient vector which, in this case, is assumed to be aligned with the positive z-axis. The gradient vector magnitude is then given by q = γδG/2π where γ is the gyromagnetic ratio, δ is the diffusion pulse duration and G is the diffusion gradient strength. The considered experiment is the pulsed gradient spin echo (PGSE) sequence with pulse duration δ small compared to ∆. Finally, after defining the dimensionless wave-number κ = πqu, the complex signal is given by

Özarslan et al. [2008b]

E

([ζ 1 , ζ 2 ], κ) = e -κ 2 F (ζ 2 ) -F (ζ 1 ) 2(ζ 2 -ζ 1 ) (5.14)
where ζ 1 and ζ 2 are the dimensionless coordinates of the voxel corresponding to z 1 and z 2 via eq. (5.12), and F (ζ) is

F (ζ) =ζ + 1 √ π e -(ζ-iκ) 2 + e i4κζ i4κ [1 -erf(ζ + iκ)] -ζ -i κ + 1 4κ erf(-ζ + iκ).
(5.15)

To better understand the situation one can look at the isolines drawn, for the only sake of exemplification, around the voxel in fig. 5.1a. These indicate that diffusion in the surrounding of the voxel is not symmetric due to the presence of the plate.

In other words, spins in motion are probing a point asymmetric environment and the diffusion sensed within the voxel, that could be thought as a measuring region, is influenced by this "boundary" effect. This concept was used by Özarslan et al. [2008b] to explain, for instance, the situation in which the voxel of interest, that contains nothing but cerebrospinal fluid, is located in the proximity of the cortical surface of the brain. Here, eq. (5.14) has been used to calculate the predicted complex signal attenuation along the direction perpendicular to the infinite plate.

As it can be seen in fig. 5.1b, the signal attenuation is complex-valued. From Figure 5.2: Circumference of infinitesimal thickness imaged with a 2 × 2 grid of voxel. The voxel in the first quadrant is highlighted for explanation convenience.

Spins are assumed to diffuse only on the circumference. In the rest of the space spins are thought to be either static or not present. The image is adapted from the work of Özarslan et al. [2009b].

studies of eq. ( 5.14), it can be seen that the imaginary component tends to vanish as the considered voxel size, z 2 -z 1 , or the distance from the plate, z 1 , increase.

Although the resolution in a clinical setup is not compatible with the observation of this phenomenon, in experimental setup was possible to demonstrate effects on the magnitude attenuation predicted by eq. (5.14).

The understanding of the generation of a complex attenuation can be better achieved with the example presented in the next section. This example, firstly presented in the work of Özarslan et al. [2009b], is functional to obtain an empirical intuition of why an imaginary signal component can occur.

Diffusion on a circumference

This section describes an example of diffusion on a circumference that well renders the concept of an asymmetric ensemble average propagator (EAP). After describing the mathematical derivation as summarized from the work of Özarslan et al. [2009b], a newly introduced intuitive explanation will be provided to get familiar with the concept of asymmetry and its influence on the complex signal attenuation via the Fourier relationship.

Let us consider a circumference of infinitesimal thickness where the spins only lies on the circumference itself, the rest being completely empty or with no motion due to diffusion. Then, let us imagine that a MRI image is acquired in a 2 × 2 matrix, each voxel defining a measuring region that covers one quarter of the circumference as depicted in fig. 5.2. The idea is to define the EAP and the signal attenuation of the upper right voxel that is highlighted in the figure, which probes the diffusion coming from a curving fiber geometry. As in the work of [START_REF] Liu | Characterizing non-gaussian diffusion by using generalized diffusion tensors[END_REF], spins are thought of being able to move on the whole circumference during the acquisition time, thus allowing leakage. The formal mathematical derivation of the EAP in such a scenario is explained by Özarslan et al. [2009b]. For this thesis the whole formulation has been re-implemented with the purpose of exploration and understanding. Indeed, the understandingThanks to Nathanaël Foy for the brainstorming sessions. of this example made possible the extension to the concepts therein in the proposed examples, that apply to different geometries and scenarios, introduced in chapters 6 and 7. For this reason, the adopted procedure is now summarized. The apparent propagator for a displacement r at a voxel location X is given by [ Özarslan et al., 2009b] P (X, r) = ρ(X -r)

R 3 ρ(r )K(r ; X, ∆) dr K(X -r; X, ∆) (5.16)
where the initial spin density for the circular geometry is ρ(r 0 ) = δ(r 0 -a) 2πr 0 (5.17) with r 0 and φ 0 being the polar coordinates of r 0 , and a the radius of the circumference. The propagator K( • ) is

K(r 0 ; r 1 , V ) = δ(r 1 -r 0 ) πr 1 1 2 + ∞ m=1 e -m 2 V 2 cos m(φ 1 -φ 0 ) (5.18)
where again r 1 and φ 1 are the polar coordinates of r 1 , and V is a dimensionless variable that can be thought as an equivalent of the diffusion time. Indeed, for V → 0 the propagator tends to a delta function, whereas for V → ∞ it tends to a uniform distribution over the circumference. In this thesis, the EAP expressed by eq. (5.16) has been evaluated numerically, truncating the infinite sum in eq. ( 5.18) to an integer m max . The corresponding signal attenuation has been calculated by exploiting the inverse of eq. (4.11). Images in fig. 5.3 show the EAP obtained with eq. ( 5.16) for different diffusion times (V ∈ [0.5, 5000]). Note that the obtained EAP is point asymmetric, and intuitively the amount of asymmetry increases with the diffusion time. The EAP for a circumference is calculated via eq. ( 5.16). By taking its Fourier transform a complex signal is obtained, and from this a magnitude signal attenuation can be calculated. The inverse Fourier transform of the magnitude signal leads do a different EAP, which is symmetric. Colormaps of EAPs and complex signal have been clipped, with respect to the actual extent, for visualization purposes.

remains closed to the center of the displacement space. This phenomenon will be better observed in the theoretical experiment proposed in chapter 7.

Asymmetry and complex signal

The example discussed in the previous section, that reproduces the work of Özarslan et al. [2009b], should be considered as a toy example that helps visualizing how an asymmetric geometry can affect the EAP and the signal attenuation. Although it remains theoretical, due to the peculiar assumptions therein, it sheds light on the intuition behind the generation of the EAP for a specific a geometry.

Here, it is proposed to stand on the teaching given the previous example to realize that when the EAP is asymmetric the corresponding signal is complexvalued. Indeed, under NPA, the signal attenuation can be calculated as the Fourier transform of the EAP, according to eq. (4.11) (regardless the sign convention). The first row of fig. 5.5 shows the result of this procedure, obtained numerically, where it is clear that the outcome of the transformation has a real and an imaginary part. It is interesting to see what would happen in the case only the magnitude of the signal was considered, as normally done in diffusion MRI. Indeed, by taking the inverse transform of the magnitude, which is positive real and symmetric, one finds a real valued EAP that is different from the original asymmetric one. This is shown in the second row of fig. 5.5. Note how both of the calculated EAPs, corresponding to a infinite diffusion time regime, report the circular geometry. However, it appears clear the the EAP computed from the magnitude signal is symmetric.

The observation reported in this section is at the heart of what proposed in the next chapters. However, it has to be understood that this observation is uniquely derived from the properties of the Fourier transform, which in turn is applicable under NPA. Therefore, whether the complex signal attenuation could actually be measured in a MRI sequence it is subject to a different type of considerations, partly reported in section 5.2, which is beyond the scope of this thesis. Therefore, it follows a list of the assumptions considered by this circumference example, and by those proposed after in this thesis. These are:

diffusing spins can enter and exit the measured voxel contributing to the final signal, i.e. leakage occurs;

the spin density can be considered uniform within the geometry interested by contiguous boundaries, and constant over time;

a complex-valued signal attenuation can be measured, phase artifacts apart, with a PGSE sequence and under the previous listed conditions;

the Fourier relationship holds, i.e. validity of the narrow pulse approximation.

If at least one of the previous assumptions does not hold, then the validity of the examples proposed in the following has to be questioned. Moreover, even under the previous assumptions, EAP asymmetry could be observed in vivo admitted that asymmetric microstructural environments exist and that their contribution survives the averaging phenomenon that necessarily happens in MRI, due to the large volume of the tissue probed by a voxel. Nevertheless, according to the previous assumptions, in the next chapter it is proposed a plausible example of microstructure, consisting of compressed sinusoidal axons, which shape induces intrinsic asymmetry in the EAP.

Conclusion

This chapter presented an overview of the interpretation of the complex diffusion MRI signal, and tried to bring clarity with regard to the assumptions, hypothesis, and controversies in the field. Nevertheless, the next two chapters present the contributions of this thesis on further exploring the relationship between asymmetries of the diffusion process and intrinsically asymmetric geometries, or peculiar experimental conditions. Particularly, chapter 6 introduces the derivation of the EAP for tortuous axons where the intrinsic asymmetry of the geometry leads to asymmetric EAPs. Then, chapter 7 proposes a different paradigm to induce diffusion asymmetry. This is exploited for the generation of a theoretical phantom on which analysis of the EAP's asymmetric properties can be performed. Later on, this theoretical investigation will leave place to the exploitation of the complex signal to perform denoising, the other contribution reported in chapter 8, where eventual asymmetric properties will be neglected.

Chapter 6

The axonal compression model

Overview

This chapter further elaborates the concept of EAP asymmetry introduced in the examples of the previous chapter, under the assumptions reported in section 5.5.

Particularly, it proposes a new axonal geometry that could induce asymmetry in the underlying diffusion process, and a method to quantify it. The chapter also proposes a mathematical derivation of the EAP in such a geometry, that is also supported by an empiric explanation for the sake of intelligibility. Finally, simulated results are presented together with considerations about the limitations of the proposed approach.

Keywords tortuosity; compression; intrinsic asymmetry; axon.

Introduction

White matter (WM) tissue may present different microstructural organizations depending on its nature and pathological condition. For instance, axons composing the tissue can be straight as shown in fig. 6.1 left. However, axons may appear sinusoidal due to pathological conditions such as those caused by cervical cord injury. In fact, in these cases not only the shape can be sinusoidal, but longitudinal elongation or compression can result in changes of undulation amplitude as studied by [START_REF] Nilsson | The importance of axonal undulation in diffusion mr measurements: a monte carlo simulation study[END_REF], and shown in central fig. 6.1. Moreover, to a more general extent, variations in the axonal tortuosity might appear as a consequence of the compression. Indeed, compressed axons are partially convoluted and show irregular undulation with a non-uniform tortuosity along the longitudinal direction. Diffusion MRI allows probing such a WM microstructural environment in which spins move.

Under the assumptions seen in the previous chapter, especially in section 5.5, the complex-valued signal attenuation can then be related to the average displacement probability of spins, i.e. the ensemble average propagator (EAP).

In the work of [START_REF] Nilsson | The importance of axonal undulation in diffusion mr measurements: a monte carlo simulation study[END_REF] the undulated axons are modeled as sinusoids, and a method to calculate the corresponding EAP is provided. An interesting consequence of this work was to show that in the presence of compression, such as in the situation depicted in central fig. 6.1, the undulation is a confounding factor for axonal diameter estimation from DWIs, such as in methods introduced by [START_REF] Assaf | Axcaliber: a method for measuring axon diameter distribution from diffusion mri[END_REF] and Alexander [2008]. Indeed, the presence of undulation causes a broadening of the EAP perpendicularly to the longitudinal direction along which the undulated axon develops, thus causing the overestimation of the axonal diameter when making the assumption that no undulation occurs.

This chapter contributes by extending the sinusoidal axonal model of the EAP by including compression, that consists in a variation of the axonal tortuosity along the development of the axon itself as represented in the right image of fig. α > 0, where α indicates the tortuosity variation rate. Secondly, it is noted that such a compressed undulated axon defines a point asymmetric geometry which would result in a asymmetric EAP according to the assumptions presented in the previous chapter. This time, however, the proposed compressed sinusoidal geometry, as opposed to the circular geometry discussed in section 5.4, offers no points of internal symmetry1 , thus it is here denoted as an intrinsically asymmetric geometry. This will lead in a slightly different intuition behind the derivation of the EAP, compared to the one seen for the circular geometry. Under the assumption of having the complex signal attenuation (see section 5.5), it will be shown that the amount of compression, or tortuosity variation α, can be related to the amount of asymmetry detected in the EAP. It is here proposed to quantify asymmetry as the Hellinger distance [START_REF] Hellinger | Neue begründung der theorie quadratischer formen von unendlichvielen veränderlichen[END_REF] between each EAP and its axially reflected version. The chapter also presents experiments elucidating the link between EAP asymmetry and axonal tortuosity rates, thus compression degrees. The content of this chapter is based on the work developed in Pizzolato et al. [2015b].

Tortuous axon model

This section presents the extension of the model for undulated axons, that was introduced by [START_REF] Nilsson | The importance of axonal undulation in diffusion mr measurements: a monte carlo simulation study[END_REF], to include axonal compression, i.e. a variation of the tortuosity along the longitudinal direction of the axon itself. This is illustrated in the right image of fig. 6.1. The proposed model leads to a generalization of the existing model to account for compression, but at the same time it constitutes the definition of an intrinsically asymmetric geometry that could potentially be found in the WM tissue.

The undulated axon is modeled as a sinusoid that lies in the xz-plane and evolves along z with oscillations along x. The proposed model is constituted by an axon of infinitesimal thickness that goes through two full periods of oscillations within the interval [-Z, Z] with Z > 0. Mathematically it is described as

x = f (z) = A sin 2π(z -Z) L(z, α) (6.1)
where A and L(z, α) are respectively the amplitude and the wavelength of the sinusoid, with 0 < A ≤ X and L(z, α) > 0. The tortuosity variation is achieved by considering a varying wavelength according to

L(z, α) = α(z + Z) + l, f or α > 0 (6.2)
where α is the spatial rate regulating the wavelength increment along the z-axis, or alternatively, the tortuosity variation rate. The wavelength increases progressively along the z-axis starting from a minimum value l in z = -Z. The rate α is related to the degree of intrinsic asymmetry induced in the undulated shape along the z-axis.

A value of α = 0 defines a perfect sinusoid, corresponding to the model presented by [START_REF] Nilsson | The importance of axonal undulation in diffusion mr measurements: a monte carlo simulation study[END_REF]. However, in order to distinguish between straight axons and undulated ones, in this chapter experiments will be performed considering the case where α = 0 and A = 0, thus identifying the straight axons as shown in the right image of fig. 6.1. In the next section, the EAP for such as undulated and compressed geometry will be derived. To do so, it is useful to define the curve length for the model defined in eq. (6.1). In particular the length between z and z + ∆z is given by

s(∆z; z) = z+∆z z 1 + A cos 2π(ζ -Z) L(ζ, α) 2π(l + 2αZ) L(ζ, α) 2 2 dζ (6.3)
where ∆z is a displacement along the z-axis with respect to the coordinate z.

Ensemble average propagator for tortuous axons

This section presents the derivation of the EAP for the axonal model proposed in section 6.2. The EAP for a spin displacement (∆x, ∆z) subject to a diffusion time τ is [Tanner and [START_REF] John | Restricted self-diffusion of protons in colloidal systems by the pulsed-gradient, spin-echo method[END_REF]]

P (∆x, ∆z; τ ) = ∞ -∞ ρ 0 (x, z)P (∆x, ∆z|x, z, τ ) dxdz = ∞ -∞ ρ 0 (z)P (∆x, ∆z|z, τ ) dz (6.4)
where ρ 0 is the initial spin density and P is the displacement pdf, i.e. the propagator. The second equality in eq. ( 6.4) is due to the proposed model parametrization in eq. (6.1). As suggested by [START_REF] Nilsson | The importance of axonal undulation in diffusion mr measurements: a monte carlo simulation study[END_REF], ρ 0 is computed by assuming an initial uniform spin density along the axonal model

ρ 0 (z) = cos arctan f (z) -1 (6.5)
where f (z) is the derivative of f in eq. ( 6.1). To produce an approximation for the displacement probability P in eq. ( 6.4), the propagator is decomposed as

P (∆x, ∆z|z, τ ) = P (∆x|∆z, z)P (∆z|z, τ ) (6.6)
where it is assumed that spins diffuse freely along the curve, i.e. they follow a normal distribution with variance 2Dτ where D is the self-diffusion coefficient of water. Then, the probability of a spin experiencing a net displacement along the z-axis ∆z can be approximated as

P (∆z|z, τ ) = 1 √ 4πDτ e -s(∆z;z) 2 4Dτ (6.7)
where s is the arc length of the axon between z and z + ∆ z as given in eq. ( 6.3).

Since the coordinate x on the axon is determined by z, as specified by eq. ( 6.1), the displacement density along x, provided to know ∆z and z, is [START_REF] Nilsson | The importance of axonal undulation in diffusion mr measurements: a monte carlo simulation study[END_REF]]

P (∆x|∆z, z) = δ D (f (z) + ∆x -f (z + ∆z)) (6.8)
where δ D is the Dirac delta function. Given eqs. (6.5), (6.6), (6.8) and (7.4), the EAP of a tortuous axon, according to the proposed model in section 6.2, can be obtained by computing P (∆x, ∆z; τ ) via numerical integration of eq. (6.4). Section 6.4 gives the complex signal formalization for the derived EAP, expressing the actual Fourier relationship existing between the two. Moreover it proposes the a method to quantify EAP asymmetry. However, before going to discuss that and the performed experiments, it is here proposed an intuitive and empiric explanation of the generation of the EAP for the case of a compressed undulated axon, according to the derivation proposed in the present section.

An empiric description of the EAP

Similarly to what done for diffusion in a circular geometry in section 5.4, here it is provided an empiric explanation of how the EAP corresponding to diffusion within undulated and compressed axons could be thought. However, differently from what observed for the circular geometry, here the explanation does not need considering different voxels, that in the former case where covering a considerable portion of the trajectory that could contribute to the signal. Indeed, the intrinsic asymmetry of the model proposed in section 6.2 allows considering a single voxel containing the geometry. The explanation is presented in fig. 6.2, where the case of a straight axon are "dragged and dropped" onto the axes of the displacement space, with the origin corresponding to zero displacement. The EAP can be thought as the integration average of the local propagators along the whole domain, i.e. trajectory, contained within the voxel. As an exemplification, images (e,f,g) reiterate the procedure for the other remaining three sampling points, adding at each time the new local propagator to the displacement space, to be averaged at the end of the procedure. Note that while for the straight axon the same local propagator is added at each step, in the case of the compressed axon different propagators are added resulting in a point asymmetric EAP. The Fourier transform of such a function, will lead to a signal attenuation with non-zero phase. This relationship, as well as the calculation of the EAP asymmetry are now introduced in the following section.

Diffusion signal and EAP asymmetry

As seen in chapter 4 the signal acquired in DW-MRI can be represented as function of the diffusion time τ and frequency vector q = γδG/2π where γ is the gyromagnetic ratio, δ is the diffusion pulse duration and G is the diffusion gradient vector.

The measured q-dependent signal S(q) has a complex nature and its attenuation E(q) = S(q)/S(0) is related to the EAP via a the Fourier relationship, reported in eq. (4.11), that in the considered bidimensional case is

E(q x , q z ; τ ) = ∞ -∞
P (∆x, ∆z|τ )e -j2π(qx∆x+qz∆z) d∆z d∆x (6.9) which is valid under NPA, and where q x and q z are the components of the frequency vector q. The attenuation in eq. ( 6.9) has in general a complex nature.

To measure the asymmetry of the EAP obtained via inverse Fourier transform, it is proposed to use the distance between the EAP itself and its reflected version with respect to the origin. Despite different measures are possible, after defining the displacement vector r = (∆x, ∆z), it is proposed to adopt the Hellinger distance [START_REF] Hellinger | Neue begründung der theorie quadratischer formen von unendlichvielen veränderlichen[END_REF] because it constitutes a proper metric between probability measures

H = 1 √ 2 P (r|τ ) -P (-r|τ ) 2 (6.10)
where 0 ≤ H ≤ 1, 0 corresponding to equality and 1 to maximum inequality. In other words, a value H = 0 corresponds to the situation of complete EAP symmetry, whereas a positive value corresponds to asymmetry.

Experiments

This section presents experiments showing the relationship between EAP asymmetry and axonal tortuosity for the EAP obtained as described in section 6.3.

α =0 α =1 α =2 α =3 α =4 α =5 α =6 α =7
Us ua l EAP from m a g nitude s ig na l α =8 The EAP is calculated via evaluation of eq. (6.4) for the straight, and for compressed axons with increasing tortuosity rates α ∈ [START_REF] Marco | Elucidating Dispersion Effects in Perfusion MRI by Means of Dispersion-Compliant Bases[END_REF]8]. Axons present undulation amplitude A = 4 µm with a basal wavelength l = 50 µm [START_REF] Nilsson | The importance of axonal undulation in diffusion mr measurements: a monte carlo simulation study[END_REF]. The diffusion coefficient is set to D = 2 × 10 -9 m 2 s -1 [Hall and Alexander, 2009] and τ = 28.6 ms. The maximum observable displacement in both directions is set to X = Z = 50 µm, where the observation frame, i.e. the voxel, corresponds to the xz-plane limited to the interval [-X, X] along x and [-Z, Z] along z. As α increases, then eq. ( 6.2) generates an undulated axonal trajectory that is more and more compressed towards the frame's negative bound in -Z (see fig. 6.1 right). The observation frame is discretized over a 257 × 257 displacement grid, and the EAP is calculated for each sample via numerical evaluation of eq. (6.4). Clearly this corresponds to an extremely high resolution of the EAP, which in turn corresponds to an extremely high maximum q-value.

The complex signal attenuation is obtained, for each considered α, as the Fourier transform of the EAP over the displacement grid according to eq. (6.9). The EAPs are then calculated as the inverse Fourier transform of both the complex and magnitude signal attenuations.

Images in fig. 6.3 show the log-transformed EAPs calculated for the straight axon (α = 0) and for different degrees of axonal compression (α > 0). The first row contains the EAPs as usually calculated from the magnitude signal. The second row shows the EAPs obtained from the complex signal. Note that no difference exists between the first and second row in the case of the straight axon. Despite the shape of the magnitude-based EAP changes as α increases, it always shows axial symmetry. On the other hand, evident asymmetry is present in the case of the complex-based EAPs for α > 0. When looking at the EAPs in the second row of fig. 6.3, the difference between the upper and lower lobes of each EAP clearly appears. It also appear clear that the degree of asymmetry between the lobes is more marked at high α, for instance the difference between the lower and upper lobes for α = 8 is more accentuated than in the case of α = 1.

The previous observation motivates to quantify asymmetry of the EAPs. This is done by calculating the proposed Hellinger distance, eq. (6.10), as function of α, as shown in fig. 6.4. In fact, it is found an incremental trend between the complexderived EAP asymmetry and the tortuosity rate α, in agreement with images in the second row of fig. 6.3.

To summarize the results, it is noted that EAP asymmetry increases with the amount of intrinsic asymmetry of the geometry indicated by the tortuosity rate, i.e. the amount of compression of the undulated axon. As specified earlier, the simulations have been carried out with a high resolution (257 × 257 grid for 100 µm displacement space) which is not feasible in real acquisitions. An interesting study would be to quantify the expected decrease of the measured asymmetry as the resolution in r-space, i.e. the q-value, decreases. In fact, lower resolutions of the images in fig. 6.3 would signify that a r-space pixel corresponds to the average of a larger portion of EAP. In other words, the EAPs would look smoother and more blurry, intuitively leading to a decrease of the differences between the upper and lower lobes (see the second row of fig. 6.3). At the same time, one could conjecture that by letting spins to diffuse a longer time τ , the amount of measured asymmetry would increase since spins have more time to interact with the surrounding geometrical environment, thus probing its eventual asymmetry. This phenomenon will be studied in the next chapter. However, assuming the validity of the presented simulation under the assumptions presented in sections 5.2 and 5.5, there are some objective factors that make the presented simulation a convenient tool for theoretical investigation rather than a theory that could be applicable in real practice. For instance, in real scenarios axons have no infinitesimal thickness. Maximum q-values achievable in an NMR experiment might not be sufficient to observe the phenomenon. In addition, it is assumed to be able to acquire a complex-valued signal, where the phase contribution is purely due to diffusion, whereas in real practice the phase is dominated by motion artifacts, as it will be shown in chapter 8. Nevertheless, chapter 7 will continue on the analysis of the asymmetric EAP by exploiting an experimental rather than geometrical paradigm, as proposed source of diffusion asymmetry.

Conclusion

This chapter extends further the concept of diffusion asymmetry by introducing a potentially existing axonal geometry that has the property of being intrinsically asymmetric, and by proposing a way of quantifying such asymmetry. Although the proposed simulation has the characteristics of a theoretical investigation rather than a prediction of the real practice, it might prompt a different understanding of diffusion processes compared to the classic point of view. Guided by the desire of measuring asymmetry in a real NMR experiment, the next chapter presents a phantom designed to generate asymmetry that could potentially be measured. However, difficulties in realizing such a phantom and the strong assumptions behind its design reduce its experimental feasibility. Nevertheless, it provides another tool for theoretical investigation, that can help understanding better the phenomenon of asymmetry.

Introduction

As seen in the previous chapters, the average spin displacement probability, i.e. the EAP, may show -under certain assumptions -asymmetry due to the characteristics of tissue geometry. However, the task to test the existence assumptions of asymmetry, and eventually measure it in experimental conditions, is non-trivial. To test whether -in the case the assumptions were proven to be valid -EAP asymmetry could be measured, it is here presented an in silico study. This consists on the concept of a phantom to measure EAP asymmetry via DW-MRI, designed to highlight the relationship between the controllable experimental parameters and said asymmetry.

Technical difficulties in building a physical phantom for DW-MRI mainly consist on designing axonal geometries at micrometer scale with controllable properties and tissue-like diffusion characteristics. For instance, at the present moment building a phantom resembling compressed axons with different degrees of compression is very challenging. An alternative solution to induce diffusion asymmetry needs to be found. An example could be to reproduce the conditions for diffusion in the proximity of an infinite plate (see section 5.3) that were feasible as described by Özarslan et al. [2008b]. Adding to this, EAP asymmetry may only be retrieved from the complex DW signal thus implying that a problem of observabilty of the phenomenon also arises, specially with reference to the high sensibility of the signal's phase to noise and bulk movement. This chapter, however, mainly discusses the physical principles underlying the proposed phantom and investigate the relationships with the controllable parameters of a Pulsed Gradient Spin Echo (PGSE) sequence, as developed in Pizzolato et al. [2016f].

In what follows, it is presented an experimental setup designed with the precise scope of inducing diffusion asymmetry in the sampled tissue, such that it is possible to predict the corresponding complex diffusion signal. The adopted strategy consists in applying a temperature gradient in a spinal cord tract, along the direction of the fibers, to obtain a spatially localized diffusion coefficient. Henceforth, this is denoted as a temperature-induced diffusivity gradient. The diffusion coefficient will be forced to vary monotonically with the position coordinate, i.e. it will increase as the spatial location gets closer to the highest temperature position. In this way, at each point along the fiber, a water spin moving towards the highest temperature position will experience higher diffusivity values compared to a particle moving in the opposite direction. In general, this would impose skewness of the local propagator, that once integrated over the whole domain would result in an asymmetric EAP.

Although the ultimate purpose of the work reported in this chapter is to assess the feasibility of the temperature phantom, the main contribution consists on the definition of a simplified framework that allows to study the possible generation of EAP asymmetry only based on experimental conditions, i.e. independently from a specific geometry such as that proposed in chapter 6. Once again, the proposed example relies on the assumptions reported in section 5.5.

The chapter first describes the proposed phantom and the related assumptions. A temperature gradient applied to the phantom will be considered, and the corresponding spatial diffusivity profile induced in the phantom will be calculated. The EAP and the signal attenuation are then derived for simulated voxels located along the longitudinal direction, that is the direction aligned with the fibers. The quantification of the induced asymmetry in the EAP will be discussed as the experimental conditions and acquisition parameters variate, particularly highlighting the relationship with the temperature-induced diffusivity gradient. Finally, the feasibility of the proposed technique and its inherent limitations will be discussed.

Phantom design and experimental setup

The simulation considers a temperature gradient applied to a longitudinal spinal cord tract of length l, where the longitudinal direction corresponds to that parallel to the fibers. The temperature gradient is generated in the spinal cord by heating up the two extremities at two different temperatures maintained constant during time: one extremity at low temperature T L and the other at high temperature T H . The schematic representation of the experimental setup is illustrated in fig. 7.1. The corresponding temperature profile along the spinal cord tract is consequently obtained, as discussed later in section 7.3.

The temperature difference between the extremities produces a profile T (z) of varying temperature along the longitudinal direction z. The spatially-localized temperature T (z) in the spinal cord tract is used to obtain the corresponding diffusivity value D(T ). The diffusivity values are computed along the spinal cord in The DW-MRI acquisition is simulated bidimensionally, similarly to the derivation shown in chapter 6. The spinal cord tract is discretized in adjacent spatial units corresponding to voxels distributed along the longitudinal direction (fig. 7.1). The tissue underlying each voxel is considered as being composed of straight fibers with infinitesimal thickness. Water particles diffusing within these fibers are subject to the locally observed temperature-induced diffusivity gradient. Therefore the simulated signal takes into account the local values of the diffusion coefficient.

The simulation takes into account voxels size d vox (m) and PGSE parameters such as maximum diffusion gradient strength G max (T /m), pulse duration δ (s) and separation between pulses ∆ (s). Finally, the EAP asymmetry is calculated for each voxel and studied as function of the applied temperature gradient. The following section provides details about the assumed temperature gradient along the longitudinal direction of the spinal cord tract.

Applied temperature gradient

This section presents the assumptions made for the temperature profile applied along the spinal cord tract. The latter is considered as a homogeneous rod extending longitudinally to the fibers' direction. Due to the long time required for the MR acquisition, stationary conditions are assumed. Therefore no transient dynamic is considered for the calculation of the resulting temperature gradient which will then be constant. Hence, the temperature gradient along the rod is obtained as the steady-state solution of the heat equation when the two extremities are maintained at two different constant temperatures, rendering a linear temperature profile

T (z) = z • T H -T L l + T L (7.1)
where T L is the temperature in z = 0, T H is the temperature applied in z = l, and l is the length of the spinal cord tract as shown fig. 7.1. Henceforth, temperature will be expressed in Celsius degrees ( • C).

The following section reports the calculation of the diffusivity profile along the spinal cord tract, for different tissue types, given the temperature profile calculated with the presented eq. (7.1).

Induced diffusivity profile

In order to calculate the propagator, it is important to determine the local value of the diffusion coefficient, which varies spatially according to the applied temperature profile. However, the diffusion coefficient depends also on the characteristics of the liquid in which particles diffuse and on the geometrical properties of the particles themselves. These dependencies are well described by the Stokes-Einstein equation:

D(T ) = k B (T + 273.15) 6πRη (7.2)
where k B is the Boltzmann constant, R the radius of the spherical particle and η the dynamic viscosity. Assuming free diffusion of water particles, the radius can be approximated by applying the inverse formula of eq. ( 7. is also a function of the temperature and cannot simply be a constant. Therefore, to calculate η 25 f w it is here employed the equation that describes the relationship between free water viscosity and temperature given by Al-Shemmeri [2012],

η f w (T ) = 2.414 • 10 -5 • 10 247.8 T +273.15-140.0 (7.3)
which is accurate to within 2.5% from 0 • C to 370 • C. Using eq. ( 7.3) the predicted free water viscosity at 25 • C is η 25 f w = 0.891 • 10 -3 Kg/ms. The particle's radius is then approximated by inverting eq. (7.2), obtaining R ≈ 1.066 -10 m.

Diffusivity in tissues is lower than that of free water, considering the same temperature conditions. Indeed, dynamic viscosity increases in tissue according to the microstructural properties therein. Assuming a constant temperature of 37.9 • C in the living brain tissue [START_REF] Kozak | Using diffusion mri for measuring the temperature of cerebrospinal fluid within the lateral ventricles[END_REF], the dynamic viscosity in CSF (η 37.9 csf ), GM (η 37.9 gm ) and WM (η 37.9 wm ) can be calculated from eq. ( 7.2) knowing the corresponding diffusivity values, which can be found to be D 37.9 csf = 2.9 • 10 -9 m 2 /s, D 37.9 gm = 0.89 • 10 -9 m 2 /s, and D 37.9 wm = 0.73 • 10 -9 m 2 /s as given by [START_REF] Helenius | Diffusion-weighted mr imaging in normal human brains in various age groups[END_REF].

Having the viscosity values for the different tissues at 37.9 • C, it is possible to calculate an empirical scaling factor between free water viscosity, as expressed by eq. ( 7.3), and tissue viscosity. Therefore, viscosity values could be found to correspond to η csf ≈ 1.09 • η f w , η gm ≈ 3.54 • η f w and η wm ≈ 4.32 • η f w . Note that the scaling factor calculated for gray matter is in agreement with the value assumed in the work of [START_REF] Nicolas | Water diffusion in q-space imaging as a probe of cell local viscosity and anomalous diffusion in grey and white matter[END_REF]. These factors are adopted to scale eq. ( 7.3) to calculate the dynamic viscosity at all the temperatures for the different tissue types. By substituting η in eq. ( 7.2) with the opportunely scaled version of eq. ( 7.3), it is possible to compute the diffusivity profiles along the spinal cord tract, as shown in fig. 7.3. It is here pointed out that the diffusion coefficients used to obtain the scaling factors are the apparent diffusion coefficients for the corresponding tissues. Particularly the apparent diffusion coefficient in a tissue can be seen as the product λ • D * where D * is the diffusion coefficient of the media, that depends on temperature and viscosity, and λ is the tortuosity coefficient which solely depends on geometry. For instance, in WM the tortuosity changes with the considered diffusion direction, i.e. the diffusion coefficient along the direction parallel to fibers is higher than D 37.9

wm . However, the latter value is considered with the awareness that it renders a lower bound for the diffusivity profile along fibers. Moreover, the diffusivity profile of CSF and GM may only be considered within the fibers for the sake of representing 1D diffusion in the corresponding tissues. Now that the spatial diffusivity profile along the phantom is known, it is possible to calculate the EAP for each adjacent voxel, and obtain the corresponding complex signal attenuation.

EAP derivation

This section presents the derivation of the Ensemble Average Propagator (EAP) when a temperature-induced diffusivity gradient is applied to the spinal cord tract. The obtained EAP is used to recover the complex DW signal via Fourier relationship. From the complex signal the original EAP can be recovered exploiting the inverse relationship and its asymmetry can be quantified.

The EAP is obtained in 2D by considering diffusion within straight fibers with infinitesimal thickness. Fibers are considered to be aligned along the longitudinal direction of the spinal cord tract. Assuming that particles diffuse freely along the fiber, their displacement follows a normal distribution with variance 2D(T, η)τ , where τ is the diffusion time and D(T, η) is the diffusion coefficient as function of temperature and dynamic viscosity. Then, the probability of a particle experiencing a net displacement, ∆z, along the z-axis can be approximated by

P (∆z, τ, T, η) = 1 4πD(T, η)τ e -∆z 2 4D(T,η)τ (7.4)
which depends on the temperature and, consequently, on the location of the particle within the spinal cord tract. Note that since both the temperature and the dynamic viscosity are function of the particle location z, eq. ( 7.4) can be rewritten as

P (∆z|z, τ ) = 1 4πD(z + ∆z)τ e - ∆z 2 4D(z+∆z)τ (7.5)
where the dependency of the diffusivity on the location is explicit. The fact that the diffusion coefficient is not constant but varies monotonically with the location is the cause of the asymmetry of eq. (7.5), which will be reflected in the EAP.

The EAP that will be considered for calculating the signal is that accounting for the ensemble of particles within the whole voxel. To do so, the length of the spinal cord tract l is discretized in locations z i each 1µm. At each location z i eq. (7.5) is evaluated over an observation frame large enough to observe the Gaussian decay, thus obtaining a local propagator. Particularly the local propagator is calculated from a maximum negative displacement -∆z max to a maximum positive displacement ∆z max with ∆z max = (6D 37.9 f w τ )

1 2 . The final EAP for each voxel is obtained by numerical integration of the local propagators calculated for the locations z i within the voxel's limits, considering voxels of size d vox (= 1 mm in the rest), corresponding to the numerical implementation of

P (∆z; τ ) = z∈voxel ρ 0 P (∆z|z, τ ) dz (7.6)
where ρ 0 is the constant initial density of particles. Note that it is assumed that the density is constant, which is an approximation. The resolution of the calculated EAP, i.e. the minimum observable displacement r, is calculated from the maximum q-value as r = 1/q max . Sequence parameters such as δ and ∆ are chosen in agreement with plausible real values to account for the chosen τ . The EAP is used to calculate the complex signal via Fourier transform under NPA with eq. (6.9).

The final goal is to quantify the asymmetry of the EAP calculated for each voxel. To do so, the EAP is re-computed from the complex signal and the EAP asymmetry is calculated as the Hellinger distance between each EAP and its axially reflected version according to eq. (6.10), thus obtaining a value H = 0, corresponding to symmetry, or a value 0 < H ≤ 1 corresponding to asymmetry.

Results

Results show that the amount of EAP asymmetry increases as the temperature gradient increases. Particularly, acquisitions with increasing T H (i.e. 45, 55, 65, 75, 85 • C) while keeping T L fixed (i.e. 25 • C) will lead to an increasing EAP asymmetry. This relationship can be predicted with the developed model diffusivity profile is quasi-linear, is constant. The only actual differences between voxels have to be ascribed to the slight non-linearity of the temperature-induced diffusivity gradient (fig. 7.3).

EAP asymmetry increases with the temperature-induced diffusivity gradient. The most interesting result is that, despite EAP asymmetry across voxels remains substantially constant, the amount of asymmetry increases with the temperature gradient (for a given dynamic viscosity or tissue type). For instance, the right side of fig. EAP asymmetry increases with the diffusion time. When particles can diffuse for a longer time they can probe a wider range of diffusivity values. Therefore, the differences between the EAP values calculated for positive and negative displacements are greater. However, the asymmetry dependency on the diffusion time seems to follow a saturative profile as shown in the right side image of fig. 7.6.

The selected dynamic viscosity influences results. Depending on the considered type of tissue (dynamic viscosity) the results can be more or less relevant. Indeed, for a given temperature gradient, considering the dynamic viscosity of free water diffusion leads to a greater diffusivity gradient than that obtained considering white matter tissue (fig. 7.3 left). The considered tissue type therefore affects the resulting EAP as illustrated in fig. 7.7.

Conclusion

The designed theoretical phantom, admitted that the underlying assumptions are met, helps realizing an important limitation in the measurement of EAP asymmetry. Indeed, looking at the values of H in figs. 7.5 and 7.6 one can see that they are very low. This is even more relevant when considered that experiments do not involve any noise. Therefore the measurement of asymmetry, with the pre-sented phantom, would impose many challenges about the actual observability of the potential phenomenon.

Nevertheless, the developed phantom could be used as a tool to investigate asymmetric properties of the EAP, under the conditions specified in section 5.5, as a function of the controllable parameters of an experimental acquisition, and could give important insights.

This chapter concludes the theoretical investigation over the complex signal attenuation carried out also in chapters 5 and 6. These chapters contributed in clarifying the relationship between complex signal and EAP asymmetry, and in introducing new examples and quantification paradigms where asymmetric diffusion might occur. In the next chapter, the complex signal is exploited in a totally different manner, neglecting the possible contribution of diffusion to the imaginary part of the signal itself. The complex diffusion-weighted images (DWIs) will then be processed to perform denoising, thus overcoming problems that typically arise when using magnitude DWIs instead.

Overview

The chapter takes into account the complex nature of DWIs and exploits it to avoid computing the magnitude DWIs which are affected by a bias, e.g. the Rician bias, that is related to the non-Gaussianity of the noise distribution. This is achieved by means of phase correction. In this chapter the effects of phase correction on the signal bias reduction are studied. It is also proposed an analysis of the effects of phase correction in improving the scalar metrics obtained with signal representations such as DTI and MAP-MRI. After this, the strengths and limitations of the current methodological framework are revised. As a consequence, it is proposed an improved phase correction framework that overcomes these limitations. The result is a phase correction procedure that allows estimating dynamically the amount of correction required, and that takes into account the spatial variability of the noise, which is proper of DWIs. This is validated in silico with experiments based on a Human Connectome Project dataset, and in vivo on acquired data.

Keywords noise floor; non-Gaussianity; Rician bias; MAP-MRI; DTI;

Introduction

Diffusion-Weighted Magnetic Resonance Imaging (DW-MRI) is inherently a low Signal-to-Noise Ratio (SNR) technique [Jones and [START_REF] Derek | squashing peanuts and smashing pumpkins": How noise distorts diffusion-weighted mr data[END_REF]. More diffusion weighting -globally encoded by a larger b-value -leads to lower signal intensities and consequently to a poorer SNR. In such a low SNR regime, the magnitude of the complex DW signal, as function of q/b-value, can be dominated by a bias, namely noise floor, which is due to the non-Gaussian distribution of the noise, as described in section 8.2. The noise floor causes a signal overestimation, as shown in fig. 8.1, that is more important at high b-values and when diffusion is less restricted, i.e. when the signal decay is fast due to a higher apparent diffusivity, as explained in section 4.3. This introduces a bias that leads to the distortion of the estimated quantitative diffusion metrics, for instance those described in section 4.4. where Gaussian diffusion assumptions -or an exponential description as function of the b-value -might still well approximate the signal attenuation along a restricted direction. However, some diffusion MRI techniques require the acquisition of Diffusion-Weighted Images (DWIs) at relatively high b-values [START_REF] Özarslan | Mean apparent propagator (map) mri: a novel diffusion imaging method for mapping tissue microstructure[END_REF][START_REF] Daniel | A general framework for experiment design in diffusion mri and its application in measuring direct tissue-microstructure features[END_REF][START_REF] Zhang | NODDI: practical in vivo neurite orientation dispersion and density imaging of the human brain[END_REF], where the Noise Floor affects the signal estimation and consequent parameter calculations. Indeed, at higher b-values methods like MAP-MRI, described in section 4.4.3, render a more suitable signal representation since the Gaussian diffusion assumption is released. For this very reason, one interested in studying the diffusion signal with high diffusion weightings 1 should consider using this kind of signal representation. However, in this high diffusion-weighting regime the presence of noise floor is more accentuated, as shown in fig. 8.1, thus it is expected the scalar metrics computed from MAP-MRI, i.e. the RTOP, RTAP, RTPP, MSD, and QIV described in section 4.4.4, to be affected to a great extent. This chapter contributes to clarifying the effect of the noise distribution, and that of the noise floor, on the diffusion signal attenuation and on the scalar metrics computed from DTI and, for the first time, MAP-MRI representations (q-space metrics). This is achieved in different ways. From a signal point of view, it is proposed a simulation in which the complex signal for diffusion in the proximity infinite plate, seen in section 5.3, is reconstructed with SHORE from its magnitude, with the proposed maximum-likelihood estimation (MLE) in eq. ( 4.27), or directly 1 The MGH Human Connectome Project releases datasets with maximum b-value of 10000 s/mm 2 . in the complex domain. In a more realistic scenario, both simulations and real data will be used to see the effect of the non-Gaussian noise distribution on signal and metrics. At the same time, it will be presented a method to remove the noise bias from the data, namely phase correction (PC), and its efficacy will be characterized and quantified. The presented analysis is based on the work developed in Pizzolato et al. [2016d]. The PC method, originally proposed by [START_REF] Matthew A Bernstein | Improved detectability in low signal-to-noise ratio magnetic resonance images by means of a phase-corrected real reconstruction[END_REF], consists of a strategy for removing the noise floor from the magnitude DWIs by estimating the true phase of the corresponding complex DWIs. This chapter contributes by quantifying the effect of a total variation (TV) based PC method on DTI and q-space metrics and by showing the inherent limitations. In fact, PC is based on image filtering of the complex DWIs, but at the present time it is not indicated what amount of filtering, e.g. how much regularization, needs to be performed. In addition, the spatial variability of noise in the DWIs is not taken into account. The chapter also contributes by proposing a suitable strategy to overcome such limitations.

The next section describes the issue of the non-Gaussian noise distribution of magnitude DWIs. Section 8.3 will provide on overview of signal reconstruction based on magnitude and complex signal. Then, section 8.4 will introduces the PC method whose efficacy will be evaluated in many ways in the experimental part reported in section 8.5, where also limitations of the procedure will be highlighted. Section 8.6 presents the proposed improvements about PC. The proposed strategy is validated in silico and, in section 8.7, on real data.

The Rician noise distribution of MRI data

The term diffusion-weighted image (DWI) indicates the image obtained from the MRI scanner where the contrast, e.g. the structure of the imaged brain, is localized in the image within a frame of reference that can be related to that of the external world, i.e. the so-called image space. However, the image is actually acquired in the frequency space, the so called k-space, where each sample measured by a coil is independent from the others. The main source of random fluctuations in these collected samples is thermal noise, whose variance depends on the following relationship

σ 2 thermal ∝ 4k B T R ef f B W (8.1)
where as usual k B is Boltzmann's constant, T the absolute temperature of the resistor, R ef f is the resistance of the coil loaded by the scanned object, and B W is the bandwidth of the noise-voltage detecting system [START_REF] Aja | A review on statistical noise models for magnetic resonance imaging[END_REF]. Henkelman [1985] disclosed that under some assumptions, particularly the fact that the noise affects equally to all the frequencies, the k-space signal coming from a coil can be modeled as an additive complex white Gaussian noise

s(k) = a(k) + n(k; 0, σ 2 (k)) (8.2)
where a(k) is the underlying true signal, n( • ) indicates the additive noise, and s( • ) is the noisy signal. Looking at eq. ( 8.2) one can note that the noise variance varies with the frequency position. Later on, section 8.6 will describe a proposed strategy to include the eventual knowledge of the non-stationarity of noise within a phase correction procedure. However, if the noise in the signal is considered stationary then σ 2 (k) = σ 2 so that the noise can be considered as a complex bivariate Gaussian distribution [START_REF] Mark | Measurement of signal intensities in the presence of noise in mr images[END_REF] as

n(k; 0, σ 2 ) = n r (k; 0, σ 2 ) + j • n i (k; 0, σ 2 ) (8.3)
where the superscripts "r,i" represent the real and imaginary parts of the signal. The complex DWI, i.e. the image space signal that depends on the position p, is obtained via inverse Discrete Fourier Transform (iDFT) of s(k), considering that data is sampled on a Cartesian lattice. If no further processing is performed, the noise in the image is still Gaussian according to

DWI(p) = A(p) + η(p; 0, σ 2 (p)) (8.4)
where again σ 2 (p) = σ 2 in case of stationarity. The noise distribution is depicted in fig. 8.2. Furthermore, note that eq. ( 8.4) has a general meaning, and it is not just restricted to DWIs. The noise distribution of the magnitude image M = |DWI(p)| is non-Gaussian, and it depends on the number of coils used for the acquisition, the adopted reconstruction strategy, and on the acceleration strategy [START_REF] Aja | A review on statistical noise models for magnetic resonance imaging[END_REF]. For a single coil reconstruction the magnitude follows a Rician distribution

p(M | σ, |A|) = M σ 2 exp - M 2 + |A| 2 2σ 2 I 0 |A|M σ 2 Θ(M ) (8.5)
where I 0 ( • ) is the zeroth order modified Bessel function of the first kind, and Θ( • ) is the Heaviside step function. In eq. ( 8.5) the dependence on the spatial position p has been dropped for the sake of clarity. The review of the noise distributions of magnitude DWI obtained for different acquisition setups is beyond the scope of the present chapter. In fact, the phase correction method analyzed and extended in this chapter can be applied to complex DWIs coming from each single coil. However, for the sake of experiments presented in section 8.7, it is worth mentioning that when data is acquired with SENSE [Pruessmann et al., 1999], such as in the case of the real data used in section 8.7, the distribution of the magnitude image still follows a Rician distribution [START_REF] Aja | A review on statistical noise models for magnetic resonance imaging[END_REF].

The next section performs a comparison between a 1D signal reconstruction based on the signal magnitude, accounting for the distribution reported in eq. (8.5), and the reconstruction based on the corresponding complex signal.

In-voxel magnitude and complex reconstructions

Unlike the rest of this chapter, this section presents a comparison between magnitude and complex-based diffusion signal reconstruction within a voxel that is based on the work developed in [START_REF] Pizzolato | Magnitude and complex based diffusion signal reconstruction[END_REF]. This means that it is not assumed to process DWIs, but rather to look at the hypothetical complex signal attenuation measured in a single voxel. To do so, it is considered the signal attenuation predicted for diffusion in the proximity of an infinite plate given in eq. (5.14). Admitted that the noise distribution of the magnitude data points is known, a maximum likelihood estimation (MLE) could be employed to remove the noise floor bias, without the need of actually performing a reconstruction in the complex domain. The objective of the comparison is to study whether it is equivalent to consider a MLE reconstruction and a complex based one.

The magnitude attenuation is reconstructed from either the magnitude or complex samples of eq. (5.14), where Gaussian complex noise was added on the real and imaginary parts with σ = SNR -1 . The number of samples used is denoted as N s , where the number is doubled when complex samples are considered. The reconstruction is performed with SHORE using N = 6 bases to represent either the magnitude, the real, or the imaginary parts. In the first case the actual noise variance σ 2 was used to perform the MLE reconstruction indicated in eq. (4.27). In the second case the magnitude was reconstructed from the complex samples, according to eqs. (4.22) and (4.23). For the sake of illustration also a standard LLS-based reconstruction from the magnitude samples, i.e. eq. ( 4.24), was performed. The root sum of squares error of the residuals with respect to the noise-free magnitude attenuation was computed for the three techniques: RSS M for the LLS technique, RSS M L for MLE, and RSS C for the complex-based reconstruction. Figure 8.3a shows a peculiar example of the trend of the RSS values as function of the SNR. Note that while, in this case, RSS C remains lower, the MLE outperforms the standard LLS reconstruction from the magnitude samples only for SNR < 5. Indeed, for SNR ≥ 5 the Rician distribution well approximates a Gaussian [START_REF] Pajevic | Parametric and non-parametric statistical analysis of dt-mri data[END_REF]Basser, 2003, Salvador et al., 2005], thus not constituting a bias for LLS estimation2 . Figure 8.3b reports an analysis of the best reconstruction as function of the number of samples, along rows, and of the SNR, along columns. It reports the difference RSS M L -RSS C where red corresponds to positive values, thus a better performance of the complex reconstruction, and blue corresponds to negative values, thus a better performance of the MLE reconstruction from the magnitude samples. From the image it can be deduced that the reconstruction from the complex samples is generally recommended unless there are few samples (below 9). However, MLE requires the knowledge of the noise variance, which is actually difficult to estimate. Moreover, results might also reflect the fact that the magnitude estimated from the complex samples adds two error contributions from two different LLS fitting procedures, i.e. the real and imaginary fittings.

The example proposed in this section, although it is based on the complex attenuation obtained under the very specific assumption reported in section 5.3, clarifies what are the inherent issues in using the magnitude signal. In the rest of the chapter, the "in-voxel" representation will no longer be studied, since it relates to an ideal case in which the signal attenuation is only due to diffusion. Therefore, the rest of the chapter focuses on the characterization and processing of complex DWIs. Particularly, the next section describes phase correction, that is used to change the signal's noise distribution to restore Gaussianity.

Phase correction

Phase correction (PC) is a procedure initially proposed by [START_REF] Matthew A Bernstein | Improved detectability in low signal-to-noise ratio magnetic resonance images by means of a phase-corrected real reconstruction[END_REF] that processes a complex DWI to produce a real-valued DWI that contains the image content affected by additive Gaussian noise, as opposed to the magnitude DWI that is affected by Rician noise as seen in section 8.2. In DW-MRI the images are typically acquired slice by slice in 2D, thus the components of the position vector p of eq. (8.4) are x and y. Consequently, the complex DWI can be written as

DWI xy = rDWI xy + j • iDWI xy (8.6)
where r and i indicate the real and imaginary parts. PC aims at the estimation of the true, i.e. noise-free, phase of eq. ( 8.6) to be used to "correct" the image. In fact, if DWI xy is a good estimation of the phase, then the phase-corrected image is obtained via complex rotation ) where the Rician bias, e.g. the noise floor, is absent. The effectiveness of phase correction clearly depends on the quality of the phase estimation. Recent phase corrections for noise floor removal consist on filtering the real and imaginary images, i.e. the rDWI and iDWI, to obtain a low-frequency version of the DWI's phase, which is used to complex-rotate the rDWI and iDWI such that the former contains signal plus Gaussian distributed noise, and the latter only noise (which will be discarded). The filtering is typically performed via a convolution procedure as in [START_REF] Douglas E Prah | A simple method for rectified noise floor suppression: Phase-corrected real data reconstruction with application to diffusion-weighted imaging[END_REF], [START_REF] Sprenger | Real valued diffusion-weighted imaging using decorrelated phase filtering[END_REF] or via a non-linear procedure as in the work of [START_REF] Karla | Nonlinear phase correction for navigated diffusion imaging[END_REF], [START_REF] Eichner | Real diffusion-weighted MRI enabling true signal averaging and increased diffusion contrast[END_REF]. However, the correct estimation of the low-frequency phase depends on the correct choice of the convolution kernel (and its size) or regularization parameter. Therefore, the effectiveness of phase correction on signal debiasing and diffusion parameters estimation, such as DTI and MAP-MRI metrics, needs to be assessed. In this chapter it is implemented a total variation method, known to better preserve discontinuities in the images [START_REF] Eichner | Real diffusion-weighted MRI enabling true signal averaging and increased diffusion contrast[END_REF]. Particularly, for each image u 0 defined on coordinates x ∈ X, y ∈ Y , the filtered image u is found as the minimizer of

inf u∈C λ X,Y (u 0 -u) 2 dxdy + X,Y |∇u|dxdy (8.8)
where λ is the regularization parameter that expresses the attachment to the data. In this thesis, eq. (8.8) has been developed in house by adopting a numerical implicit scheme that does not require any additional parameter and that is suitable for a Gauss-Seidel iteration. The estimated real and imaginary images are then given by rDWI xy = (u)

iDWI xy = (u) (8.9)
and the estimated phase is

DWI xy = tan -1 iDWI rDWI xy (8.10)
that is used to perform the complex rotation in eq. (8.7). The following section presents the experimental part of the chapter that quantifies the effect of phase correction on the signal and on the scalar metrics calculated from DTI and MAP-MRI. The analysis also contributes by highlighting that the selection of the regularization parameter λ, not solved in literature, is paramount to benefit from phase correction. This will be solved later on in section 8.6.

Characterization of phase correction

The present section contributes with the characterization of the effects of phase correction in terms of signal debiasing and, newly, on DTI and MAP-MRI (q-space) metrics. The analysis is based on the work developed in Pizzolato et al. [2016d].

After a methodological part, the analysis is divided in three main experiments grouped in two parts, depending on whether the influence of phase correction is quantified on signal or on scalar metrics. The performed in silico experiments are based on synthetic phantoms and on HCP3 data.

Simulation and diffusion signal reconstruction

The complex DWIs have in all cases been created by generating a synthetic phase image, Φ xy , associated with a magnitude image, M xy . The phase images are created in order to mimic the outcome of subject movements. indeed, the phase image is created as a bidimensional sinusoidal wave oriented along the direction v = (v x , v y ), with frequencies f x , f y , and initial shifts φ x , φ y

Φ(x, y) = π • sin 2π v x ||v|| f x x w x + φ x + 2π v y ||v|| f y y w y + φ y (8.11)
where w x , w y are scale parameters: in this case they correspond to the width of the image along the corresponding direction, i.e w x = card(X) and w y = card(Y ). Eventually, constant phase patches are added. Assuming to have the ground-truth images of magnitude M xy and phase Φ xy , the latter resulting from eq. (8.11), then

rDWI xy = M xy • cos(Φ xy ) + η r xy iDWI xy = M xy • sin(Φ xy ) + η i xy and η r xy , η i xy ∈ N (0, σ 2 ). (8.12)
The noise is added with a value of σ calculated according to the DW-MRI convention

σ = x,y ρ(x, y)M b=0 xy card[ρ(X × Y )]SNR 0 (8.13)
where SNR 0 is defined on the magnitude image without diffusion weighting M b=0 xy , and ρ ∈ {0, 1} is a mask defined on the pairs (x, y), e.g. a mask of the tissue-related signal like the brain mask. The Rician magnitude |DWI| xy and the phase DWI xy are calculated from the real and imaginary parts in eq. (8.12).

The data used for the experiments is a HCP brain dataset corrected for eddy currents. From this, DWIs of interest are selected corresponding to b ∈ [0, 1000, 3000] s/mm 2 , i.e. the DWIs are grouped based on their diffusion-weighting. Other experiments use Phantomαs [START_REF] Caruyer | Phantomas: a flexible software library to simulate diffusion mr phantoms[END_REF] to obtain the ground-truth magnitude images, M xy . This software requires input with a geometrical description of tissue structures and fiber bundles. In this chapter it is used the well known geometry produced for the HARDI reconstruction challenge 2013 4 . Here, the DWIs are generated for a 3-shells scheme with b ∈ {1000, 2000, 3000} s/mm 2 , 51 samples per shell, with samples uniformly distributed within and among shells according to the electrostatic repulsion strategy developed by [START_REF] Caruyer | Design of multishell sampling schemes with uniform coverage in diffusion mri[END_REF].

The DTI and MAP-MRI signal reconstructions, on which the calculated scalar metrics are based, have been performed as indicated in sections 4.4.1 and 4.4.3. In particular the DTI fitting has been performed with eq. ( 4.20), constraining the signal to be positive, and the MAP-MRI reconstruction has been performed with analytic Laplacian regularization [Fick et al., 2016] combined with positivity constraints on the recovered EAP [ Özarslan et al., 2013b]. Indeed, the phase-corrected real DWIs can contain negative values because the final expected noise distribution has zero-mean and is Gaussian, and the noise floor in the signal should be absent. Therefore, positivity constraints on the signal or EAP are required to obtain a physically meaningful reconstruction.

Influence on diffusion-weighted signal

This first part of the analysis presents two experiments with the objective of quantifying the effect of phase correction on signal debiasing. This is achieved via synthetic experiments that are created by processing real data from a HCP dataset.

First experiment: effect on signal noise distribution

The goal of this experiment is to quantify in vivo the actual DWIs intensities and to see the effect of a given SNR 0 on the actual SNR and, consequently, on the noise distribution. At the same time, the debiasing effect of phase correction is visualized.

To do so, a HCP dataset is clustered to obtain typical signal values at b-value 1000 and 3000 s/mm 2 . For each b-value, the pixels of the corresponding DWIsincluding all the directions -are clustered with k-means to divide the signal intensities into 4 clusters. The centroid of each cluster is used to define background, low, medium, and high mean signal values respectively. Based on these, it is created a ground-truth synthetic magnitude image -M xy in eq. (8.12) -composed of three circles each containing, from left to right, low, medium and high signal respectively, as shown in fig. 8.5. Outside the circles, background signal was added. A synthetic phase was generated and the noisy complex DWI consequently created. After calculating the average b = 0 signal (S(0) avg = 758 a.u.) in the HCP dataset, noise was added as in eq. (8.12) in low SNR regime: SNR 0 = 10. Figure 8.5 shows, for each b-value, the noisy magnitude |DWI| xy and the estimated phase-corrected real part (DWI pc xy ) (λ set to 0.75 after visual inspection). In addition, an effective SNR map is present along with histograms of the magnitude and phase-corrected real signals for each circle.

It is concluded that in both cases the phase-corrected real image presents more contrast with the background compared to the magnitude one. This is more evident at low SNR values -left circle at b = 1000 s/mm 2 , left and central circles at b = 3000 s/mm 2 -that are more likely with high b-values. The (DWI pc xy ) shows darker colors, i.e. lower signal intensities, as highlighted by the histograms: the magnitude (green line) has a Rician distribution for low SNRs (typically below SNR = 5) whereas the estimated real part (blue line) always shows a Gaussian distribution, thus including negative signal intensities. It is pointed out that since this is an experiment grounded on real data, the centroid of the clusters -especially at low signal values -are based on Rician data and might overestimate the actual (noise-free) ones. This means that the Rician bias in histograms (green line) might be an underestimation of the true one. The next experiment will perform a less classical comparison between magnitude and phase-corrected data, and aims at the prediction of brain areas that will be most likely affected by noise floor -or more generally by Rician bias -based on statistical metrics.

Second experiment: bias prediction in the brain

This experiment aims at the prediction and visualization of the Rician bias, i.e. the distance from Gaussianity, and of the effect of phase correction in a brain. The synthetic data is based on a HCP dataset with the objective of creating a mean ground-truth magnitude DWI, M xy .

Particularly, the mean b = 0 image S(0) xy is calculated, for a slice of interest, by averaging the 40 non-diffusion-weighted images in the dataset. Since at b = 0 the SNR is very high the averaging procedure is not biased, i.e. the noise distribution can be considered Gaussian. Then, all the DWIs corresponding to b = 1000 s/mm 2 are selected to perform DTI reconstruction to obtain the mean diffusivity map, MD xy , where for each voxel MD is calculated according to eq. (4.34). At this point, the ground-truth magnitude DWI is obtained at any b-value by extrapolating with (8.14)

Although eq. (8.14) assumes Gaussian isotropic diffusion, this phantom represents an average description of a magnitude DWI. After generating a synthetic phase image, as described in section 8.5.1, the noisy complex DWI is calculated, for each b ∈ {1000, 2000, 3000} s/mm 2 , as in eq. (8.12). Figure 8.6 shows the b = 0 magnitude and the phase used for the phantom (left column). Then, the Rician magnitude |DWI| xy and the phase-corrected real part (DWI pc xy ) (λ = 0.75) are calculated, where noise was added with SNR 0 = 10. In order to compare the bias of the two, it is necessary to calculate an unbiased version of the magnitude DWI, that is affected by Gaussian noise. For this reason, an additional magnitude image, |DWI| G xy , is created by adding Gaussian noise to M xy with the same SNR 0 . This image constitutes a reference for Gaussianity. A total of 1000 noise occurrences are generated and, for each pixel of the images, the signal intensities histogram (as in fig. 8 xy (which is used as reference). The hypothesis is that, for each pixel, the phase-corrected signal distribution should be closer to that of the reference Gaussian magnitude image, than the Rician magnitude. The discrete Hellinger distance [START_REF] Hellinger | Neue begründung der theorie quadratischer formen von unendlichvielen veränderlichen[END_REF] bias occurs in a larger number of pixels. Conversely, the phase-corrected real part does not show a clear change and the amount of bias remains steady as the b-value increases.

Influence on DTI and q-space metrics

This experiments quantify the bias on the estimated DTI and MAP-MRI metrics for the Rician magnitude, as well as the debiasing power of phase correction, by looking at the change in the distributions of such metrics compared to the Gaussian noise case. In other words, the aim is to quantify the distance from Gaussianity of the metrics, and the effect of phase correction on them. Moreover, the experiment investigates the effect of the regularization parameter λ of the total variation filtering in eq. (8.8) on which phase correction is based. Since it is required to know the ground-truth values of metrics, the experiment has been performed considering synthetic complex DWIs whose magnitude images are generated with Phantomαs [START_REF] Caruyer | Phantomas: a flexible software library to simulate diffusion mr phantoms[END_REF], as described in section 8.5.1.

In order to create complex DWIs the tool has been extended to include phase images.

Particularly, when creating the DWI corresponding to the gradient direction g = (g x , g y , g z ), the synthetic phase image, computed with eq. (8.11), is oriented towards v = (g x , g y ). Different slices in the volume, i.e. the DWIs stacked along the z axis, have a phase shift between them to mimic a movement phase artifact. Figure 8.7 shows the original noisy data (Rician magnitude, phase, noisy real and imaginary parts) and the one after phase correction, for a reference slice (b = 1000 s/mm 2 ).

The experiment takes into account different values of SNR 0 and different acquisition schemes, and considers all of the metrics calculated with DTI and MAP-MRI reconstructions. A summary is provided in table 8.1. For each SNR 0 ∈ {10, 20, 30}, the Rician magnitude data and the Gaussian reference DWI image are generated. The phase-corrected real DWI is also computed. However, this experiment also investigates the effect of the regularization parameter λ of the total variation filtering in eq. (8.8). Therefore, for each SNR 0 six phase-corrected datasets, with λ ∈ {0.25, 0.5, 0.75, 1, 2, 5}, are calculated. For each combination of SNR 0 and type of data -Rician, Gaussian, and the six phase-corrected ones -the acquisitions based on single-shell scheme (at b-value 1000, 2000 and 3000 s/mm 2 ) are reconstructed with DTI, whereas multi-shell data, i.e. the aggregation of the three single-shell schemes, is reconstructed with both DTI and MAP-MRI. Calculated metrics based on DTI are: mean diffusivity (MD), principal diffusivity (PD), and fractional anisotropy (FA). The q-space metrics based on closed formulas derived on MAP-MRI are: the return to origin (RTOP), axis (RTAP), and plane (RTPP) probabilities, the mean squared displacement (MSD), and the q-space inverse variance (QIV).

Metrics distributions are calculated within a mask generated on the groundtruth synthetic dataset. Particularly the mask is generated by considering only the voxels where RTOP ∈ [0.5e6, 0.7e6]. The range was chosen based on visual inspection, since it allows to extract a mask of fibers, as shown in fig. 8.8. DTI and q-space metric histograms based on the Rician magnitude and on the reference Gaussian data are generated for each SN R 0 considering voxels within the mask. Similarly, phase-corrected metric histograms are generated for each considered value of λ. Note that λ has a great influence on the phase correction results. Particularly, a large λ implies strong attachment to data, resulting in a poor phase correction since the estimated low-frequency phase is very similar to the original noisy one, which means that DW I xy -DW I xy ≈ 0 in eq. (8.7). Indeed, the blue histograms (phase-corrected) in the first row of fig. 8.9 almost entirely overlap the green ones (Rician magnitude data). As the attachment to data decreases (from top to bottom), the blue phase-corrected histograms move towards the (red) Gaussian based distributions, visually reducing the distance from Gaussianity. As in second experiment, the distance from Gaussian metrics is quantified using Hellinger's formula in eq. (8.15), where histograms this time relate to metrics and not to signal intensities. Figure 8.10 illustrates the variation of the H distance for the phase-corrected data as function of λ, for each acquisition setup, reconstruction method (DTI, MAP), and diffusion metric. In each image, the dashed lines represent the distance of the metric calculated on Rician magnitude data from the corresponding Gaussian one, whereas the solid lines report the distance for metrics calculated on phase-corrected data, which varies with λ (abscissa). Color codes indicate the SNR 0 value. It can be observed that phase correction leads to metric distributions that are closer to the Gaussianity (H distance close to 0) than the Rician magnitude ones, for specific ranges of λ. In general, phase correction debiases the metric distributions up to a great extent. The improvement over the Rician magnitude is clearly correlated with the combination of acquisition scheme -especially the maximum b-value -and SNR 0 as also indicated by the signal intensities experiments illustrated in figs. 8.5 and 8.6. Indeed, at high b-values the signal is low -especially along the less restricted diffusion direction -which, in combination with a poor SNR 0 , causes the effective SNR to fall well below 5 where a Rician distribution diverges from a Gaussian one. Therefore, the best value of λ (highlighted with a dot in fig. 8.10) also depends on these factors. The best λ seems depend on the considered metric. For instance, at SNR 0 = 30 the best λ for RTPP is 0.75 whereas for RTAP is 2. This can be associated to the fact that metrics that are highly related to signal measured along the less restricted diffusion direction, i.e. low intensity signal, such as PD and RTPP, benefit more than others of phase correction. A quantitative summary of the results is available in table 8.2.

The influence of the phase correction amount

Phase correction does not always lead to bias reduction

Looking at fig. 8.10 it is particularly interesting to note that in some cases, as for DTI at b = 1000 s/mm 2 , too much filtering (small λ) causes the phase-corrected metric distributions to be more distant from Gaussianity compared to those based on the Rician magnitude (dashed lines). In fact, the problem of the choice of λ is not solved in the literature, and empiric values based on visual inspection are normally used. Moreover, the present method does not take into account the spatial variability of the noise variance. Therefore, the next section proposes an improvement over this situation by disclosing a method that computes automatically the correct value of λ and takes into account spatial variability of noise.

Dynamic phase correction with noise variance map

The experimental part in section 8.5.3, and especially fig. 8.10, discloses that phase correction has a great debiasing potential (see table 8.2) but at the same time it is not useful unless the estimated phase is optimal. In other words, phase correction is not effective -and can lead to more bias than that of the magnitude DWIif the regularization is not properly performed. In fact, leaving the choice of λ to empiricism, i.e. visual inspection, can lead to altered image content, and is inefficient. Moreover, the phase correction strategy based on the total variation functional in eq. (8.8) does not take into account the spatial variability of the noise variance as expressed in eq. ( 8.4), that is when σ 2 xy = σ 2 (p) = σ 2 . This section contributes by proposing an improvement over the current technique, to take into account both an automatic choice of the regularization parameter, which is dynamically updated while estimating the phase, and the spatial variability of the noise variance.

Dynamic TV regularization

The optimal phase correction is the one that allows estimating a phase image, DWI xy , that after performing the complex rotation in eq. (8.7) leads to a phasecorrected imaginary image, {DWI pc xy }, that only contains Gaussian noise. The phase image is calculated from a denoised version of the original complex DWI xy via the total variation (TV) approach reported in eq. (8.8). It is here remarked that since the TV filtering is applied in the complex domain, then additive white Gaussian noise can be considered on the real and imaginary parts, according to eq. (8.3). Therefore, if an estimation of the noise variance, σ2 , is available, the discrepancy criterion predicts that the optimal denoising satisfies

i -i 0 2 2 = N 2 σ2 (8.16)
where i is the denoised image, i 0 the original one, and N 2 is the total number of pixels. In order to find λ satisfying eq. (8.16) [START_REF] Chambolle | An algorithm for total variation minimization and applications[END_REF] proposed an iterative rule which is practically implemented as

λ n = λ n-1 i -i 0 2 N σ (8.17)
where n is the TV iteration and where λ 0 = 0.7079 σ + 0.6849 σ2 (8.18) as suggested by [START_REF] Getreuer | Rudin-osher-fatemi total variation denoising using split bregman[END_REF]. As suggested in the same reference, the λ obtained after one TV minimization, i.e. after one complete cycle of iterations, can be used as input for the next one. This leads to monotonic convergence of λ n as function of cycles, as shown in the left side example in fig. 8.12. However, it is here noted that in practice the rule in eq. ( 8.17) can be applied to iterations rather than cycles, meaning that only one TV cycle can be performed. A similar consideration was also made by [START_REF] Chambolle | An algorithm for total variation minimization and applications[END_REF]. Although also in this case convergence is reached, it is here shown that monotonicity is not guaranteed, as reported by the right plot of fig. 8.12. On the other hand, this strategy leads to fast convergence. The two strategies lead to very similar results. In both cases, the rule in eq. (8.17) was considered at each iteration/cycle separately for the real and imaginary parts, averaging the results for λ n . Experiments proving the effectiveness of the proposed solution follow.

Synthetic validation

To show the performance of the dynamic phase correction, synthetic experiments are performed based on the same phantom adopted in the second experiment of section 8.5.

Spatial variability of the noise variance

In order to consider the spatial variability of the noise variance, eq. (8.8) is adapted as where w(x, y) is a weighting function that can be obtained with a MRI noise map.

In particular, the weighting function corresponds to the variance map normalized such that the sum of its elements is equal to the number of pixels in the image, N 2 , or in the selected mask. In other words, the weighting function is opportunely scaled such that it does not add energy to the system. The solution in eq. ( 8.20) can be plugged together with the dynamic selection of the regularization parameter λ seen in section 8.6.1. Particularly, supposing to have the knowledge of the mean variance in the image σ2 , this can be used in eq. (8.17) while using the map w(x, y) at the same time. Some experiments, that show the difference between considering the noise stationarity or varying, follow.

Synthetic validation

In this experiment a synthetic noise standard deviation map was created with average standard deviation σ = 50.46 and maximum standard deviation σ max = 70.58. The map was obtained by applying a radial gradient originating from the center, and is shown in the left image of fig. 8.13. The variability of the standard deviation is better appreciated in the central image of fig. 8.13 which reports the standard deviation map subtracted by the average. Pixels where the standard deviation is larger than the average are reported in white colors, whereas the opposite situation leads to blue colors. Both stationary and spatially varying phase corrections were applied to the data, feeding σ = 50.46 to both of the procedures. For the spatially varying phase correction, additionally, the weithing map w(x, y) was calculated based on the standard deviation map σ(x, y) shown in fig. 8.13.

At convergence, both of the phase correction procedures could estimate the average standard deviation via eq. ( 8.19), σP C = 50.9 and σwP C = 51.6. Additionally, when using the noise map it was possible to estimate the maximum standard deviation σwP C max = 71.98. However, as expected, the stationary phase correction left unchanged the spatial variability of the noise standard deviation, whereas by using the proposed approach it was possible to render the noise variance stationary. Indeed, the right image in fig. 8.13 shows the difference between the locally estimated standard deviation on the phase corrected imaginary images without taking into account spatial variability and with the proposed approach. The estimation was performed with a convolutional filter of radius 10. It can be seen that this right image mimics the one in the center. In the following section, the dynamic phase correction with spatial noise variability is validated in vivo.

Real data validation

This section reports the application of the proposed dynamic phase correction with spatial noise variability on real data. The adopted dataset5 is composed by 10 shells, b ∈ {56, 222,501,890,1390,2002,2725,3560, 4505, 5562} s/mm 2 antipodically acquired along 7 collinear directions, and with 7 b 0 images. A DWI of the dataset is reported in fig. 8.14 for illustration. The acquisition was performed with SENSE [START_REF] Klaas P Pruessmann | Sense: sensitivity encoding for fast mri[END_REF], and the sequence parameters were G max ≈ 44 mT /m, ∆ = 50 ms, and δ = 32.8 ms. The dataset includes one MRI noise map, on which it is possible to estimate the noise variance locally. Due to the low angular resolution, only DTI signal reconstruction was performed, obtaining the PD, MD, FA and color FA maps. In particular, the reconstruction was performed for magnitude DWIs and for dynamically phase-corrected real DWIs with and without considering the spatial variability of noise.

The complex noise map is reported in fig. 8.15, which also reports -in the right image -the calculated variance map to be used to take into account the spatial variation of the variance itself. This variance map is generated by locally estimating the standard deviation in the real and imaginary parts of the noise map with a convolutional standard deviation estimator of radius 10. The two maps are averaged and an average standard deviation is calculated ( σ ≈ 181) to be fed to the phase correction procedures. Then the two maps are used to generate the unique variance map. This is used to compute the map w(x, y) as explained in section 8.6.2, which is used in eq. ( 8.20). In the second part of the validation, the magnitude and phase corrected DWIs -with and without noise variabilityare reconstructed with DTI, as seen for the synthetic experiments. The objective is to characterize the effect of the different proposed versions of dynamic phase correction on the scalar metrics. Images in fig. 8.16 report the maps calculated for the three types of data: the Rician magnitude, and the two phase-corrected DWIs datasets. It appears clear the underestimation of FA, PD, and MD when DTI is performed on magnitude DWIs, as observed previously in the literature [START_REF] Eichner | Real diffusion-weighted MRI enabling true signal averaging and increased diffusion contrast[END_REF]. Indeed, when using magnitude DWIs the noise floor leads to overestimation of the signal, especially along the less restricted direction [Jones and [START_REF] Derek | squashing peanuts and smashing pumpkins": How noise distorts diffusion-weighted mr data[END_REF]. PD among all seems to be the most affected as seen in the region of the ventricles. This could be explained by the fact that PD expresses the value of the apparent diffusivity precisely along the less restricted direction. However, this has an influence on the estimation of the DTI tensor, and consequently on all the other metrics. In fact, images in fig. 8.17 show the difference maps between metrics calculated on the magnitude DWIs and those calculated on the dynamically phase corrected DWIs while accounting for noise variability. These difference maps confirm that without correction the values of the metrics are smaller.

In order to compare the results based on dynamic phase correction with (wPC) and without (PC) using the information about the noise variance variability, difference maps have been calculated in fig. 8.18. These maps show differences that are consistent with the noise maps in fig. 8.15. In particular, compared to wPC, PC renders larger values (blue) close to the center of the image and lower values at the borders (white).

Discussion

The proposed phase correction method could significantly affect the obtained diffusion metrics by reducing the noise bias. However, phase correction is not the only technique to remove bias from DWIs. Indeed, Rician maximum likelihood methods are available. However, an advantage of phase correction is that it can operate a filtering procedure in the complex domain, where noise is Gaussian. Therefore, there is no need of making assumptions about the similarity of pixels within the kernel region, or to use non-local means techniques [START_REF] Wiest-Daesslé | Rician noise removal by non-local means filtering for low signal-to-noise ratio mri: applications to dt-mri[END_REF][START_REF] Descoteaux | Impact of rician adapted non-local means filtering on hardi[END_REF]. To this regard it should be pointed out that phase correction does not aim at denoising but rather at bias correction. Indeed, the corrected real DWIs are still affected by (Gaussian) noise which can be removed via a subsequent processing. However, its use as a post-processing tool -that operates on the multi-coil reconstructed DWIs -should be compared in the future with the aforementioned Rician magnitude-based methods. On the other hand, phase correction has the advantage of being applicable before multi-coil reconstruction of the DWIs, that is exploiting the complex signal coming from each coil independently, which opens for a more general usage that is insensitive to the noise transformation operated by the chosen multi-coil reconstruction.

Phase correction relies on a filtering procedure which is a delicate choice. However, in this chapter it was used a total variation functional, as initially proposed by [START_REF] Eichner | Real diffusion-weighted MRI enabling true signal averaging and increased diffusion contrast[END_REF], which seems a good choice for preserving discontinuities in the image, that a conventional Laplacian filtering might destroy. However, a criteria for choosing the amount of phase correction was necessary to avoid possible mistakes due to an empirical choice. The proposed approach solves this problem and additionally allows accounting for the spatial distribution of the noise variance, which is important in all MRI modalities. Accounting for the noise spatial variability can be very important also when the phase correction is performed on single images acquired from each coil, before multi-coil reconstruction.

Results in section 8.7 show that phase correction has a great impact on the calculated metrics, where differences compared to magnitude-based estimates reach a large extent. As a consequence it is expected that the use of the proposed dynamic phase correction, possibly taking into account the noise spatial variability, can have a major impact on clinical studies and on diffusion post-processing. For instance, in tractography a more correct value of PD, used to infer the principal propagation direction of a tract, can play a game-changing role.

Conclusion

This chapter contributes by making an analysis of role of the complex diffusion MRI signal in denoising, with a particular focus on phase correction. The chapter proposes a study on the effects of phase correction on removing the Rician bias and on the computed DTI and q-space (MAP-MRI) metrics. It is shown that phase correction can be very effective in debiasing signal and metrics but only in the case the regularization, on which it is based, is properly performed. It is then proposed a strategy to improve the regularization by finding an optimal regularization parameter, i.e. the dynamic phase correction. Moreover, since noise in DWIs is generally non-stationary, it is also proposed to take into account the noise variability. The proposed approach is then validated on synthetic and real data, showing that it leads to large differences in the estimated metrics compared to results based on magnitude DWIs, differences that are supposingly justified by a reduced or absent signal bias. Chapter 9

Perfusion basics and topics

133

Overview

What is the characterization of perfusion in vivo? The early studies on circulation in organs by [START_REF] Stewart | Researches on the circulation time in organs and on the influences which affect it[END_REF] and [START_REF] Wv F Hamilton | Studies on the circulation[END_REF] were later formalized by [START_REF] Meier | On the theory of the indicator-dilution method for measurement of blood flow and volume[END_REF] in the indicator-dilution theory. This theory describes how to measure blood flow and volume of a specific organ by injecting a tracer -the indicator -in the vascular system and then measuring its dilution below the organ itself. This elegant theory was transposed to measure perfusion with DSC-MRI, and obtain maps of the cerebral blood flow, volume, and mean transit time. As the chapter reads, it appears clear the central role of the so-called transport function, and especially of its alter ego, the residue function, r(t), which describes the residual fraction of tracer in the voxel during time, a tissue-specific characteristic. After, the chapter focuses on the open challenges to be faced in the estimation of the residue function in vivo, namely time-delay and dispersion, challenges that will be mathematically formalized and tackled in the next chapters.

Keywords residue function; response function; delay; dispersion.

The aim of perfusion imaging in the brain

Perfusion imaging via Dynamic Susceptibility Contrast MRI (DSC-MRI) provides useful information in neurological diagnosis, characterization, and grading of brain tumors (such as gliomas), and ischemic strokes [START_REF] Shiroishi | Principles of T 2 *-weighted dynamic susceptibility contrast MRI technique in brain tumor imaging: Principles of T 2 *-Weighted DSC-MRI[END_REF].

This information is supplied to physicians in the form of hemodynamic/perfusion parameters brain maps, such as the cerebral blood flow (CBF), blood volume (CBV), and mean transit time (MTT). The parameters can be calculated by processing the concentration time-curves, obtained with eq. (3.26), for each voxel according to a well-known set of formulas that are part of the indicator-dilution theory pioneered by [START_REF] Meier | On the theory of the indicator-dilution method for measurement of blood flow and volume[END_REF], and that will be discussed later in more detail. However, for each voxel the processing normally considers a measured arterial and tissue concentration time-curves, C a (t) and C ts (t), which result from the intravenous injection of a tracer. The estimation of the perfusion parameters is closely related to the estimation of the time-dependent residual amount of tracer -the residue function r(t) -in the voxel. The main goal, in DSC-MRI, is therefore the knowledge of the full set of unknowns {CBF, CBV, MTT, r(t)}. However, as it will be more clear after, a complete understanding of the tissue perfusion requires the extension of this set of unknowns. Issues like the time-delay, τ , between the measured arterial and tissue concentration, and the dispersion of the actual arterial concentration pose severe questions in perfusion imaging that is the purpose of this thesis to solve.

Figure 9.1: The vascular bed, with injection point "P" and sampling point "Q", selected from the original drawing of [START_REF] Meier | On the theory of the indicator-dilution method for measurement of blood flow and volume[END_REF].

Indicator-dilution theory

The indicator-dilution theory, formalized by [START_REF] Meier | On the theory of the indicator-dilution method for measurement of blood flow and volume[END_REF] and later on applied to perfusion MRI, defines the relationship between important vascularrelated parameters namely the volume V (ml) of indicator -the tracer -in the vascular system, the flow F (ml/s) and the average time < t > (s) it takes for a particle of indicator entering the system to leave it definitely. The theory relies on two basic assumptions: the first, called representative behavior assumption, states that the flow of indicator particles is representative of that of the blood, and the second, called stationarity of flow, states that the distribution of the transit times of the particles through the system does not change with time.

Imagine a vascular bed such as that drawn by [START_REF] Meier | On the theory of the indicator-dilution method for measurement of blood flow and volume[END_REF] and reported in fig. 9.1, where an artery transports blood to a network of capillaries that, in turn, reunify into a venous vessel. The study of the perfusion of this vascular bed -which is assiociated to an organ of interest or the tissue's parenchymacan be carried out by injecting a concentrated bolus1 of indicator at the inlet "P" and measuring its concentration during time at the outlet "Q". If C(t) (mmol/L) is the concentration of indicator measured at the system's outlet, then after an instantaneous injection of q 0 units of indicator the rate at which it leaves the system at time t is F • C(t) (units/s), therefore the fraction of indicator leaving the system per unit time, known as the transport function, is given by [START_REF] Meier | On the theory of the indicator-dilution method for measurement of blood flow and volume[END_REF]

h(t) = F • C(t) q 0 (9.1)
which can be regarded as the probability density function of the transit times t of the system assuming

∞ 0 h(t)dt = 1 (9.2)
which implies that all the fluid entering the system will leave. Clearly the mean transit time can be computed as the average

< t >= ∞ 0 t • h(t)dt. (9.3)
Since the rate at which the fluid enters and leaves the system is F , then the volume of fluid leaving in a time dt is given by F • dt. However, at a specific time t only the fraction of particles that have a transit time equal to t, that is h(t), will leave the system and contribute to the concentration C(t) measured at that instant. Therefore, the rate per unit time, or flux, at which the particles having a transit time t will leave the system is given by F • h(t) • dt, and the volume of such particles is obtained by multiplying this rate by the transit time t required for them to be cleared out of the system dV = t • F • h(t) • dt. Hence, the total volume of the indicator in the fluid (blood) is the sum of all of the volumes for each possible transit time, that is [START_REF] Meier | On the theory of the indicator-dilution method for measurement of blood flow and volume[END_REF]]

V = ∞ 0 t • F • h(t)dt = F ∞ 0 t • h(t)dt (9.4)
where by using eq. ( 9.3), the integral can be substituted giving the relationship

V = F • < t > (9.5)
which is the result of the Central Volume Theorem found by [START_REF] Stewart | Researches on the circulation time in organs and on the influences which affect it[END_REF] and discussed by [START_REF] Wv F Hamilton | Studies on the circulation[END_REF].

The transport function is the frequency density function of transit times. Hence, the cumulative frequency function H(t) can be defined as

H(t) = t 0 h(θ)dθ (9.6)
which is dimensionless and where lim t→∞

H(t) = 1. (9.7)
The total fraction of tracer that has left the system from time zero up to time t is H(t), therefore the fraction of the amount of tracer remaining in the tissue is

r(t) = 1 -H(t) (9.8)
where r(t) is the residue function2 and, since H(0) = 0, r(0) = 1. Given its nature, at any time t the residue function r(t) can only assume positive or zero values, therefore non-negativity is a desired property while recovering it, according to the theory explained in the next section.

Perfusion convolution in the brain

This section describes the relationships between the perfusion parameters CBV, CBF and MTT in the brain, and the role of the residue function in their estimation. Particularly, in perfusion MRI we consider as known (with some limitations that will be clearer later) the input concentration to the vascular system, i.e. the arterial concentration C a (t), and the output, i.e. the tissue concentration C ts (t). As anticipated by [START_REF] Kenneth | Equations for measuring blood flow by external monitoring of radioisotopes[END_REF], C ts (t) is expressed as the convolution between C a (t) and the unknown residue function r(t) multiplied by the cerebral blood flow

C ts (t) = ρ k H CBF t 0 C a (θ)r(t -θ)dθ (9.9)
where ρ is the brain tissue density that is needed to provide the correct flow units, and k H accounts for the difference in hematocrit between large vessels H LV such as arteries, and small vessels H SV such as capillaries. In fact, the indicator is linked to the plasma (usually to the albumine) and, especially in small vessels, the flow of plasma is smaller than that of heritrocites which are the particles of interest for transportation of oxygen. Therefore, to obtain volume and flow values related to the whole blood, hematocrit should be taken into account as suggested by [START_REF] Meier | On the theory of the indicator-dilution method for measurement of blood flow and volume[END_REF], according to [START_REF] Katrin A Rempp | Quantification of regional cerebral blood flow and volume with dynamic susceptibility contrast-enhanced mr imaging[END_REF][START_REF] Østergaard | High resolution measurement of cerebral blood flow using intravascular tracer bolus passages. part i: Mathematical approach and statistical analysis[END_REF][START_REF] Knutsson | Absolute quantification of perfusion using dynamic susceptibility contrast MRI: pitfalls and possibilities[END_REF]]

k H = 1 -H LV 1 -H SV .
(9.10) Equation 9.9 expresses the aim of interpreting the arterial input function (AIF), C a (t), as a superposition of consecutive ideal boluses C a (θ)dθ injected at time θ, where for each ideal bolus the tissue concentration, C ts (t), at time t is proportional to C a (θ)r(t -θ) [START_REF] Calamante | Measuring cerebral blood flow using magnetic resonance imaging techniques[END_REF].

In order to estimate the CBF it is necessary to perform a deconvolution to obtain the response function R(t) = CBF • r(t) and then CBF = R(0) since r(0) = 1. However, this estimation can be corrupted by practical issues when processing DSC-MRI data. These are issues related to the estimation of the AIF, namely delay and dispersion, that will be discussed in detail after.

The CBV is expressed as the relative amount of tracer in the voxel with respect to the arterial reference and is calculated as the ratio between the integrals of the tissue and arterial concentrations, according to

CBV = k H ρ t 0 C ts (t)dt t 0 C a (t)dt (9.11)
which is normally expressed in milliliters per 100 grams of tissue (ml/100g), or in % by omitting the density and the hematocrit constant, as it will be the case in the rest of this thesis. An early example of the contrast in the brain is reported in fig. 9.2. The mean transit time (MTT) describes the average time it takes for a particle of tracer entering the vascular system, underlying a voxel, to leave it definitely. According to the central volume theorem expressed in eq. ( 9.5) it can be expressed as the ratio of the other two parameters as MTT = CBV CBF .

(9.12)

However, despite the general validity of eq. (9.12), [START_REF] Robert M Weisskoff | Pitfalls in mr measurement of tissue blood flow with intravascular tracers: which mean transit time?[END_REF] show that in DSC-MRI practice, this equation does not correspond to eq. ( 9.3) but rather to

MTT = ∞ 0 t • h(t)dt ∞ 0 h(t)dt (9.13)
which is the correct formulation of the dependence between the mean transit time and the transport function h(t) [START_REF] Calamante | Measuring cerebral blood flow using magnetic resonance imaging techniques[END_REF].

The use of the central volume theorem offers another way to compute the cerebral blood volume [Vonken et al., 1999b, Bjørnerud and[START_REF] Bjørnerud | A fully automated method for quantitative cerebral hemodynamic analysis using DSC-MRI[END_REF] CBV = R(t)dt (9.14) which is dependent, this time, on the recovered response function, and it is therefore sensitive to the adopted deconvolution strategy. Hence, it is useful to employ eq. (9.14) rather than eq. (9.11) when the goal is performing comparisons to rank different techniques according to the error in recovering CBV.

It appears straightforward to stress the importance of performing a good estimation of the response function, R(t), since the estimation of the perfusion parameters, CBF, CBV and MTT, and of the residue function, r(t), is heavily dependent on that. A whole literature about what is the right shape of R(t), or better the shape of the residue function r(t), has flourished during the past. Although the actual shape reconstructed in vivo is affected by phenomena such as dispersion, it is worth summarizing the main models used to represent the residue function.

The residue function

The knowledge of the shape of the response function R(t) is fundamental to obtain a good estimation of the cerebral blood flow, since CBF = R(0) as previously shown.

In the absence of artifacts, such as dispersion, the shape of R(t) corresponds to that of the residue function r(t). Knowing the actual shape could potentially lead to the characterization of the transport function h(t) and help inferring properties of the tissue perfusion. However, the actual shape is unknown and speculations about it can be made via experimental assessments, simulations or pharmacokinetic assumptions.

The most popular formulation of the residue function is the mono-exponential one, which considers the capillary bed as one single well-mixed compartment [START_REF] Jacquez | Compartmental analysis in biology and medicine[END_REF][START_REF] Lassen | Tracer kinetic methods in medical physiology[END_REF][START_REF] Bassingthwaighte | The cardiovascular system[END_REF] r exp (t) = e -t MTT (9.15) which was initially suggested by [START_REF] Seymour | Measurement of regional circulation by the local clearance of radioactive sodium[END_REF] and studied, to perform deconvolution of eq. (9.9), by Østergaard et al. [1996]. [START_REF] Lassen | The blood flow of the cerebral cortex determined by radioactive krypton85[END_REF] extended this model by assuming that the brain was made of a number of homogeneous (or first order) compartments arranged in parallel, hypothesis supported by [START_REF] Jacquez | Compartmental analysis in biology and medicine[END_REF], and proposed a bi-exponential model that accounts for fast and slow flowing capillary components

r bi-exp (t) = f • e -τ F t + (1 -f ) • e -τ S t (9.16)
where τ F and τ S are the fast and slow time-rates respectively, and f specifies the relative weight of the fast component. Zierler [1965] advocated against the use of the mono-exponential model, providing evidence that it does not entirely describe experimental results. Other expressions motivated by the empirical observation were proposed such as the Fermi [START_REF] Axel | Tissue mean transit time from dynamic computed tomography by a simple deconvolution technique[END_REF] and the Lorentzian [START_REF] Calamante | Quantification of bolus-tracking MRI: Improved characterization of the tissue residue function using Tikhonov regularization[END_REF]] models respectively

r lor (t) = 1 1 + πt 2 MTT 2 (9.17) r f er (t) = 1 + e -µ κ 1 + e t-µ κ (9.18)
where in eq. ( 9.18) µ and κ are shape parameters. This relation aims to model a situation in which the mixing of the indicator with blood is slower than that modeled by the exponential decay. One other expression proposed by [START_REF] Mouridsen | Bayesian estimation of cerebral perfusion using a physiological model of microvasculature[END_REF], called the Vascular Model (VM), models the capillary bed transit times with a continuous gamma distribution leading to

r vm (t) = ∞ t θ α-1 e -θ/β β α Γ(α) dθ α, β > 0 (9.19)
where α and β are shape parameters and Γ is a gamma function.

In a recent work, [START_REF] Mehndiratta | Modeling the residue function in DSC-MRI simulations: Analytical approximation to in vivo data: Modeling Residue Function for DSC-MRI[END_REF] performed a comparison of the models on data from normal subjects and atherosclerotic patients. They show that the shape of r(t) is best described in vivo by a bi-exponential model, a finding that is in agreement with previous literature, such as the work of Park and Payne [2013].

Regardless of how the residue function is represented, to recover it there exist some physiological and technical issues that have to be accounted for. Just to mention, the simple assumption that r(0) = 1 and R(0) = CBF is no longer valid in the presence of dispersion. Indeed, if dispersion is not taken into account, the estimated r(t) is no longer monotonically decreasing. The next section describes the technical details of perfusion processing and the related non-ideal behaviors such as delay and dispersion.

Practical issues and non-ideal behaviors

In order to formulate the input-output problem of the convolution integral expressed in eq. ( 9.9), it is important to define what are these input and output. In the processing, each tissue-voxel is considered to have a concentration C ts (t) which is the result of the convolution of the response function R(t) and the arterial input, C a (t). However, with DSC-MRI it can only be measured a concentration time-curve per voxel, therefore the actual arterial input is unknown. Physiological knowledge of the brain suggests that the tracer injected into the vascular system passes in some well-defined arteries before reaching the brain tissue. For this, reason the most common procedure to define C a (t) consist on the identification of these arterial areas -such as the middle cerebral artery (MCA) [START_REF] Porkka | Arterial input function measurement with mri[END_REF], Rosen et al., 1991] shown in fig. 9.3 -and the selection of the signal therein. This technical choice3 for the arterial input is the source of non-ideal behaviors with respect to the theory described in section 9.3, which are related to a physiologic non-representation of the actual arterial input by the selected global C a (t). These issues, addressed in this thesis, are the aforementioned time-delay and the bolus dispersion. For completeness, here is reported also the complication related to the physiological phenomenon of the recirculation of the tracer.

Time-delay

Once the bolus of tracer -injected in the vascular system -reaches the brain through an artery such as the MCA, from there it starts being transported within the tissue via the microcirculation. Clearly, the bolus arrival time measured in a voxel will be delayed compared to that measured at the location of the artery that "feeds" it, and the delay will be, in principle, longer as the considered voxel is further away from this location. This can be visualized in the right image of fig. 9.3, where the decay of the tissue signal -measured at the location indicated by the green dot in the left map -seems to be posterior to that of the selected AIF, that corresponds to the yellow dot. However, the actual bolus arrival time is difficult to predict due to the complex vascular organization of the brain tissue. In general, the measured tissue concentration, C ts (t), is said to be delayed compared to the selected AIF, C a (t), of an unknown time-delay τ [START_REF] Calamante | Delay and Dispersion Effects in Dynamic Susceptibility Contrast MRI Simulations Using Singular Value Decomposition[END_REF]]. In the case of wrong choice of the AIF the delay can also assume negative values. The presence of a non-zero delay affects the validity of the convolution expressed eq. ( 9.9) due to incorrect causality assumptions, and this must be handled.

Bolus dispersion

The second main issue related to the choice of the AIF is the dispersion of the bolus of tracer injected into the vascular system. Indeed, due to the nature of the vascular network the bolus may naturally spread along the way from the AIF location to the voxel of interest, where the tissue concentration C ts (t) is measured [START_REF] Østergaard | Cerebral blood flow measurements by magnetic resonance imaging bolus tracking: comparison with [15o] h2o positron emission tomography in humans[END_REF][Østergaard et al., , Calamante et al., 2000]]. Conversely, dispersion can appear in the case of pathology such as a steno-occlusive disease, as studied via simulations by [START_REF] Calamante | Estimation of bolus dispersion effects in perfusion MRI using image-based computational fluid dynamics[END_REF]. This is likely the case reported in fig. 9.3, where the dispersed signal corresponding to the green dot, sampled in the hypoperfused infarcted region of interest (ROI), shows a spread profile. It is extremely complicated to appreciate dispersion just by looking at the signal or at the corresponding concentration time-curve. As it will be clearer after, this thesis contributes to it with the introduction of a shape descriptor, namely the dispersion time, to characterize dispersion.

In DSC-MRI, data can generally contain voxels with no dispersion, dispersed voxels within healthy tissue, and voxels where dispersion is the result of a pathological scenario. Dispersion is related to macrovascular phenomena and constitutes a confounding factor for the estimation of the microvascular response function, R(t), and tissue perfusion parameters. This thesis also contributes with a framework to perform perfusion deconvolution that, additionally to the estimation of the delay, is compliant with dispersion.

Recirculation

The experimental acquisition time of the signal S(t) is sufficiently long to include an undesired part due to the recirculation of the tracer within the voxel-region. As visible in the signals of fig. 9.3, after the decay corresponding to the first passage of the tracer, the signal intensity does not return to the baseline value, i.e. the steady intensity before the signal decay. The perfusion theory, therefore, holds only for the first passage of the tracer [START_REF] Meier | On the theory of the indicator-dilution method for measurement of blood flow and volume[END_REF]. A solution for this consists in correcting for this bias, typically with a gamma-variate function [START_REF] Starmer | Computer computations of cardiac output using the gamma function[END_REF]Clark, 1970, Calamante et al., 1999].

Challenges: Delay and Dispersion

This thesis tackles the technical issues related to recovering the response function R(t) and the perfusion parameters in the presence of delay and dispersion, by proposing a delay/dispersion-compliant framework to perform the deconvolution of eq. ( 9.9), and by integrating it with existing techniques that characterize dispersion to ameliorate their performance. In fact, on one side the delay complicates the deconvolution, and on the other one dispersion leads to the misestimation of the perfusion parameters. However, for a long time the macrovascular effect of dispersion, although acknowledged, was often neglected for the sake of finding the actual microvascular perfusion described by an ideal R(t).

The next chapter introduces dispersion mathematically, and gives a wide view over its effects on the response function. This will lead to the definition of effective, i.e. potentially dispersed, response function and perfusion parameters. After an accurate phenomenological observation of these effects, the chapter will introduce a method to characterizing dispersion effects with the dispersion time, τ δ .

More specifically, the contributions of this thesis are structured in the following chapters:

chapter 10 Presents the mathematical description of dispersion, and the proposed phenomenological interpretation of it; it also reports the effects of dispersion on the response function and perfusion parameters, and introduces the proposed dispersion time to characterize them; the chapter will clarify the reason why, in the presence of dispersion, the solution of the deconvolution eq. ( 9.9) leads to an effective response function and perfusion parameters;

chapter 11 proposes a description of the effective response function by means of Dispersion-Compliant Bases (DCB): these are proposed in a non-linear and linearized versions; the chapter also reports details about state-of-the-art techniques, which have been implemented in this thesis to perform comparisons; after proposing DCB deconvolution, it is also proposed a method to recover dispersion-free perfusion parameters, i.e. CBF and MTT;

chapter 12 presents experiments to validate in silico and in vivo the performance of the proposed DCB deconvolution framework; it shows that DCB has higher precision in estimating the effective response function and perfusion parameters compared to the state-of-the-art; it will also be shown that the use of the framework as a pre-processing step improves the accuracy of dispersionfree estimates of the perfusion parameters when adopting a model to describe dispersion;

chapter 13 proposes a method developed to infer in vivo the correct model to adopt in order to obtain perfusion parameters that are free from the bias due to dispersion, i.e. CBF and MTT; the chapter also provides evidence of the importance of performing perfusion deconvolution with the least amount of assumptions about the possible undergoing dispersion process, which is the case of the proposed DCB approach.

A summary of the contributions is given in chapter 14.

Overview

What is dispersion, how is it formalized, why is it fundamental to take it into account? How does the response function change, in the presence of dispersion, into an effective one? How do the perfusion parameters CBF and MTT, are redefined as effective perfusion parameters, CBF * and MTT * ? This chapter answers these questions and, in addition, and proposes ways of characterizing dispersion in order to exploit it rather than just considering it as an artifact, which is often done in the literature. This chapter further sets the motivations behind the development of a "dispersion-compliant" framework for perfusion deconvolution, which will be introduced in the next chapter.

Keywords effective response function; effective perfusion parameters; vascular transport function dispersion time.

Confounding factor or useful insight?

The vision forwarded in this thesis is that dispersion can be a good indicator of pathological conditions related to steno-occlusive diseases. Indeed, dispersion is often regarded as a technical issue related to the non-representation of the actual arterial input by the selected arterial input function (AIF). The presence of dispersion constitutes a bias for the estimation of perfusion parameters. For this reason, it has been previously tackled at a fundamental level with acquisition methods [START_REF] Alsop | Defining a local input function for perfusion quantification with bolus contrast mri[END_REF][START_REF] Calamante | Defining a local arterial input function for perfusion mri using independent component analysis[END_REF][START_REF] Sourbron | A tracer-kinetic field theory for medical imaging[END_REF]] that try to avoid the possibility of dispersion 1 . Instead, the non-representation issue could be exploited, for instance, to describe the mismatch between the response function of healthy and pathological tissues: an AIF selected in a healthy tissue is not representative of the actual arterial input of a pathological tissue, thus the recovered response function in the latter will be dispersed and far from ideal. Still, the majority of works in the literature neglect the bias due to the presence of dispersion. However, some methods aim at the minimization of dispersion effects from a computational point of view [START_REF] Willats | Improved deconvolution of perfusion MRI data in the presence of bolus delay and dispersion[END_REF][START_REF] Willats | Minimising the effects of bolus dispersion in bolus-tracking MRI[END_REF]. Some others aim at the restoration of the dispersion-free perfusion parameters [START_REF] Zanderigo | Nonlinear Stochastic Regularization to Characterize Tissue Residue Function in Bolus-Tracking MRI: Assessment and Comparison With SVD, Block-Circulant SVD, and Tikhonov[END_REF][START_REF] Mouannes-Srour | Correction for arterial-tissue delay and dispersion in absolute quan-titative cerebral perfusion DSC MR imaging[END_REF], whereas the work of Mehndiratta et al. [2014a] suggests a model-based quantification of dispersion.

The characterization of dispersion, as it will be clearer in this chapter, starts with a good estimation of the response function R(t) which, due to the possible presence of dispersion is denoted with a " * ", i.e. R * (t), to indicate an effective response function. However, obtaining a good estimation revealed to be an herculean challenge due to the low signal-to-noise (SNR) of available DSC-MRI data combined with the need of producing highly non-linear models for dispersion quantification. Moreover, as this thesis contributes to clarify, in a particular voxel-region dispersion can be absent or present. In the latter case, dispersion can have different "intensities", and can have different nature, i.e. different possible models that explain the phenomenon. Existing methods employ heuristics or assume models which can condition the estimation of R * (t) via deconvolution, thus biasing results. For this reason, an important purpose of this thesis is the development of a robust framework to estimate and characterize R * (t) in the presence or absence of dispersion, and without any formal modeling of dispersion itself.

This chapter presents the mathematical definition of dispersion. Then, based on phenomenological observation, it proposes a method to characterize it with a single parameter: the dispersion time. This is a basic but fundamental step towards the characterization of dispersion without the use of models or heuristics, and it will be of help for dispersion quantification as it will be clearer later on in this chapter.

Dispersion formalism

Dispersion has been widely studied in literature [START_REF] Calamante | Delay and Dispersion Effects in Dynamic Susceptibility Contrast MRI Simulations Using Singular Value Decomposition[END_REF], 2002[START_REF] Calamante | Estimation of bolus dispersion effects in perfusion MRI using image-based computational fluid dynamics[END_REF][START_REF] Daniel | An introduction to computational diffusion mri: the diffusion tensor and beyond[END_REF][START_REF] Calamante | Bolus dispersion issues related to the quantification of perfusion MRI data[END_REF][START_REF] Ko | Reexamining the quantification of perfusion MRI data in the presence of bolus dispersion[END_REF][START_REF] Willats | The role of bolus delay and dispersion in predictor models for stroke[END_REF], Chappell et al., 2013]. It can be characterized by a Vascular Transport Function (VTF): this is the probability density function VTF(t) of the vascular transit times t. In fact, the bolus of tracer may undergo dispersion along the way to the voxel of interest, where the tissue concentration C ts (t) is measured. This effect is mathematically described as a convolution between the selected AIF, C a (t), and the Vascular Transport Function [START_REF] Calamante | Delay and Dispersion Effects in Dynamic Susceptibility Contrast MRI Simulations Using Singular Value Decomposition[END_REF]. Therefore, the effective, i.e. potentially dispersed, arterial input, C * a (t), to the voxel of interest is

C * a (t) = C a ⊗ VTF(t) (10.1)
that, when inserted in the convolution integral of eq. (9.9)2 , leads to

C ts (t) = CBF • C * a ⊗ r(t) (10.2)
which can be rewritten as where

C ts (t) = CBF • {C a ⊗ VTF(t)} ⊗ r(t) = CBF • C a ⊗ {VTF ⊗r(t)} (t) = C a ⊗ {VTF ⊗ [CBF • r] (t)} (t) = C a ⊗ {VTF ⊗R(t)} (t) = C a ⊗ R * (t) (10.3) ⊗ ⊗ C a (t)
R * (t) = R ⊗ VTF(t) (10.4)
is the effective dispersed tissue response function, which corresponds to the original only in the absence of dispersion, i.e. R * (t) = R(t) ⇐⇒ VTF(t) = δ(t).

As stated earlier, in the presence of dispersion the function VTF(t) is the probability density function of the occurrence of the vascular transit time t. The shape of VTF is an expression of the vascular dynamic and is unknown. In the literature there are different models for it, called dispersion kernels (DKs). These are the exponential, EDK [START_REF] Calamante | Delay and Dispersion Effects in Dynamic Susceptibility Contrast MRI Simulations Using Singular Value Decomposition[END_REF], the lognormal, LNDK [START_REF] Willats | Modelling the bolus dispersion from DSC-MRI data[END_REF][START_REF] Calamante | Arterial input function in perfusion MRI: a comprehensive review[END_REF], and the gamma, GDK [Chappell et al., 2013],

VTF(t) =            EDK(t, β) = βe -βt LNDK(t, µ, σ) = 1 tσ √ 2π e -(ln(t)-µ) 2 2σ 2 GDK(t, p, s) = s 1+sp Γ(1+sp) t sp e -st (10.5)
where β, µ, σ, p, s are shape parameters. Figure 10.1 clarifies the role of the VTF in the forward modeling of tissue perfusion. Particularly, dashed curves correspond to the non-dispersed case, and solid curves to the presence of dispersion with color Table 10.1: Dispersion kernels and related shape parameters for the case of low, medium and high dispersion as given by Mehndiratta et al. [2014a].

low medium high EDK β = 1 β = 1/2 β = 1/4 LNDK µ = -1, σ = 1 µ = -0.15, σ = 0.75 µ = 0.59, σ = 0.78 GDK p = 1, s = 2 p = 3, s = 1 p = 5, s = 0.5
codes indicating the influence of the different DKs in eq. ( 10.5). Together with the "nature" of dispersion, i.e. the DK model of VTF, dispersion levels can be encoded. Table 10.1 reports the shape parameters in case of low, medium, and high dispersion as given by Mehndiratta et al. [2014a]. This section clarifies that, in the presence of dispersion, data support the effective response function R * (t) instead of R(t). As it will be specified in the next section, one of the immediate consequences of this is the fact that the actual cerebral blood flow cannot be estimated as the peak of the effective response function, i.e. R * (0) = CBF if VTF(t) = δ(t). Instead, an effective one, CBF * , can be calculated as the maximum of R * (t). The next section explains the dispersion effects on the effective response function and perfusion parameters, and propose to summarize them with the dispersion time.

Dispersion effects

The content of this section is based on the analysis of dispersion effects published in Pizzolato et al. [2016b]. Similarly to what shown for the transport function, h(t), in section 9.2, for a specific vascular transport function (VTF) it can be specified a vascular mean transit time MTT v , from the site of the measurement of the AIF, C a (t), to the voxel of interest [START_REF] Calamante | Delay and Dispersion Effects in Dynamic Susceptibility Contrast MRI Simulations Using Singular Value Decomposition[END_REF].

In order to illustrate the effect of MTT v on the measured perfusion parameters, we can derive the effective response function R * (t) when VTF(t) = EDK(t, β). In fact, in this case it is MTT v = 1/β: a low value of β corresponds to a pronounced dispersion effect. Then, in the case of a bi-exponential residue function, in eq. ( 9.16), the effective response function is calculated according to eq. ( 10.4) as

R * (t) = -β (β -τ F )(β -τ S ) [(f (τ F -τ S ) + β -τ F )e -βt + f (τ S -β)e -τ F t + (f (β -τ F ) -β + τ F )e -τ S t ] • CBF (10.6)
which is valid for β > 0 and depends on the original mean transit time MTT = f /τ F + (1 -f )/τ S . The influence of dispersion on the effective response function is shown in fig. 10.2. Here curves are generated with eq. (9.16) -multiplied by CBF -in the absence of dispersion (purely decaying function in dashed line), and with eq. (10.6) in case of dispersion. The images show that with dispersion the response function is nonmonotonic, i.e. it has an increasing-decreasing profile, and that its maximum does not correspond to the beginning of the curve. Particularly, the top image shows two important dispersion effects: the decrease of the peak and the broadening of the shape. These effects seem to increase with MTT v , i.e. 1/β, in the direction of the arrow. Indeed, the effective response function tends to a purely decaying curve (a bi-exponential) for MTT v → 0, i.e. β → ∞. In this image MTT is constant. However, in the bottom image MTT v is constant while MTT varies. It shows that, for a specific value of MTT v > 0, the dispersion effects are enhanced when MTT is high. The figure also illustrates CBF, i.e. the peak of the dispersion-free curve in the top image, and the effective blood flow CBF * , i.e. the peak of a dispersed curve. In these circumstances, the effective cerebral blood flow, CBF * , becomes CBF * = max{R * (t)} (10.7) and, consequently, the effective mean transit time, MTT *

MTT * = CBV CBF * (10.8)
which is the equivalent of eq. ( 9.12) derived from the central volume theorem discussed in section 9.3.

What if dispersion is neglected?

Many techniques in the literature, the majority indeed, neglect the bias due to the presence of dispersion. Therefore they retrieve the effective perfusion parameters but they treat them as if they were dispersion-free. Alternatively, these techniques force the solution of the response function to comply with an ideal one, i.e. monotonically decaying, thus distorting the relationship between the sought solution and the nature of acquired data. It is to be understood that, as it will be better exemplified in the next chapter, recovering the response function is matter of deconvolution, which is ill-posed and actually allows for a camel to pass through the eye of a needle3 . As shown in fig. 10.2, the presence of dispersion causes a drop of the peak of the response function: if this is not taken into account then CBF is underestimated, and MTT is consequently overestimated [START_REF] Calamante | Delay and Dispersion Effects in Dynamic Susceptibility Contrast MRI Simulations Using Singular Value Decomposition[END_REF].

Dispersion time

A central phenomenological observation that can be done looking at fig. 10.2 is that the position of the peak of the response function, with respect to the beginning of the curve itself, is sensitive to the aforementioned dispersion effects. The beginning of the curve corresponds to the time delay, τ , existing between the selected AIF and the considered tissue concentration, as explained in section 9.5.1. However, due the commutativity of the convolution, it is possible to "transfer" the delay to the response function as follows

C a (t -τ ) ⊗ R * (t) = C a (t) ⊗ R * (t -τ ) (10.9)
where the delay determines the beginning of the function as represented in fig. 10.2. It can be noted that the time between τ and the peak of the response function, t max , is sensitive to the aforementioned dispersion effects, i.e. the joint contributions of dispersion and mean transit time. After defining the time of maximum as

t max = argmax t {R * (t)} (10.10)
Algorithm 1 Calculates dispersion-free parameters CBF and MTT (10.11) which corresponds to the time-to-peak of R * (t). A positive value of dispersion time (τ δ > 0) reveals the presence of dispersion whereas a value equal to zero (τ δ = 0) implies the absence of dispersion. We propose to detect and characterize dispersion effects by calculating the dispersion time τ δ directly on the effective response function. Since this parameter depends only on the recovered R * (t), its estimation is general and does not assume a specific model for the underlying VTF.

Require: R * (t), CBV, τ τ δ = argmax t {R * (t)} -τ if τ δ > 0 then R(t) = R * ⊗ -1 VTF(t) CBF = R(0) else CBF = R * (0) end if MTT = CBV / CBF the dispersion time, τ δ , is defined as τ δ = t max -τ

A serviceable aid for quantification

The idea behind dispersion time, despite its simple derivation, is fundamental in this thesis for the purpose of dispersion characterization. A large value of τ δ corresponds to pronounced dispersion effects, with all the physiological implications of it. However, in this thesis the dispersion time assumes also a procedural importance for dispersion detection and quantification. Indeed, dispersion can be quantified by using a model for VTF, such as those in eq. (10.5), but this is delicate since dispersion is not always present. As it will be clearer in the development of the thesis, the dispersion time offers a way to detect dispersion and to take the decision of quantifying it only where dispersion is actually present. The algorithm is described in algorithm 1. This reduces the risk of overfitting non-linear parameters of complicate models. However, the quantification of τ δ is challenging because its estimation is subordinated to that of the delay, τ , and time of maximum, t max . These, in turn, are extracted from the recovered response function, which is notoriously difficult to estimate especially in the proximity of the discontinuity, i.e. the beginning of the curve. For this reason, it is of paramount importance to recover the response function, R * (t), as robustly as possible. This further adds on the motivation of developing a new, robust, and dispersion-compliant deconvolution framework for perfusion MRI. This will be presented in the next chapter. The next section reports the proposal for another feature of R * (t) which is connected to dispersion and allows its characterization without the full knowledge of the VTF.

Dispersion index

Another feature of the response function, in the presence of dispersion, is the broadening of its shape. In fact, as dispersion effects increase, the shape of the curves in fig. 10.2 deviates from a purely decaying one. This deviation can be quantified as the difference between the integrals of the decreasing and increasing parts of the curve, normalized by the total area. This is the dispersion index, δ,

δ = ∞ τ +τ δ R * (t)dt - τ +τ δ τ R * (t)dt ∞ τ R * (t)dt (10.12)
where 0 < δ ≤ 1, under the assumption that the area of the decreasing part of curve is larger than that of the increasing one. In this case, absence of dispersion is indicated by δ = 1, whereas the presence is indicated by a lower value. It can be noted that the denominator of eq. ( 10.12) corresponds exactly to CBV according to eq. (9.14). As analyzed in Pizzolato et al. [2016b] 4 , the value of the index δ, despite the physiological implications that are beyond the scope of this thesis, can be interpreted in the light of τ δ and MTT * . Henceforth, for the purpose of building a perfusion deconvolution framework that is dispersion-compliant, only the dispersion time will be considered.

Conclusion

Dispersion is a complex phenomenon that cannot be neglected and that, instead, can be of help in characterizing pathological scenarios. The use of a parameter such as the dispersion time allows detecting and featuring dispersion effects, and deciding algorithmically whether it is the case to proceed with model-based quantification. The next chapter proposes to perform deconvolution with Dispersion-Compliant Bases (DCB) that approximate the response function, even in the presence of dispersion, without making assumptions about the presence, amount or nature of dispersion itself. It also describes two deconvolution techniques selected from literature, to be compared with DCB, that are "compliant" with the presence of dispersion: these are the oSVD [START_REF] Wu | Tracer arrival timing-insensitive technique for estimating flow in mr perfusion-weighted imaging using singular value decomposition with a block-circulant deconvolution matrix[END_REF] and CPI+VTF [Mehndiratta et al., 2014a] methods. The proposed deconvolution framework, object of the main contribution of this thesis, will be presented in two version, based on non-linear and linear formulations of DCB.

Overview

The bias due to the effects of dispersion is neglected by the majority of the techniques for DSC-MRI deconvolution. However, there are some that are "compliant" with it or, eventually, explicitly aim at its quantification. After presenting the motivations behind the development of a deconvolution framework, the chapter presents the gold-standard techniques in the field, and then introduces the proposed Dispersion-Compliant Bases (DCB). These are capable of representing the effective response function. At the end, it is also presented a method to recover dispersion-free CBF and MTT.

Keywords estimation; bases; models.

Motivations

Several techniques have been proposed to perform perfusion deconvolution. The solution for the response function can be sought by adopting the mono/bi-exponential models from pharmacokinetic [START_REF] Jacquez | Compartmental analysis in biology and medicine[END_REF][START_REF] Bassingthwaighte | The cardiovascular system[END_REF][START_REF] Lassen | The retention of [99mtc]-d, l-HM-PAO in the human brain after intracarotid bolus injection: a kinetic analysis[END_REF] or the vascular model [START_REF] Mouridsen | Bayesian estimation of cerebral perfusion using a physiological model of microvasculature[END_REF]. Other approaches adopt the truncated singular value decomposition tSVD [Østergaard et al., 1996] and its variants [START_REF] Wu | Tracer arrival timing-insensitive technique for estimating flow in mr perfusion-weighted imaging using singular value decomposition with a block-circulant deconvolution matrix[END_REF][START_REF] Smith | Removing the effect of svd algorithmic artifacts present in quantitative mr perfusion studies[END_REF][START_REF] Bjørnerud | A fully automated method for quantitative cerebral hemodynamic analysis using DSC-MRI[END_REF], Tikhonov regularization [START_REF] Calamante | Quantification of bolus-tracking MRI: Improved characterization of the tissue residue function using Tikhonov regularization[END_REF], non-linear stochastic regularization [START_REF] Zanderigo | Nonlinear Stochastic Regularization to Characterize Tissue Residue Function in Bolus-Tracking MRI: Assessment and Comparison With SVD, Block-Circulant SVD, and Tikhonov[END_REF], maximum entropy principle [START_REF] Charter | Maximum entropy and its application to the calculation of drug absorption rates[END_REF]Gull, 1987, Drabycz et al., 2005], maximum likelihood estimation [Vonken et al., 1999a], fully bayesian methods [START_REF] Boutelier | Bayesian hemodynamic parameter estimation by bolus tracking perfusion weighted imaging[END_REF], wavelets [Wirestam andStåhlberg, 2005, Connelly et al., 2006], exponential approximation [START_REF] Stephen | Deconvolution for DCE-MRI using an exponential approximation basis[END_REF], ARMA modeling [START_REF] Batchelor | ARMA regularization of cardiac perfusion modeling[END_REF], and control point interpolation [START_REF] Mehndiratta | A control point interpolation method for the nonparametric quantification of cerebral haemodynamics from dynamic susceptibility contrast MRI[END_REF]. Some of these approaches allow for time-delay estimation/insensitivity and/or for non-negativity of the estimated response function but do not account, at least explicitly, for dispersion.

However, there have been attempts to quantify dispersion effects by adopting a model of the vascular transport function (VTF), to be estimated directly during a non-linear deconvolution process. Unfortunately, the shape of the VTF is not known and different models support as different concentration time-curves, C ts (t), as demonstrated by fig. 10.1 in the previous chapter. The assumption of one single model might be reductive and can bias the results, which is likely in an ill-posed deconvolution procedure. This motivates the choice of developing a deconvolution framework that makes no assumptions about the dispersion kernel (DK), be it exponential, gamma, lognormal, or other. The choice falls into a signal-driven procedure based on Dispersion-Compliant Bases (DCB). The purpose is, for some aspects, similar to the well-known family of techniques that make use of the truncated singular value decomposition (SVD) deconvolution developed by Østergaard et al. [1996] and [START_REF] Wu | Tracer arrival timing-insensitive technique for estimating flow in mr perfusion-weighted imaging using singular value decomposition with a block-circulant deconvolution matrix[END_REF]. In fact, both SVD and the proposed DCB are able to recover the effective response function without assumption about VTF, but the solution is sought in very different ways and with much different results.

The SVD technique suffers from some well-studied limitations [START_REF] Bjørnerud | A fully automated method for quantitative cerebral hemodynamic analysis using DSC-MRI[END_REF], such as strong oscillations in the solution for R * (t), lack of nonnegativity assumptions, SNR-dependent optimal settings, underestimation of CBF, and consequent overestimation of MTT. Nevertheless, its ease of use, intelligibility, low computational cost, and acceptable results make of it the most widespread technique for perfusion deconvolution in DSC-MRI. Among the attempts of estimating the effective response function with the desired properties, a notable one is the technique called Control Point Interpolation with Vascular Transport Function (CPI+VTF) developed by Mehndiratta et al. [2014a], which shows comparably better performance than SVD. As the name suggests, the method performs deconvolution under the assumption that the VTF model is known, which is assumed to be a gamma dispersion kernel (GDK) as presented in eq. (10.5). Despite the method reveals to be more effective than SVD, it shares the issues related to the choice of a specific model to represent the VTF. Indeed CPI+VTF, as it will be clarified later, aims at estimating the dispersion-free CBF and MTT, which can be conditioned by a wrong assumption of the dispersion kernel. The following sections describes in details the adopted version of SVD, that is the block-circulant SVD with oscillation index (oSVD), and the aforementioned CPI+VTF.

Deconvolution with oSVD

The block-circulant SVD deconvolution with oscillation index (oSVD) developed by [START_REF] Wu | Tracer arrival timing-insensitive technique for estimating flow in mr perfusion-weighted imaging using singular value decomposition with a block-circulant deconvolution matrix[END_REF] is based on its ancestor truncated version, tSVD, proposed by Østergaard et al. [1996]. It is worth then starting with describing the latter. The tSVD approach represents the convolution of eq. (9.9) in a discretized form. In fact, by assuming the tissue concentration, C ts (t), and the AIF, C a (t), sampled on an equally spaced time grid of size M such as t 1 , t 2 , . . . , t M with ∆t = t i+1 -t i the convolution can be discretized as

C ts (t j ) = ∆t j i=0 C a (t i )R * (t j -t i ) (11.1)
where here R * is used to underline that this is an effective response function. Equation 11.1 can be written in matrix form as

       C ts (t 1 )
C ts (t 2 ) . . .

C ts (t M )        = ∆t        C a (t 1 ) 0 • • • 0 C a (t 2 ) C a (t 1 ) • • • 0 . . . . . . . . . . . . C a (t M ) C a (t M -1 ) • • • C a (t 1 )               R * (t 1 ) R * (t 2 ) . . . R * (t M )        (11.2) 158CHAPTER 11. DISPERSION-COMPLIANT BASES (DCB) DECONVOLUTION
and by incorporating ∆t in the matrix the problem can be formulated as

c t = Ar (11.3)
where c t is the M × 1 vector containing the tissue concentration samples, A is the M × M convolution matrix, and r is the vector containing the unknown samples of the response function. An obvious strategy for recovering the residue function would then be to adopt the raw deconvolution solution

r = A -1 c t (11.4)
which unfortunately is very sensitive to noise and can thus lead to unstable solutions. The tSVD approach tackles the problem of reducing the sensitivity to noise of convolution by adopting the singular value decomposition

A = USV T (11.5)
where S is a rectangular diagonal matrix with non-negative values known as singular values. Usually the singular values are sorted in descendant order along the diagonal where the smaller ones correspond to those giving oscillations or to high frequencies.

The method consists in selecting a threshold, p SVD , which is usually chosen in a range [5%, 20%] of the maximum singular value: the singular values below this threshold are truncated. After truncation the tSVD deconvolution is r = VWU T c t (11.6) where W corresponds to 1/S along the diagonal if the corresponding values of S are above p SVD , and zeros elsewhere. The chosen threshold highly influences the deconvolved solution and several techniques have been proposed. [START_REF] Liu | Cerebral blood flow measurement by dynamic contrast mri using singular value decomposition with an adaptive threshold[END_REF] proposed a method to perform a pixel-by-pixel thresholding based on the SNR evaluated at the maximum of the measured concentration time curve C ts (t). More robust methods such as the generalized cross validation (GCV) and the L-curve criterion [START_REF] Per | Analysis of discrete ill-posed problems by means of the l-curve[END_REF][START_REF] Per | The use of the l-curve in the regularization of discrete ill-posed problems[END_REF], are discussed by [START_REF] Sourbron | Deconvolution of dynamic contrast-enhanced mri data by linear inversion: Choice of the regularization parameter[END_REF]. [START_REF] Koh | Assessment of Perfusion by Dynamic Contrast-Enhanced Imaging Using a Deconvolution Approach Based on Regression and Singular Value Decomposition[END_REF] proposed to automatically select p SVD via a regression approach finding the transition points on the Picard plot. The literature on the subject is actually extensive. In order to tackle the time-delay between the AIF and the tissue concentration, as discussed in section 9.5. 1, Wu et al. [2003] proposed a block-circulant version of eq. (11.6), namely cSVD, which is insensitive to delay thanks to the use of a circulant convolution matrix A c . More precisely C ts (t) and C a (t) are extended by zero-padding up to a length L ≥ 2M to avoid aliasing. Then the new L × L squared circulant convolution matrix has the following entries

A c i,j = ∆t • C a (t i-j+1 ) for j ≤ i ∆t • C a (t L+i-j+1 ) for j > i . (11.7)
where i, j indicate the row and column index respectively.

This formulation of the convolution allows estimating a response function vector r c , of size L × 1, that is not compromised by the presence of delay, such that

r c = V c W c U T c c c t . (11.8)
However, when using cSVD the discontinuity in the correspondence of t = 0 causes leakage frequencies to be amplified, giving rise to spurious oscillations that dominate the deconvolved signal. Still Wu and colleagues proposed to limit oscillations with the oSVD technique, which modulates the threshold p SVD according to an oscillation index1 , OI

OI = 1 ∆t L CBF * L-1 k=2 |R * (k) -2R * (k -1) + R * (k -2)| (11.9)
such that p SVD is updated iteratively until the response function's OI falls below a user-defined value.

The oSVD technique is straightforward, and allows removing the influence of noise, thus stabilizing the deconvolution while obtaining a delay-insensitive estimation of the effective response function. However, it renders negative samples of the response function (which are non-physiological), does not allow estimating the delay, and always depends on a user-defined threshold which is an heuristic. In addition, as previously reported, it is known to underestimate the cerebral blood flow in specific scenarios. To this regard an intuitive analysis is here provided. For a given choice of p SVD , the singular values are truncated to remove the corresponding high frequencies. However, in the case of a dispersion-free scenario, as reported in the left image of fig. 11.1, the discontinuity at the time-delay (5 seconds) is a high frequency change, which is clearly cut by oSVD: see the ground-truth curve in dashed line and the reconstructed one in solid line. Indeed, the method implements a compromise between cutting oscillations -with a high frequency cutoff -and estimating the peak of R * (t): the two objectives would require, indeed, opposite solutions. On the other hand, when the ground-truth data is smoother, as in the dispersed case depicted in the right image of fig. 11.1, then the discontinuity at the time delay is less hard. As a consequence, the high frequency cutoff mainly cuts oscillations due to noise and "leaves untouched" the profile of the curve. Hence, in the presence of dispersion, which leads to smooth profiles of the response function, oSVD performs a good estimation. Note that the peak of the dispersed function is better recovered than that of the dispersion-free one. Indeed, the high frequency cutoff penalizes the recovery of a sudden change such as in the case of the dispersion-free response function.

CPI and CPI+VTF

The CPI+VTF deconvolution [Mehndiratta et al., 2014a] is a method based on Control Point Interpolation that overcomes the limitation of its ancestor CPI [START_REF] Mehndiratta | A control point interpolation method for the nonparametric quantification of cerebral haemodynamics from dynamic susceptibility contrast MRI[END_REF] by allowing the estimation of the effective response function in case of dispersion. The CPI technique uses consecutive control points opportunely constrained to represent the monotonic residue function r(t). Particularly, it uses η control points, some of which have "freedom of movement" both in time and amplitude. These points form virtual nodes to be interpolated using piece-wise cubic splines [START_REF] James M Hyman | Accurate monotonicity preserving cubic interpolation[END_REF]. A representation of this is shown in the left side of fig. 11.2

The first control point is fixed so that t 1 = 0 and r(t 1 ) = 1 to meet the definition r(0) = 1. The last point is set to the last point of the sampling grid t η = t M whereas its amplitude is let evolve, as for previous points. Said exceptions apart, the amplitude, r(t i ), and time position, t i , of the i-th control point are related to the previous as

r(t i ) = r(t i-1 ) • Ω i t i = t i-1 • ζ i .
(11.10)

The method implements the following empirical constraints for a total amount of 2η -1 parameters to be non-linearly estimated, including the delay. The model relies on an initial guess for the control points and CBF given by a mono-exponential deconvolution, based on the model reported in eq. ( 9.15). The parameters Ω i and ζ i are estimated on a logarithmic scale to render the model more linear with respect to them.

CBF ≥ 0 Ω i ≤ 1.5 ζ i ≥ ∆t/3 (11.11) gamma GDK t ⊗ t R * (t) = Control Point Interpolation R (it gives CBF,MTT)
The CPI method was previously compared to others [START_REF] Mehndiratta | A control point interpolation method for the nonparametric quantification of cerebral haemodynamics from dynamic susceptibility contrast MRI[END_REF]] revealing the best performance. However, it is designed for dispersion-free data, since it forces r(t 1 ) = 1, and it cannot handle the presence of dispersion. To overcome this issue Mehndiratta et al. [2014a] proposed a dispersion-compliant version of CPI, i.e. CPI+VTF, that simply convolves the CPI solutions of the residue function, multiplied by CBF, with the gamma dispersion kernel (GDK) of eq. (10.5), chosen as a model for VTF based on arguments given by Chappell et al. [2013]. A representation of the approach is shown in fig. 11.2. This method allows for quantification of dispersion by estimating the parameters p and s of the gamma model, while at the same time estimating the dispersion-free CBF, MTT, and the residue function r(t). However, the main limitation remains the assumption of a specific model for the VTF. In fact, chapter 13 of this thesis brings evidence that different brain regions support different models of VTF among those reported in eq. (10.5), i.e. EDK, LNDK. Secondly, here it is pointed out that in practice CPI+VTF hardly manages to represent a dispersion-free scenario because GDK(t, p, s) → δ(t) ⇐⇒ p = 0, s → ∞ (11.12) which is unfeasible. In fact, as pointed out by Mehndiratta et al. [2014a], CPI works better in dispersion-free scenarios.

Dispersion-compliant deconvolution

The motivation behind the development of a dispersion-compliant deconvolution procedure is to obtain a representation of the effective response function in both the absence and presence of dispersion while avoiding the explicit use of models, and overcoming the issues of the oSVD approach. The objective is to obtain smoother solutions for R * (t), without oscillations, while allowing for the estimation of the delay and the implementation of non-negativity constraints. The present section is based on work in Pizzolato et al. [2015aPizzolato et al. [ , 2016a]]. The dispersion-compliant deconvolution relies on efficient ways of representing the effective response function. In what follows, several dispersion-compliant representations are reported and compared. Among these, the main proposition consists on the Dispersion-Compliant Bases (DCB), which are firstly presented in their general non-linear version and then in their linearized version.

Non-linear Dispersion-Compliant Bases (DCB)

The proposed DCB approximate R * (t) with a sum of N bases where each basis is constituted as a sum of an exponential term and an exponential derivative term, with respect to the time rate, respectively weighted by two different constants a n and b n . In order to let the bases render the most general solution, two different time rates α n and β n , respectively for the exponential and derivative terms, are chosen. Figure 11.3 illustrates the basis. The sum of these bases can then represent any shape in the range between a multi-exponential decay and its derivative terms.

The proposed approximation has the following expression

R * DCB (t) = Θ(t) N n=1
(a n e -αnt + b n te -βnt ) (11.13) where Θ(t) is the Heaviside step function with no half-maximum convention, i.e. Θ(t) = 1 ∀t ≥ 0, and where the parameters a n , b n , α n and β n are unknowns to be estimated. Note that the time-rates α n , β n enter non-linearly in the basis formulation. The time-delay τ between C a (t) and C ts (t) can be explicitly considered in the response function approximation because the convolution is commutative as shown in eq. ( 10.9). Therefore, when taking into account the time-delay eq. ( 11.13) becomes

R * DCB (t) = Θ(t -τ ) N n=1
[a n e -αn(t-τ ) + b n (t -τ )e -βn(t-τ ) ] (11.14) where τ is an additional non-linear parameter to be estimated.

Deconvolution

As previously seen, the convolution problem in eq. (9.9) can be discretized by assuming that the tissue, C ts (t), and the arterial, C a (t), concentrations are sampled on an equally spaced time grid t 1 , t 2 , . . . , t M of size M , with ∆t = t i+1 -t i , as shown in eq. (11.1) The convolution is then formulated according to eq. ( 11.3), i.e. c ts = Ar, where A is the M × M convolution matrix containing the samples of the arterial input concentration, c ts contains the M samples of the tissue concentration and r contains the M unknown samples of the effective response function.

If a representation G for the response function is defined, where G maps the parameters to the response function according to eq. ( 11.13), then the deconvolution problem aims to find the related set of unknown parameters. The representation in eq. (11.13) 11.16) where the matrix has dimension M × 2N , and the vector 2N × 1. Therefore, the deconvolution problem incorporating eq. ( 11.13) can be solved as pL , pNL = argmin

G(p N L ) =        e -α 1 t 1 t 1 e -β 1 t 1 e -α 2 t 1 • • • t 1 e -β N t 1 e -α 1 t 2 t 2 e -β 1 t 2 e -α 2 t 2 • • • t 2 e -β N t 2 . . . . . . . . . • • • . . . e -α 1 t M t M e -β 1 t M e -α 2 t M • • • t M e -β N t M        p L = (a 1 , b 1 , a 2 , . . . , b N ) T (
p L ,p N L ||c ts -AG(p N L )p L || 2 2
(11.17) using a gradient-descent method (see the implementation section below).

The deconvolution problem incorporating eq. (11.14) includes the time-delay as additional non-linear parameter to be estimated. The time-delay, τ , is taken into account by considering the circular formulation of the convolution as specified in eq. ( 11.7)2 . Also in this case, the measured arterial and tissue concentration timecurves of length M are extended by zero-padding up to a length L ≥ 2M to avoid aliasing. The convolution problem of eq. (11.15) is reformulated as 11.18) where G τ (p N L ) is the L × 2N design matrix, that depends on τ and p N L , extended to the circular time sampling grid, and where the superscript " c " is dropped on vectors for clarity. The delay, τ , is estimated via grid search over a range [τ min , τ max ] seconds (chosen by the user), with a time step τ s ≤ ∆t. The estimated vectors of parameters pL , pNL can again be obtained via gradient-descent methods when the estimated delay τ is that minimizing

c ts = A c G τ (p N L )p L (
||c ts -A c G τ (p N L )p L || 2 2 among all τ ∈ [τ min , τ max ].
The use of bases like DCB offers the advantage of computing the analytic Jacobian matrix that can be supplied to the non-linear routine speeding up the estimation of the parameters. More details about the implementation of DCB deconvolution are now provided.

Implementation

The maximum basis order N is chosen after performing some synthetic experiments such that it renders a good compromise between explaining the variety in the data, i.e. the different dispersion kernels (exponential, gamma, and lognormal), and robustness to noise. The idea is to find the minimum order that allows minimizing the estimation error of the peak of the effective residue function, i.e. CBF * , in the presence of noise. To do so, effective response functions, R * (t), are synthetically generated by convolving the bi-exponential r(t) of eq. ( 9.16) with each of the various dispersion kernels: gamma (GDK), exponential (EDK), and lognormal (LNDK). The bi-exponential r(t) is generated with the median normal tissue parameters given by [START_REF] Mehndiratta | Modeling the residue function in DSC-MRI simulations: Analytical approximation to in vivo data: Modeling Residue Function for DSC-MRI[END_REF]: f = 0.97, τ F = 0.68 and τ S = 0.05. The dispersion kernel parameters are set to values corresponding to low, medium, and high dispersion, as reported in table 10.1. Secondly, for each generated R * (t)resulting from the combination of dispersion kernel and intensity -100 DCB fittings are performed for each order N ∈ {2, 3, 4, 5, 6, 7, 8}, corresponding to as many noisy realizations in the case of high noise regime, with SN R = 10. The order N = 5 is the minimum order reducing the CBF * estimation error below threshold.

The initial guess for DCB parameters is based on oSVD results to approximate a mono-exponential decay, particularly p L = (CBF oSV D , 0, . . . , 0) T and α n ,

β n = 1/M T T oSV D ∀n ∈ [1, N ].
DCB deconvolution is solved considering the minimization of the l 2 reconstruction norm with respect to the tissue concentration C ts (t). For this thesis the technique is implemented using in-house software written in MATLAB 3 , and the nonlinear optimization is performed using the trust-region-reflective algorithm within the lsqnonlin function.

Linearization motivations

The non-linear formulation of the Dispersion-Compliant Bases (NL-DCB) in eqs. (11.13) and (11.14) is general, and will be validated in the next chapter following the results in Pizzolato et al. [2016a]. However, NL-DCB shares with CPI+VTF the complication of having non-linear parameters: this maks the deconvolution slow and introduces the risk that the optimization procedure falls in a local minimum. In addition, NL-DCB shares with oSVD the lack of non-negativity constraints. The motivation behind a linearization of the deconvolution problem is to make it computationally cheap while preserving all of the desired qualities for the recovered response function:

• general formulation (capability of approximating any response function);

3 MATLAB Release 2015a, The MathWorks, Inc., Natick, Massachusetts, United States. • physiologically relevant and smooth recovered response function;

• delay estimation;

• non-negativity constraints.

The following sections 11.4.3 to 11.4.5 introduce three linearization strategies that will be compared in section 11.4.6.

Dispersed Exponential Approximation (DEA)

The response function in the absence of dispersion is substantially well described by a multi-exponential decay as argued in section 9.4. For this reason [START_REF] Stephen | Deconvolution for DCE-MRI using an exponential approximation basis[END_REF] proposed to represent the response function via exponential approximation (EA)

R(t) = Θ(t) N n=1
k n e -λnt (11.19) with k n linear parameters to be estimated and λ n time-rates to be fixed a priori.

Indeed [START_REF] Schwartz | Étude des sommes d'exponentielles réelles[END_REF] demonstrated that under certain conditions the approximation is non-parametric, according to the following theorem.

Theorem (Müntz). Given the time-rates {λ n } ∞ n=1 positive and satisfying lim n→∞ λ n = ∞, then the span of {e -λnt } is dense in L p [0, ∞) for 1 ≤ p < ∞ if and only if the series of the time-rate reciprocals diverges, i.e.

∞ n=1 1/λ n = ∞.
Following the theorem Keeling and colleagues proposed to fix the time-rates according to the harmonic law λ n = n/T, ∀n = 1, . . . , N (11.20) where T = t M is the length of the time observation interval of the acquired C ts (t). However, they did not take into account the presence of dispersion, therefore they aimed at obtaining a response function such that

• R(t) ≥ 0, t ≥ 0;

• dR(t)/dt ≤ 0, t ≥ 0; which leads to the impossibility of fitting an increasing-decreasing function such as R * (t) in the presence of dispersion.

In this thesis, the EA method is extended to account for dispersion, namely dispersed exponential approximation (DEA) by dropping the non-increasing condition on the derivative of the response function. Secondly, as pointed out by [START_REF] Batchelor | ARMA regularization of cardiac perfusion modeling[END_REF], the choice of the time-rates according to eq. ( 11.20) seems to be not flexible and tends to require high basis order N . Here it is proposed to scale the time-rates differently, and precisely according to a previous estimation of the mean transit time based on oSVD, such that the time-rates are more appropriate for the time range of the response function, which is expressed as

R * (t) = Θ(t) N n=1
k n e -nt P • MTT oSVD (11.21) where P is a parameter used to scale the denominator. The choice of this parameter will influence the approximation performance per number of basis functions. As an indication, here P = 4.

Deconvolution

The deconvolution accounting for eq. ( 11.21) requires, this time, a linear estimation.

The forward convolution problem can be expressed as 11.22) where G is the dispersed exponential approximation design matrix, and p L is the linear coefficient vector of parameters such that 11.23) where the matrix has dimension M × N , and the vector N × 1. The sought solution can be found with linear least-squares. However positivity can be enforced via quadratic programming such that pL = argmin 11.24) which can be solved with standard routines. In the next section, a generalization of the original EA technique, which however can recover with dispersed response functions, is reported.

c ts = AGp L (
G =         e - t 1
P • MTT oSVD e - 2t 1
P • MTT oSVD e - 3t 1
P • MTT oSVD
• MTT oSVD • • • e - N t M P • MTT oSVD         p L = (k 1 , k 2 , . . . , k N ) T (
Gp L ≥0 ||c ts -AGp L || 2 2 (
168CHAPTER 11. DISPERSION-COMPLIANT BASES (DCB) DECONVOLUTION

ARMA

After the release of the EA deconvolution, [START_REF] Batchelor | ARMA regularization of cardiac perfusion modeling[END_REF] proposed a generalization based on auto-regressive-moving-average (ARMA) modeling of the convolution. This corresponds to representing the discretized convolution as

C ts (m) = K k=0 b k C a (m -k) + J j=0
a j C ts (m -j) (11.25) which leads to an impulse response function of the associated system in the span of {t i-1 e -λ j t } q i=1 where q is the multiplicity of the corresponding root p j = e -λ j of the characteristic polynomial (11.26) After finding the roots p j , the time-rates λ j = -ln p j are known and eq. (11.19) can be applied to find the linear coefficients. However, here the root-finding step is skipped and the solution is obtained as the impulse response of the system corresponding to the estimated ARMA coefficients.

P (z) = 1 - J j=1 a j z -j .

Linearized Dispersion-Compliant Bases (L-DCB)

The idea is to linearize the general solution of NL-DCB, reported in eq. (11.13), to obtain a L-DCB. Note, that the basis in that equation includes an exponential derivative term which is corresponding to a pole with multiplicity 2 in the ARMA formulation. The first step to linearize eq. (11.13) is to reduce the number of exponentials such that 11.27) where each basis is still a sum of an exponential decay, for b n = 0, and its derivative term with respect to the time-rate, for a n = 0. However, the time rates α n still enter non-linearly in the formulation. Similarly to the derivation of DEA in section 11.4.3, the time-rates are chosen harmonically and according to a previous estimation of the mean transit time based on oSVD (11.28) and the delay can be taken into account as seen for eq. (11.14), precisely (11.29) which constitutes the representation of the effective response function via linear dispersion-compliant bases.

R * (t) = Θ(t) N n=1 (a n + b n t)e -αnt (
R * (t) = Θ(t) N n=1 (a n + b n t)e - nt P • MTT oSVD
R * (t) = Θ(t -τ ) N n=1 [a n + b n (t -τ )]e - n(t-τ ) P • MTT oSVD

Deconvolution

The representation of R * (t) reported in eq. ( 11.28) leads to a design matrix and vector of parameters as 11.30) of size M × 2N and 2N × 1 respectively. In the case where eq. ( 11.29) is considered, the estimation of the vector of linear coefficients can be performed considering the convolution as c ts = A c G τ p L (11.31) where G τ is the L × 2N design matrix depending on τ extended to the circular time sampling grid, and where c ts is the zero-padded tissue concentration vector. As for the non-linear version of DCB, the delay (τ ) is estimated via grid search over a range [τ min , τ max ] seconds (chosen by the user), with a time step τ s ≤ ∆t. The estimated vector of parameters pL is obtained as pL = argmin (11.32) where the estimated τ is such that it minimizes

G =         e - t 1
P • MTT oSVD t 1 e - t 1
P • MTT oSVD e - 2t 1
P • MTT oSVD • • • t 1 e - N
• MTT oSVD • • • t M e - N t M P • MTT oSVD         p L = (a 1 , b 1 , . . . , b N ) T (
G τ p L ≥0 ||c ts -A c G τ p L || 2 2
||c ts -A c G τ (p N L )p L || 2 2 among all τ ∈ [τ min , τ max ].

Comparison between linear techniques

In order to compare the three linear techniques, the root mean squared error (RMSE) of the effective response function, R * (t), fitting is computed. Response functions were generated in the presence of dispersion with gamma, exponential and lognormal kernels according to eq. ( 10.5) with dispersion intensities specified by values in table 10.1. 6 similarly reports a direct comparison between L-DCB and DEA, this time also including the estimation of CBF * and MTT * (second and third columns). Results indicate a higher stability for L-DCB, especially for SNR ≤ 30, for all of the comparison criteria. These results might be explained in the light of a more suitable shape of the basis in the case of L-DCB, which is closer to the shape of a dispersed response function as shown in fig. 10.2. Moreover, the derivative term te -αt arises from the natural representation of the convolution by means of an ARMA model, and in this sense L-DCB generalizes the dispersed exponential approximation while, at the same time, guaranteeing robustness to noise.

Finally, in the noise-free RMSE comparison it appears clear that L-DCB fitting error reaches zero for 12 parameters, corresponding to maximum basis order N = 6. This will be the adopted order when using the technique.

Dispersion-free perfusion parameters from DCB

Once estimated the effective residue function, R * (t), with either of the DCB techniques, only the effective perfusion parameters, CBF * and MTT * , can be estimated.

These correspond to the dispersion-free parameters, CBF and MTT, in the absence of dispersion. In order to detect the presence of dispersion it is proposed to calculate the dispersion time, τ δ , according to eq. (10.11). As stated in algorithm 1, if τ δ > 0 then a VTF model is used to deconvolve R * (t) to obtain R(t). To this purpose, although many other options are possible, in this thesis the inherent CPI+VTF model is directly fitted to the R * (t) obtained via the proposed DCB deconvolution. This two steps procedure -DCB deconvolution and CPI+VTF model fitting -is illustrated in fig. 11.4.

Conclusion

This chapter presents the perfusion deconvolution methods implemented for comparisons, such as oSVD, CPI+VTF, and ARMA, and introduces the proposed DEA, and DCB techniques (both linear and non-linear). The presented techniques allow estimating the effective response function, even in the presence of dispersion, without assumptions about the underlying dispersion. As an exception, CPI+VTF adopts a gamma dispersion kernel (GDK) as model for the vascular transport function. However, this last technique is used whenever a dispersion-free response function, R(t), needs to be recovered, as explained in section 11.5.

The next chapter 12 reports an evaluation of the dispersion-compliant bases, in their non-linear formulation, which reveal higher performance compared to the gold-standard oSVD, CPI, and CPI+VTF techniques in both dispersion-free and dispersed scenarios. The evaluation will take into account both in silico and in vivo experiments. The evaluation will also bring evidence that dispersion-free CBF and MTT are better recovered with a previous DCB deconvolution followed by a CPI+VTF model fitting (two steps procedure), rather than with a CPI+VTF deconvolution directly.

Chapter 13, reports an in vivo evaluation of the dispersion kernels. In fact, by using a combination of linear DCB and newly implemented variants of the CPI+VTF technique, evidence will be brought to show the limitations of assuming a specific VTF model in the brain tissue. LAB1 . The estimation of the parameters for CPI+VTF and DCB deconvolutions was solved considering the minimization of the l 2 reconstruction norm with respect to the tissue concentration C ts (t), as seen for instance in eq. (11.17). The non-linear optimization was performed using the trust-region-reflective algorithm within the lsqnonlin function, which improves performance over Levenberg-Marquardt in case the initial guess is far from the solution.

In silico data

This section describes the simulation of the MR arterial and tissue signals and concentrations used for the experiments in silico. The signals are obtained from the simulated concentration curves when CBF ,CBV, MTT, time-delay τ , dispersion kernel (GDK,EDK,LNDK), and dispersion level (none, low, medium, high as specified in table 10.1) are set. An illustration of the relationship between the generated arterial and tissue concentrations is provided in fig. 10.1.

The arterial concentration is calculated as the gamma-variate function [START_REF] Starmer | Computer computations of cardiac output using the gamma function[END_REF] C a (t) = 0 for t ≤ t 0 γ 0 (t -t 0 ) ν • e -(t-t 0 )/ξ for t > t 0 (12.1) with corresponding parameters -reported in table 12.1 -set according to the findings in [START_REF] Meijs | Analysis of perfusion mri in stroke: To deconvolve, or not to deconvolve[END_REF]. The tissue concentration C ts (t) is obtained by convolution between the arterial concentration C a (t) with the effective response function R * (t), according to eq. (11.1), here reported for clarity

C ts (t j ) = ∆t j i=0 C a (t i )R * (t j -t i )
where R * (t) = CBF • r biexp (t) in the absence of dispersion or R * (t) = CBF • r biexp ⊗ VTF(t) in case of dispersion where the VTF model is selected among the dispersion kernels shown in eq. ( 10.5). The parameters f, τ F , τ S of the residue function in eq. (9.16) were set according to the optimal median values given by [START_REF] Mehndiratta | Modeling the residue function in DSC-MRI simulations: Analytical approximation to in vivo data: Modeling Residue Function for DSC-MRI[END_REF]. The time-rates τ F , τ S were varied together in percentage to render the desired value of MTT = f /τ F +(1 -f ) /τ S . The concentration time-curves are then converted into the corresponding signal intensities via eq. (3.25). The transverse relaxivity κ is fixed according to Østergaard et al. [1996]. For the arterial and tissue signals, the pre-injection baseline signal S 0 and the echo-time T E are chosen according to [START_REF] Wirestam | Wavelet-based noise reduction for improved deconvolution of time-series data in dynamic susceptibility-contrast MRI[END_REF]. Signals are obtained to simulate a repetition time T R = 1s and corrupted by Gaussian noise with zero mean and standard deviation σ = S 0 /SN R before being reconverted into concentration timecurves. Table 12.1 reports the values of the parameters used for the experiments. Presented results are obtained for SN R = 50 [START_REF] Boutelier | Bayesian hemodynamic parameter estimation by bolus tracking perfusion weighted imaging[END_REF], Mehndiratta et al., 2014a], which corresponds to an effective baseline SNR on the tissue signal around 16.7.

Table 12.1: Parameters used for in silico experiments.

t 0 γ 0 ν ξ f (κ • T E) a S a 0 (κ • T E) ts S ts 0 T R 10s

In vivo data

The MR perfusion data was acquired from a 65-years-old female with left M1 distal occlusion three hours after onset. A total of 10 ml of Gadolinium was administrated with rate 3 ml/s, with 20 ml of saline chase at the same speed. A 256 × 256 × 15 volume with 60 time samples was acquired with a 1.5T Siemens device using a Spin-Echo sequence with T E = 52 ms and T R = 1.5 s. The ratio between arterial and tissue transverse relaxivities is κ arterial /κ tissue = 0.04 [START_REF] Boutelier | Bayesian hemodynamic parameter estimation by bolus tracking perfusion weighted imaging[END_REF][START_REF] Rohrer | Comparison of magnetic properties of mri contrast media solutions at different magnetic field strengths[END_REF].

In order to give an idea of data quality, an estimation of the SNR is provided. This is based on the baseline signals, i.e. before tracer arrival, and is expressed as SNR = E[S a 0 ]/σ ts 0 , where E[S a 0 ] is the expectation of the selected arterial baseline signal and σ ts 0 is the calculated standard deviation of the baseline tissue signal in each voxel. Figure 12.1 shows the SNR map for a slice of interest (left) and the histogram of the SNR values for the whole acquired volume (right). The mean SNR is above 50.

In silico experiments and results

This section presents the ensemble of experiments performed to validate the proposed method based on DCB deconvolution and the related results. It first presents a comparison in the absence of dispersion between DCB, oSVD, and the state-ofthe-art CPI (section 12.4.1). After, the influence of dispersion on the results obtained with non-parametric methods is studied, particularly showing the drawbacks of oSVD (section 12.4.2). The section continues with an extensive comparison of DCB with oSVD and CPI+VTF deconvolutions, including all the dispersion kernels and levels (section 12.4.3). Finally, as an application of DCB deconvolution, it is here proposed and validated its use as a pre-processing step to improve the

Dispersion-free comparison with oSVD and CPI

This section reports the comparison between DCB with oSVD and CPI deconvolution methods in the absence of dispersion and delay. Here the comparison is with CPI because it is designed for dispersion-free data and performs better than its dispersion-compliant evolution CPI+VTF [Mehndiratta et al., 2014a]. In fig. 12.2 CBF (first row) and MTT (second row) estimates are compared for two scenarios. In the left images (a,b,c,d ), the ground-truth values ("GT" tag in the bottom) are fixed to render a low cerebral blood flow (CBF = 10 ml/100g/min in a,c) and high mean transit time (M T T = 21s in b,d ) regime. In the right images (e,f,g,h) the GT values (CBF = 30 ml/100g/min in e,g and M T T = 4s in f,h) comply with a normal tissue scenario. On top of each column it is reported the perfusion parameter that is varied in the x-axis. A summary of the results follows.

Results

The results in fig. 12.2 show that DCB (blue) and CPI (yellow) perform overall better than oSVD (purple). For instance, results in normal tissue scenario obtained with oSVD confirm the well-known CBF underestimation (e,f ) and MTT overestimation (g,h). However when the actual SNR is low, oSVD (also aided by the optimal setting described in section 12.1) helps the estimation. Indeed, in the left side of the image (a) oSVD renders the best estimates of CBF, whereas DCB has the highest standard deviation. This is due to the very low effective signal-to-noise ratio -low CBF, low MTT and SNR = 50 -which penalizes the higher flexibility of DCB compared to CPI and to the robustness to noise of oSVD. On the other hand, oSVD has the worst performance in MTT estimation.

In general, CPI and DCB have similar performance. However, in normal tissue scenario DCB reduces the standard deviation in CBF estimation as MTT increases (e), and shows the best recovery of MTT as the CBF increases (h). Moreover, in the low CBF and high MTT regime (a,b,c,d ) a pathological scenario can be identified for low CBF and prolonged MTT, that is the right side of images (a,c) and the left side of images (b,d ). In this condition, DCB shows the best attachment to the ground-truth (dashed lines). Overall, DCB deconvolution performs better or at least comparably to CPI on dispersion-free data, whereas both of the techniques improve results compared to oSVD. This result is relevant for DCB: in fact, it is more general than CPI since it is designed to handle both dispersion-free and dispersed scenarios.

Estimation error as dispersion increases

Influence of dispersion on DCB and oSVD results

This section compares the two non-parametric approaches to deconvolution, DCB and oSVD, in the presence of increasing levels of dispersion. This experiment employ the EDK dispersion kernel varying the vascular mean transit time (M T T v = 1/β) in range M T T v ∈ [START_REF] Marco | Elucidating Dispersion Effects in Perfusion MRI by Means of Dispersion-Compliant Bases[END_REF]10]s. The ground-truth R * (t) was calculated with eq. (10.6). Signals are generated with time-step ∆t = 0.1s and then down-sampled to obtain a repetition time T R = 1s. The tested delay values are 0, 1, 2, 3, 4, 5 seconds [Mehndiratta et al., 2014a] but for DCB the delay is researched in a broader range [-2, 7]s with τ s = 1s. Other tested parameters are M T T ∈ {4, 8, 12, 16}s and CBF ∈ {15, 30, 45, 60} ml/100g/min. The comparison is performed by calculating the mean error and standard deviation over 1000 noisy realizations. In fig. 12.3, the left and the central images report the relative errors for CBF * and MTT * , and the right image reports the absolute error for the time of maximum t max of the effective response function R * (t). The right image additionally reports the absolute error of the dispersion time τ δ calculated with DCB deconvolution.

Results

The DCB results show overall an improved stability compared to oSVD. The oSVD estimation error for all of the tested parameters in fig. 12.3 increases with the vascular mean transit time M T T v , i.e. the dispersion level. On the contrary DCB performance for CBF * improves with M T T v (fig. 12.3, left). The results related to the time parameters in fig. 12.3 (right) show that average absolute error for the estimation of t max falls below 1s which is sensibly lower than oSVD error at any M T T v . The image also shows the absolute error for the dispersion time τ δ obtained with DCB deconvolution. The dispersion time is a parameter derived from the joint estimation of the time of maximum t max and the time delay τ , therefore the error on τ δ accounts for both of the error contributions.

The error on τ δ estimation is also indicative of the performance of DCB deconvolution in characterizing dispersion. Note that the DCB mean absolute error falls approximately within 1s which corresponds to the temporal resolution, i.e. the repetition time T R, used for the experiment. The oSVD technique scores the best performance in CBF * estimation at low M T T v , likely because of its implementation as discussed in section 12.1, but renders a higher error for MTT * .

Results globally show that the proposed deconvolution with DCB is more robust, with respect to dispersion, than with oSVD. Moreover, the results give an indication of the bias in the estimation of the effective parameters with oSVD as dispersion increases. The next section reports a more extensive experiment, to compare all of the dispersion-compliant techniques (DCB, oSVD, and CPI+VTF), accounting for the other dispersion kernels in eq. ( 10.5), i.e. GDK and LNDK, and for the dispersion levels encoded in table 10.1.

Comparison with oSVD and CPI+VTF in dispersed scenario

This experiment compares the performance of DCB, oSVD and CPI+VTF in estimating the parameters CBF * , MTT * , t max , τ , τ δ and the fitting error. Ground-truth signals are generated as described in section 12.4.2. However, this experiment takes into account dispersion. Particularly, the data is generated taking into account all of the three dispersion kernels (GDK, EDK, and LNDK), as described in section 12.2 and illustrated in fig. 10.1, and the three dispersion levels (low, medium, and high) reported in table 10.1. The DCB, oSVD, and CPI+VTF deconvolutions are performed for 100 different noisy realizations for each combination of CBF, MTT, τ , dispersion kernel and dispersion level. Results in figs. 12.4 and 12.5 are presented in the form of error boxplots 2 .

2 Each boxplot is a comprise of lower and upper whiskers within which all of the results not considered as outliers are contained. The bottom and top horizontal sides of the box indicate the 25 th and 75 th percentile of the results respectively. The central horizontal line of the box indicates the median whereas the height of the notch-like concavity in its correspondence describes the confidence region: boxes related to different techniques with non-overlapping confidence regions along the vertical axis show medians with 5% significant difference. The percentage of outliers is reported on top of the upper whisker.

In the first part of the experiment, the estimates of CBF * , MTT * , t max , τ , τ δ , and the fitting error are computed when the ground-truth data is generated with all of the dispersion kernels but at different dispersion levels. Images in fig. 12.4 report the error boxplots for each estimated parameter and for each, when applicable, deconvolution technique adopted (oSVD, DCB, CPI+VTF). Each row from the first to the fifth corresponds to a different estimated parameter whereas the sixth row reports the reconstruction error with respect to the effective response function. The first column starting from the left of fig. 12.4, framed within a black rectangle, shows error boxplots obtained when considering results at all of the dispersion levels. Columns from the second to the fourth are related to low, medium and high levels of dispersion respectively.

In the second part of the experiment, the estimates are computed when the ground-truth data is generated considering all of the dispersion levels but with different dispersion kernels. Columns from the first to the third of fig. 12.5 show results when the underlying vascular transport function (VTF) corresponds to GDK, LNDK and EDK respectively.

Results

Deconvolution with DCB performs globally better than the compared techniques in the presence of dispersion, as shown in the first column of fig. 12.4. Indeed, for the tested parameters MTT * , t max , τ and τ δ the median of the blue boxplot is lower than the comparing ones, with 5% significant difference. Also the whiskers and horizontal sides (25 th and 75 th percentiles) generally compare favorably for DCB. An exception to the trend is the great performance of oSVD for CBF * estimation, noticed also in Mehndiratta et al. [2014a], which is counterbalanced by a poorer performance w.r.t. MTT * . However, here the oscillation index for oSVD is chosen optimally (see section 12.1), which explains the result at least in part. Overall DCB scores the best results even if it generally leads to a higher percentage of outliers (number on top of each boxplot). CPI+VTF performs comparably to oSVD for the estimation of MTT * and outperforms it with respect to t max .

Results obtained with the three techniques show different trends as the dispersion level increases. The trend for the same technique also changes based on the considered parameter, i.e. each row of fig. 12.4. For instance, the time of maximum t max estimation performance of both oSVD and CPI+VTF improves with the amount of dispersion, whereas DCB results look more stable and show lower or comparable errors. DCB shows a similar stability also for the effective mean transit time MTT * .

Results for delay τ and dispersion time τ δ (fourth and fifth row of fig. 12.4) look similar as expected. Indeed, for both of the parameters the median error of DCB (blue boxplots) is lower than that of CPI+VTF (red boxplots). DCB error In the second part of the experiment, the results are illustrated separately for each adopted dispersion kernel (GDK, LNDK, and EDK) while accounting for all of the dispersion levels, as shown in fig. 12.5. Again, similar considerations as before apply to oSVD results, which reflect the trade-off between CBF * and MTT * estimation. With any dispersion kernel, results for effective mean transit time MTT * and time of maximum t max show that DCB deconvolution outperforms oSVD and CPI+VTF. Interestingly, for these parameters DCB performs better than or comparably to CPI+VTF also when data is generated using the gamma dispersion kernel (GDK), which is the kernel inherently assumed by the VTF model-based technique. When considering the gamma dispersion kernel (GDK), DCB deconvolution shows a performance loss in delay τ and dispersion time τ δ estimation. On the other hand DCB remarkably improves τ and τ δ estimates for data based on lognormal (LNDK) and exponential (EDK) kernels. Fitting errors generally favor DCB deconvolution over CPI+VTF with the expected exception of the GDK case, where the difference is not significant (last row of fig. 12.5).

Overall, the median errors reveal that DCB deconvolution performs better than CPI+VTF and oSVD when considering all the dispersion kernels at any dispersion level (fig. 12.4). Results for the single dispersion kernels (fig. 12.5) show parameterspecific tendencies, such in the case of delay and dispersion time for the GDK, but generally confirm the positive performance of DCB deconvolution.

DCB as pre-processing step for CPI+VTF

The CPI+VTF technique, since it relies on a model of the VTF, offers the advantage of estimating the dispersion-free CBF and MTT directly in the deconvolution procedure. In this thesis, it is proposed instead to estimate these parameters in two steps according to conditions expressed in algorithm 1 and the pipeline shown in fig. 11.4: first, we obtain R * (t) via DCB deconvolution; second, we perform CPI+VTF model fitting directly on the recovered R * (t). Experiments are performed to show the advantage of estimating CBF and MTT with this second proposed method. The comparison is performed when the ground-truth data is generated accounting for all of the dispersion kernels at different dispersion levels (fig. 12.6), and for all the dispersion levels with different dispersion kernels fig. 12.7, as described in section 12.4.3.

Results

Results in fig. 12.6 show that both DCB and CPI+VTF globally render similar errors of the dispersion-free cerebral blood flow CBF (first row and column). Although CPI+VTF reduces the error at low dispersion, globally and at medium and high dispersion DCB shows the best performance. Indeed, the error with CPI+VTF increases with the level of dispersion. In addition, a previous deconvolution with DCB followed by the CPI+VTF model fitting considerably reduces the estimation error of the mean transit time MTT at all the dispersion levels (globally a 50% error reduction), as shown in the first column of fig. 12.6. Also in this case, the amount of error with CPI+VTF increases with the dispersion level. The comparison that considers different dispersion kernels (GDK, LNDK, EDK), shown in fig. 12.7, brings evidence of the benefit of pre-processing data with DCB for the joint estimation of CBF and MTT. Overall, the pre-processing with DCB allows a clear improvement in the recovery of dispersion-free parameters, for all dispersion levels and kernels.

In vivo experiments and results

A slice of the real dataset is manually segmented into healthy isoperfused and infarcted hypoperfused regions. The segmentation is reported in the top left corner of fig. Images in fig. 12.9 report the CBF, CBV, M T T and t max maps obtained with the three tested methods. Maps of delay τ and dispersion time τ δ for DCB and CPI+VTF (and absolute difference) are shown in fig. 12.10. In fact, these maps are not supported by oSVD deconvolution. Finally, fig. 12.11 reports the maps of actual dispersion-free cerebral blood flow CBF and mean transit time MTT obtained with CPI+VTF model fitting on the response function obtained with DCB (first column) and with CPI+VTF deconvolution directly (second column). The absolute difference maps are reported in the third column. In the case of DCB the model fitting is performed only on voxels where the dispersion time is not null (τ δ > 0), that is only for the voxels where dispersion is detected. Indeed, in non-dispersed voxels the effective perfusion parameters are also the actual dispersion-free ones, i.e. CBF = CBF * and M T T = M T T * . In figs. 12.9 to 12.11, large vessels were removed using a threshold on CBV values, i.e. voxels with CBV > 3% are set to zero [START_REF] Boutelier | Bayesian hemodynamic parameter estimation by bolus tracking perfusion weighted imaging[END_REF], whereas cerebrospinal fluid (CSF) was previously removed. No subsequent filtering or smoothing was applied.

A quantitative summary of the results, for the iso-and hypo-perfused regions shown in the top left image of fig. 12.8, is reported as histograms. Particularly, fig. 12.12 and fig. 12.13 report the histograms of the estimated CBF, CBF * , MTT and MTT * for DCB, CPI+VTF, and oSVD when applicable, for the isoperfused and hypoperfused regions respectively.

General Results

In general, all the techniques show decreased values of CBF * (first row of fig. 12.9) in the infarcted hypoperfused region (right side region in top left image of fig. 12.8) and corresponding increased MTT * . The maps shown fig. 12.9 are qualitatively similar, but oSVD shows less contrast between the regions for both of the parameters. This is quantified in the histograms reported in figs. 12.12 and 12.13. The maps of dispersion time τ δ in fig. 12.10 show that DCB detects none or little dispersion in the healthy region, whereas CPI+VTF renders dispersed response functions almost everywhere. A more direct evidence of this is given by the shape of the mean residue functions in the bottom left image of fig. 12.9 voxels with CBV > 3% were set to background to remove the contribution of large vessels whereas CSF was previously removed. Contrary to CPI+VTF, the DCB map of dispersion time displays several voxels with τ δ = 0 implying absence of dispersion (see algorithm 1). For DCB, dispersion-free voxels are mainly concentrated in the left hemisphere, contralaterally to the infarcted region.

However, the DCB-based estimates of MTT within the hypoperfused region are more compatible with a dispersion-free scenario. Indeed, after removing the bias due to dispersion, the dispersion-free MTT should not be larger than the apparent effective one, i.e. M T T ≤ M T T * . In the following, results are discussed in detail.

Cerebral Blood Flow

Maps of CBF * obtained with oSVD look smoother but seem to be visually less discriminating than those obtained with DCB and CPI+VTF (first row of fig. 12.9). Indeed, a comparison between the CBF * histograms in the first row of fig. 12.12 and fig.

12.13 reveals that the change in width of the oSVD distribution between iso-and hypo-perfused regions is smaller than with the other techniques. In fact, both DCB and CPI+VTF show a more evident change in spread with a support extending up to 100 ml/100g/min in the isoperfused Region Of Interest (fig. 12.12). The oSVD shows the lowest median and the smallest spread whereas CPI+VTF shows the highest median and largest spread. The histogram of DCB presents an intermediate situation. In general, CBF * histograms for oSVD, DCB, and CPI+VTF are in better agreement with each other in the hypoperfused ROI (fig. 12.13) than in the healthy isoperfused one.

Mean Transit Time

All of the methods reveal that MTT * increases in the hypoperfused area (third row of fig. 12.9) in accordance with the infarcted condition of the tissue. However, the DCB and CPI+VTF maps show the greatest contrast between healthy and infarcted regions. This is supported by the MTT * histograms in figs. 12.12 and 12.13 (third row), where the ratio of the medians between hypo/iso-perfused regions is approximately 8.5s/2s for DCB and CPI+VTF, and only 7.5s/4s for oSVD. Note

Time of maximum, delay and dispersion time

Maps of the time of maximum t max of the effective response function, shown in the last row of fig. 12.9, reveal to be discriminating between the iso-and hypo-perfused ROIs with all of the tested techniques. The t max values obtained with DCB visually look smaller than with oSVD, which typically overestimates (see results in figs. 12.3, 12.4 and 12.7). In the hypoperfused region oSVD and CPI+VTF maps look visually similar whereas the DCB map shows a reduced amount of high estimates (red areas). In the isoperfused region both oSVD and DCB show voxels where t max = 0 whereas CPI+VTF almost never renders such a value. This is reflected also in the dispersion time τ δ maps of fig. 12.10 (second row), where CPI+VTF finds τ δ > 0 almost everywhere. This phenomenon may constitute a bias for CPI+VTF because it implies that every voxel in the healthy region presents an effective response function with a dispersed shape. The mean residue function recovered with CPI+VTF, shown in the bottom left image of fig. 12.8, confirms this "dispersion-detecting" trend. Indeed, DCB renders a dispersion-free shape in the healthy tissue whereas CPI+VTF shows a dispersed profile. It is interesting to note that in the infarcted area DCB shows a smoother and more dispersed profile, which is in better agreement with oSVD than with CPI+VTF (bottom right image). Note that oSVD solution, when the response function is smooth as in the case of high dispersion, is a good solution since the high-frequency cutoff likely removes only noise and not discontinuities, as opposed to what would happen in the case of an exponential response function.

The delay τ map for CPI+VTF in fig. 12.10 also depicts no variability within the isoperfused region, whereas DCB presents positive delays more frequently. In addition, DCB seems to render lower delays in the hypoperfused ROI, which has to be ascribed to corresponding lower values of t max . In the isoperfused region, DCB seems to split the time of maximum contributions more on the delay than on the dispersion time compared to CPI+VTF. Globally, oSVD seems to overestimate t max . To summarize, DCB renders in average a dispersion-free response function profile within the healthy region, as opposed to oSVD and CPI+VTF, and a smooth dispersed profile within the infarcted one.

Dispersion-free parameters

The maps of the actual dispersion-free perfusion parameters in fig. 12.11 obtained with DCB and CPI+VTF show similar contrast between healthy and infarcted ROIs, but present substantial differences. For instance, the CPI+VTF map of CBF reports higher values than the corresponding DCB ones in the isoperfused region (see segmentation in fig. 12.8). This is confirmed in fig. 12.12: CBF histogram of DCB shows a higher density for values lower than 50 ml/100g/min, whereas with CPI+VTF the histogram spreads and presents a higher density of values above 75 ml/100g/min. On the other hand, CBF histograms (DCB and CPI+VTF) in the hypoperfused region bear similar information (fig. 12.13). The MTT maps in fig. 12.11 of DCB and CPI+VTF present differences particularly in the hypoperfused region, where the DCB map shows lower values. To clarify, in the presence of dispersion -detected by both the techniques as shown in the second row of fig. 12.10 Also in real data, DCB deconvolution managed to correctly recover voxels with zero dispersion time, which implies voxels where no dispersion occurs. These nondispersed voxels are mainly concentrated in the healthy isoperfused region, where CPI+VTF renders τ δ > 0 almost everywhere instead. Consequently, all the voxels processed with CPI+VTF have to be interpreted as dispersed, as demonstrated by the corresponding shape of the mean residue function in the bottom left image of fig. 12.8. This might be the cause of the higher CBF values, observed within the healthy region, when compared to DCB deconvolution (maps in fig. 12.11 and histograms in fig. 12.12). In fact, the pre-processing with DCB deconvolution offers the advantage of selecting the voxels affected by dispersion -when τ δ > 0 -so that we can avoid fitting a VTF model based technique (such as CPI+VTF or others) in dispersion-free voxels, thus reducing the risk of overfitting. Conceptually, if the reported DCB dispersion time map in fig. 12.10 is correct, as suggested by the synthetic experiments in fig. 12.4, it means that CPI+VTF detects dispersed shapes where there is no dispersion, which might lead to a consequent overestimation of the dispersion-free CBF (see the absolute difference map in fig. 12.11). To this regard, in this work CPI+VTF was not implemented with the original bayesian variational method as reported in [START_REF] Chappell | Variational Bayesian Inference for a Nonlinear Forward Model[END_REF], and this constitutes a limitation of our study. Indeed, here CPI+VTF is implemented without the a priori information elaborated in Mehndiratta et al. [2014a], which is designed to bias the estimation of the effective response function towards a dispersed shape with high probability for low-to-medium level of dispersion. However, it would appear that such a bias can further emphasize the "dispersion-detecting" behavior of CPI+VTF in the healthy region.

The computational cost of DCB and CPI+VTF is higher than oSVD. However for a comparable number of parameters adopted, DCB is significantly faster than CPI+VTF thanks to its analytic Jacobian. The atomic cost per iteration for each voxel is 0.014s for oSVD, 0.67s for DCB and 0.83s for CPI+VTF on an Intel Core TM i7-3840QM 2.80GHz CPU, with 32GB RAM. However, the techniques are suitable for code optimization and parallel computing.

The proposed Dispersion-Compliant Bases method (DCB) allows explaining data without assumptions about the presence, amount, or nature of dispersion. The method reconstructs plausible response functions (fig. 12.8) in healthy and pathological scenarios. Moreover, it globally renders better estimates of the effective perfusion parameters, and its use as a pre-processing step boosts the precision of dispersion-free estimates of CBF and MTT.

Conclusion

This chapter reported an extensive validation of the proposed DCB deconvolution by means of comparison with oSVD, and the state-of-the-art CPI, and CPI+VTF. The test in silico/vivo demonstrated the validity of the proposed approach which makes no assumptions about the presence, nature, and amount of dispersion. Moreover, pre-processing data with DCB deconvolution demonstrated to improve results up to 50% when estimating the dispersion-free CBF and MTT.

The final next chapter 13 reports a more investigative approach about the biology of the tissue. The aim is to infer which nature of the dispersion process is better supported by data, that is what model for the vascular transport function is better suited for a specific voxel-region. To accomplish this task, the linearized version of the dispersion-compliant bases will be adopted, together with newly developed variants of control point interpolation "+VTF": in addition to the CPI+VTF, that corresponds to CPI+GDK, the newly implemented variants +EDK and +LNDK will be employed to find the best VTF, on a voxel basis, in synergy with a L-DCB deconvolution.

Chapter 13

Unvealing the dispersion kernel with DCB

Overview

Once obtained the estimates of the effective perfusion parameters, CBF * and MTT * , one additional step is to recover the dispersion-free parameters in the voxels where dispersion is present. However, this is not an easy task because in order to do so, first it is necessary to detect the presence of dispersion, and second it is mandatory to chose a model for the vascular transport function to be used to perform the deconvolution R(t) = R * ⊗ -1 VTF(t).

As an additional contribution of this thesis, this chapter, mainly based on the work of Pizzolato et al. [2016e,c], proposes a solution to this problem. The proposed linearized dispersion-compliant bases (L-DCB), opportunely validated, are then used to perform deconvolution of C a (t) and C ts (t) in order to obtain the effective response function, R * DCB (t). The proposed dispersion-time is then used to detect the presence of dispersion and an algorithm is proposed to decide which dispersion kernel, among the exponential (EDK), gamma (GDK), and lognormal (LNDK), is the one that better describes the data for a particular voxel. To do so, variants of CPI+ technique, i.e. +EDK, +GDK, and +LNDK, are newly implemented, and their results are interpreted in the light of those obtained with DCB. The idea is that the deconvolved R * DCB (t) is general and model-independent: therefore it can be used as reference to see which VTF model better describes it.

Keywords dispersion kernel; model selection; L-DCB; VTF.

Dispersion kernel detection: motivations

There are two main reasons why it is important to detect which dispersion process the injected bolus of tracer undergoes.

The first one is biological investigation. Different shapes of the VTF reflect different distributions of the vascular transit times, as discussed in section 10.3, thus being potentially informative of the non-representative behavior of the selected healthy arterial input function (AIF) with respect to a specific pathological voxelregion. In fact, the shape of the VTF can be inherently considered a biomarker.

The second reason is methodological. Assuming a specific model for the VTF can be reductive. indeed, this can lead to a wrong estimation of CBF, MTT, and of the model parameters of the selected dispersion kernel itself. The intuitive side behind this affirmation is that if the true unknown VTF is, for instance, a LNDK, then the use of a GDK model causes an implicit approximation error. More importantly, the non-linear routine will try to estimate R ⊗ VTF(t) such that the reconstruction error falls below a specific threshold. This implies that all the approximation error due to a possible wrong choice of the VTF model will be reflected, at convergence, into the estimated R(t) and parameters Ĉ BF, M TT. While the ultimate goal is model-independent estimation of both R(t) and VTF(t), this chapter performs one step further compared to literature by proposing a strategy for in vivo VTF model-selection based on model-independent DCB deconvolution.

CPI+ techniques

The strategy adopted to perform the deconvolution R(t) = R * ⊗ -1 VTF(t) consists on using the CPI+GDK/EDK/LNDK techniques such that R(t) is estimated in a model-free fashion and the VTF is associated to a model. The variants are here implemented, similarly to what described in section 11.3, by changing the dispersion kernel, as depicted in fig. 13.1.

The techniques rely on a non-linear routine. Similarly to CPI+VTF, the parameters of the different dispersion kernels have been initialized according to the low dispersion values seen in table 10.1. The L-DCB deconvolution was performed with order N = 6.

In order to use the DCB solution for R * (t) as a reference, it is necessary to demonstrate that the performance of DCB deconvolution is superior or comparable to that of each CPI+ technique when the ground-truth dispersion kernel used for in silico experiments is that inherently assumed by the compared CPI+ technique itself. Therefore, by using the same methodology for experiments specified in section 12.2, a performance evaluation follows. Reported results show the good performance of L-DCB compared to the CPI+ techniques. This opens for the possibility of using DCB deconvolution as reference for voxel-wise VTF detection. The next section reports the proposed strategy to achieve this goal.

In vivo dispersion detection and VTF identification

It is proposed to apply DCB deconvolution on stroke MRI data to estimate the effective response function, R * DCB (t), for each voxel. The sought solution is then used as a filtering technique to reduce the illposedness of the deconvolution process. In practice, this is achieved by calculating the re-convolution Ĉts (t) = C a ⊗ R * DCB (t) where the symbol "hat" distinguishes the re-convolution from the acquired tissue concentration C ts (t). The idea is to perform deconvolution between the measured AIF, C a (t), and the DCB-filtered Ĉts (t) with the CPI+ techniques in order to verify, according to certain criteria, which of the adopted dispersion kernels, among GDK, EDK, and LNDK, better allows recovering the solution R * DCB (t) obtained with DCB. In fact, from the three CPI+ techniques the effective response functions R * GDK (t), R * EDK (t), and R * LN DK (t) are computed. The procedure is schematized in fig. 13.3. The idea is that each of the VTF model-based developed CPI+ techniques should not, in this case, be directly fitted to the response function obtained with DCB, but it should rather be let free to explain the data via deconvolution knowing that the actual convolution kernel is the one obtained, in a model-free fashion, via DCB. This procedure allows having a reference for the effective response function in vivo, such that a data-driven choice can be made despite the lack of biological validation. The choice of DCB as reference is justified by its higher performance in synthetic results and by its lack of assumptions about a VTF model: in this sense, DCB can explain data where the underlying VTF is any of the considered ones. Then, for each voxel it is selected the CPI+ technique that scores the lowest l 2 reconstruction error of Ĉts (t), and the best estimates of the effective blood flow CBF * and mean transit time MTT * using DCB results as reference. The algorithm is illustrated in fig. 13.4.

Results

Results are visualized in the form of a colormap where each voxel of the selected dataset, that corresponds to that described in section 12.3, is colored according to the detected vascular transport function. The color codes indicated in fig. 13.4 are red for EDK, green for GDK, and blue for LNDK.

In the first experiment, whose results are shown in the first row of fig. 13.5, the possibility of absence of dispersion is not considered. This means that every voxel is considered dispersed and the information coming from DCB deconvolution, such as the possibility that τ δ = 0, is not taken into account. Each dispersion kernel can, theoretically, assume a Dirac delta shape for specific limits of the parameters, such as expressed in eq. (11.12). However, this limit remain in practice only theoretical due to the impossibility of reaching those limits. Nevertheless, the images of the best for best for best for l 2 error best for "all" GDK EDK LNDK none only DKs using (second), the lowest l 2 error norm (third), or the intersection of all the previous conditions (fourth column). Images in the first row report the maps obtained when only considering the selection of the kernels, wheres maps in the second row depict a purple pixel when dispersion is not detected according to the dispersion time (τ δ = 0) calculated on the effective response function obtained via L-DCB. Not that a gamma kernel is mainly present in the hypoperfused region, whereas an exponential kernel is mainly localized elsewhere despite its presence largely correlates with the absence of dispersion: red pixels in the first row are substituted by purple pixels in the second row. first row of fig. 13.5 report the color maps of voxel-wise detected dispersion kernel where the selection is based considering the optimal reconstruction of DCB results with respect to CBF * (first column), MTT * (second column), l 2 fitting error norm (third column), and the intersection of all of these (last column).

It appears clear that the brain tissue, based on the proposed data-driven analysis, supports different models for the vascular transport function. It is interesting to note that the topographic distribution of EDK, GDK, and LNDK is very correlated to the manual segmentation of healthy/infarcted region presented in the previous chapter in fig. 12.8. Not surprisingly it is also correlated with the maps reported in fig. 12.9. Indeed, GDK (in green) is concentrated in the infarcted region where dispersion and hypoperfusion were detected. On the other hand, EDK is concentrated in all the rest of the brain tissue where none or little dispersion were detected. This can be interpreted in the light of the characteristics of the exponential and gamma models. Particularly, the convolution R ⊗ EDK(t, β) leads more easily to sharper effective response functions than R ⊗ GDK(t, p, s). This is confirmed in the second part of the experiment where, this time, the information provided by the dispersion time, such as τ δ = 0, is used.

The second part of the experiment considers a fourth category of voxels: this contains those voxels where no dispersion is detected according to τ δ = 0. The corresponding images are reported in the second row of fig. 13.5, where non-dispersed voxels are colored in magenta. This second result shows that for all the comparison criteria, i.e. the different columns of the figure, the absence of dispersion perfectly correlates with an exponential kernel, whereas the presence of dispersion is almost everywhere explained by a gamma kernel. This second affirmation finds agreement with the work of Chappell et al. [2013].

As a disclaimer, here results are the outcome of a data-driven classification based on dispersion-compliant bases deconvolution. The underlying biological interpretation is beyond the scope of this thesis. Results rather show that the use of a specific VTF model might not be justified in DSC-MRI deconvolution. In fact, different brain regions support different optimal dispersion kernels. Particularly evident is the case where no dispersion is detected via DCB deconvolution. In this case, data reveal a clear preference for an exponential kernel whereas in the presence of dispersion the preference is mainly for the gamma one, and for the lognormal.

On the use of models to characterize dispersion

Results disclosed in this chapter suggest that DSC-MRI deconvolution should be performed with the least amount of assumptions about the presence of dispersion and, even more, its nature. This translates into using approaches, like the proposed dispersion compliant bases, that avoid the use of models for describing the shape of the vascular transport function when performing deconvolution. However, the use of models is still important for dispersion quantification and, consequently, for the estimation of the dispersion-free CBF and MTT. This thesis promotes a model selection strategy where after performing deconvolution with DCB, thus avoiding the issue of a model-based bias, the presence of dispersion and the eventual VTF model are identified a posteriori.

The open aim remains the full non-parametric quantification of the residue function, R(t), and VTF. However, this is an extremely challenging task because of the difficulty in finding appropriate conditioning and constraints for the estimation.

Conclusion

This concluding chapter on perfusion MRI reports a twofold contribution. In first instance, it provides evidence in vivo of the fact that the brain tissue supports different dispersion kernels, or rather, it supports dispersed and non-dispersed voxel-regions. This underlines even further the importance of performing perfusion deconvolution in a model-free fashion and in a purely data-driven manner like in the proposed DCB technique. Secondly, this chapter proposes a strategy to avoid biasing the deconvolution procedure while, at the same time, allowing the detection of dispersion and, eventually, the identification of the dispersion kernel model to adopt for representing the vascular transport function. This opens for a hybrid DSC-MRI deconvolution where, when a voxel-region is dispersed, the correct model can be selected to recover the dispersion-free regime.

Conclusion and perspectives

211

This thesis focuses on computational aspects related to Diffusion and Perfusion Magnetic Resonance Imaging (dMRI and pMRI). These modalities are complementary to each other in brain imaging. In fact, on one side dMRI allows estimating features of the signal that are related to the tissue microstructure. On the other side, signal's features in pMRI are related to the blood passage, such as the cerebral blood flow, in the functional part of the tissue.

In this thesis, challenges in both dMRI and pMRI are tackled from a theoretical and methodological point of view, with the final objectives of introducing new signal-based features characterizing the tissue, and improving the quality of data processing.

In part I, the basics of the diffusion and perfusion signal contrasts were summarized, showing similarities in the underlying physical phenomena.

In part II, the focus was on diffusion MRI. This part has generally a twofold contribution: understanding the nature of the complex diffusion signal, and exploiting it for the interpretation of the diffusion process, such as the diffusion asymmetry, and for denoising purposes, such as phase correction.

After describing the basics of diffusion MRI in chapter 4, this part poses the question of whether the complex nature of the signal attenuation could be exploitable for inferring properties of the underlying diffusion process and tissue microstructure. Examples have been given in chapter 5 and new insights have been proposed in chapters 6 and 7, where with different paradigms, such as the axonal compression and temperature gradients, the characteristics of the potential diffusion asymmetry have been explored and characterized.

In chapter 8 the theoretical investigation has left place to the exploitation of the complex signal for denoising purposes via phase correction. With this procedure, the Rician bias that affects magnitude DWIs can be avoided. In this case, one first contribution of the thesis is to carry on a rigorous characterization and quantification of the effects of phase correction on debiasing diffusion signal and scalar metrics, such as those calculated with DTI and MAP-MRI (q-space).

Another contribution is the definition of a strategy to find the optimal amount of phase correction needed for each DWI. This helps avoiding fine tuning of parameters that could end up in a corruption rather than a correction of the DWIs. The result is a general framework to perform phase correction, which also takes into account another important characteristic of the noise, i.e. its spatial variability. The proposed framework, validated on real and synthetic data, constitutes a complete Off-the-Shelf tool for DWI pre-processing.

The contributions of the thesis, for the diffusion part, are both theoretical and methodological, posing questions and giving practical solutions to the diffusion MRI community.

Part III tackled challenges related to perfusion MRI. This part focuses on the development of a framework to perform perfusion deconvolution. Although this could be applicable to other perfusion MRI modalities, here the focus is the treatment of DSC-MRI data, accounting for the presence of bolus dispersion.

After proposing an analysis of dispersion effects in chapter 10, a method based on dispersion-compliant bases (DCB) is proposed in chapter 11 with the purpose of estimating the response function (which characterizes perfusion) also in the presence of dispersion being, in this sense, dispersion-insensitive. At the same time, the framework allows solving the issues discussed in chapter 9. In fact, the proposed DCB deconvolution, developed in chapter 11,

• allows recovering physiologically meaningful response functions;

• does not involve any explicit modeling;

• allows estimating the time-delay between the arterial and tissue concentrations;

• implements non-negativity constraints for the sought response function;

• reliably allows estimating the proposed dispersion time to detect and characterize the presence of dispersion;

• improves the performance of model-based techniques to recover dispersion-free perfusion parameters, i.e. CBF and MTT.

These achievements were validated on synthetic and real data in chapter 12 Additionally, the thesis poses some conceptual questions about whether the presence and nature of dispersion, i.e. the type of vascular transport function (VTF), can be identified in vivo. These questions are tackled by the method proposed in chapter 13, based on DCB deconvolution, that finds on a voxel basis the optimal dispersion kernel model to be used. Indeed, in the case dispersion was detected the correct model could be used to estimate dispersion-free perfusion parameters.

The contributions the perfusion part, are mainly methodological but at the same time aim at establishing a different, more data-driven, and straightforward way of processing perfusion data in MRI.

In addition to the presented contributions, this thesis poses some grounds for further exploration in both diffusion and perfusion MRI. For instance, the link between the complex diffusion signal and the tissue microstructure can be further explored. Asymmetry may be used as a biomarker, as shown in chapter 6, and efforts should point toward real applications. The proposed phase correction could be used to re-evaluate some clinical correlations between diffusion MRI metrics and pathology. Indeed, debiasing of the metrics could lead to new potential insights that might have been hidden by the bias when not performing any corrections. In perfusion MRI, the proposed DCB method can be used to infer, after deconvolution, on the correct model of dispersion, i.e the VTF. In this way, one could build dispersion-free maps of CBF and MTT by selecting, on a voxel basis, the correct model of VTF to use. In this sense the thesis opens for hybrid perfusion data processing. On this lead, since the use of basis functions revealed useful for representing the response function, in a future work the approach could be extended to represent the vascular transport function, thus completely avoiding explicit modeling for it.

Le travail effectué dans cette thèse a consisté en des contributions computationnelles liées aux modalités d'IRM basées sur des images pondérées en diffusion (IRMd) et sur des images de perfusion (IRMp) par injection de produit de contraste. Pour chacune des deux modalités, des contributions, aussi bien théoriques que méthodologiques, ont été développées pour améliorer la qualité, le traitement et l'exploitation des signaux acquis.

Dans la partie I, les bases des contrastes de diffusion et de perfusion ont été résumées, montrant des similitudes dans les phénomènes physiques sous-jacents.

Dans la partie II, consacrée à l'IRMd, l'accent a été mis sur la nature complexe du signal de diffusion, sur l'importance de l'information de phase et sur son exploitation dans le cadre du processus de débruitage par un processus de correction de phase.

Dans le chapitre 4, qui décrit les bases de l'IRM de diffusion, on s'est d'abord intéressé à la nature complexe de l'atténuation du signal et à la possibilité que la prise en compte de l'information de phase puisse être exploitée pour inférer des informations et propriétés sur le processus de diffusion sous-jacent et sur l'information de microstructure des tissus observés. Des exemples ont pu être étudiés dans le chapitre 5, et de nouvelles contributions ont été proposées dans les chapitres 6 et 7, où des paradigmes, tels que la compression axonale et les gradients de température, ont permis d'explorer et de bien caractériser l'asymétrie de la diffusion.

Dans le chapitre 8, une exploitation du signal complexe a été menée dans le cadre d'un processus de débruitage des images DWI par correction de phase. Avec cette procédure, le problème de la gestion du bruit Rician qui affecte l'amplitude des données DWIs peut être évité tout en menant une caractérisation et une quantification rigoureuse des effets de la correction de phase sur le signal de diffusion, et sur les métriques scalaires calculées avec les modèles (espace-q). Une contribution complémentaire a été la définition d'une stratégie pour trouver la quantité optimale de correction de phase nécessaire pour chaque donnée DWI. Cela permet d'éviter le réglage empirique des paramètres qui pourraient se traduire par une corruption plutôt qu'une correction des données DWIs. Le résultat est un cadre général pour effectuer la correction de phase, qui prend également en compte une caractéristique importante du bruit i.e. sa non-stationarité spatiale. Au final de cette seconde partie dédiée à l'IRMd, les contributions computationnelles proposées et validées sur plusieurs jeux de données synthétiques et réelles, constituent un cadre bien adapté pour le prétraitement des images DWI tout en y apportant des solutions pratiques dans le domaine.

Dans la partie III liée à la l'IRM de perfusion, les contributions ont porté sur le développement d'un cadre méthodologique pour effectuer la déconvolution de perfusion. Bien que cela puisse être applicable à d'autres modalités d'IRM de perfusion, l'accent a principalement porté ici sur le traitement des données DSC-MRI afin de tenir compte du biais dû à la dispersion du bolus.

Après avoir proposé une analyse des effets de dispersion dans le chapitre 10, une méthode axée sur des fonctions de base (DCB) est développée au chapitre 11 afin d'estimer les fonctions de réponse quand il y a une dispersion du bolus. On montre en particulier comment la dispersion qui peut empêcher la correcte estimation des métriques de perfusion, peut aussi donner des informations importantes sur l'état pathologique du tissu cérébral. La méthode de déconvolution développée au chapitre 11 permet en particulier de resoudre les problèmes présentés et discutés dans le chapitre 9, en ce sens qu'elle:

• permet de récupérer des fonctions de réponse physiologiquement significatives;

• n'implique aucune modélisation explicite;

• permet d'estimer le délai entre les concentrations artérielle et tissulaire;

• met en oeuvre des contraintes de non-négativité pour la fonction de réponse recherchée;

• permet d'estimer de manière fiable le temps de dispersion proposé pour détecter et caractériser la présence de dispersion;

• améliore la performance des techniques basées sur des modèles pour récupérer des paramètres de perfusion sans dispersion, i.e. CBF et MTT.

Ces réalisations ont été validées sur des données synthétiques et réelles dans le chapitre 12. De plus, la thèse pose des questions conceptuelles sur la possibilité d'identifier in vivo la présence et la nature de la dispersion, c'est-à-dire le type de fonction de transport vasculaire (VTF). Ces questions sont abordées par la méthode proposée dans le chapitre 13, basée sur la déconvolution DCB, qui trouve pour chaque voxel le modèle de dispersion optimal à utiliser. En effet, dans le cas où la dispersion a été détectée, le modèle correct pourrait être utilisé pour estimer les paramètres de perfusion sans le biais dû à la dispersion. Les contributions de la partie de perfusion sont principalement méthodologiques mais visent en même temps à établir une méthode originale et efficace axée sur les données, et plus simple pour traiter les données de perfusion.

Au final, et au delà des contributions apportées tout au long de ses chapîtres, cette thèse ouvre aussi un certain nombre de pistes à explorer en IRM de diffusion et de perfusion. Par exemple, la liaison entre le signal complexe de diffusion et la microstructure tissulaire peut être approfondie bien plus encore. Si l'asymétrie de diffusion peut être utilisée comme biomarqueur, comme cela a été montrée dans le chapitre 6, plusieurs efforts restent à mener pour arriver à une application concrète. Enfin, la correction de phase proposée pourrait être utilisée pour réévaluer certaines corrélations cliniques entre les paramètres de l'IRM de diffusion et une condition pathologique. Pour ce qui est de nos contributions en IRMp, la méthode axée sur DCB proposée peut être utilisée pour déduire, après déconvolution, le modèle correct de dispersion, c'est-à-dire la VTF. De cette façon, on pourrait construire des cartes de CBF et MTT sans le biais dû à la dispersion en sélectionnant le modèle correct de VTF à utiliser pour chaque voxel. Dans ce sens, et vu l'efficacité des fonctions de bases proposées pour représenter la fonction de transfert du tissu, on peut penser à ce que cette approche puisse aussi être étendue à la représentation de la fonction de transport vasculaire, évitant ainsi la modélisation explicite.

1

 1 Concepts of Nuclear Magnetic Resonance . . . . . . . . . 32 3.2 Spin-Echo . . . . . . . . . . . . . . . . . . . . . . . . . . . . 36 3.3 Diffusion in NMR: a Dephasing Force . . . . . . . . . . . 37 3.4 Measuring Diffusion with DW-MRI . . . . . . . . . . . . 38 3.5 Perfusion: a functional indicator for the brain . . . . . . 41 3.6 Measuring Perfusion with DSC-MRI . . . . . . . . . . . . 41 3.6.1 Effect of Paramagnetic Agents on NMR Relaxation . . . . 42 3.6.2 Diffusion: a dephasing force in perfusion . . . . . . . . . . 44 3.6.3 Dynamic Susceptibility: a Measure of Tracer Concentration During Time . . . . . . . . . . . . . . . . . . . . . . . 45 3.7 Conclusion . . . . . . . . . . . . . . . . . . . . . . . . . . . 46 31

Figure 3 . 1 :

 31 Figure3.1: Summary of NMR concepts. The first line reports the energetic interpretation, whereas the second reports the evolution of the net magnetization M in the frame of reference rotating with ω 0 . a) In the absence of a magnetic field, the spins are randomly oriented in the space, as a result no net magnetization is measured. b) When the ensemble of spins is placed in the magnetic field B 0 , their orientations distribution changes generating a non-zero M along the z positive axis; energetically this can be interpreted as if the spins in the lower energetic state, E up , outnumbered those in the higher one, E down . c) The application of a RF pulse, resonant and generating a field B 1 (which looks static in the rotating reference), causes the orientation distribution to rotate, thus M also rotates. d) A 90 • RF pulse rotates M onto the xy plane, and energetically it will appear that the two energetic states are equally represented.

Figure 3

 3 Figure3.2: Effect of field inhomogeneity on the transverse magnetization. On the left, the magnetic field and the pulsation ω is constant, thus spins preceed with equal velocity generating M xy . On the right there is magnetic field (thus a pulsation) gradient. Spins preceed at different velocities and M xy is lower and will vanish with time.

  Figure 3.3: Longitudinal and transverse relaxations.

Figure 3 . 4 :

 34 Figure 3.4: The Spin-Echo (SE) sequence described by Hahn [1950].

Figure 3

 3 Figure3.5: Effect of diffusion on the transverse magnetization. On the left, a spin of an ensemble moves within a constant magnetic filed, i.e. experiencing a constant pulsation ω = ω 0 , the ensemble generating M xy . On the right, a spin, moving within a varying magnetic field (pulsation), according to gradient, accumulates more phase (filled area) as it moves along the z axis. As a result of the different phases accumulated by the spins of the ensemble M xy is decreased.

  Figure 3.6: The Pulsed Gradient Spin-Echo (PGSE) sequence described by Stejskal and Tanner [1965].

Figure 3 . 7 :

 37 Figure 3.7: Schematic representation of perfusion and permeability. BBB is the blood-brain barrier and EES is the extracellular extravascular space.

  Figure 3.8: The lumen of a capillary (crossectional view) containing paramagnetic agent (PA) and creating a surrounding magnetic field gradient due to the susceptibility difference between the inside and the outside of the lumen itself. Water spins are diffusing within the gradient. Inspired by fig. 7 in Villringer et al. [1988].

Figure 3 . 9 :

 39 Figure 3.9: Schematic representation of the concentration-time curve of the tracer (left) and corresponding signal (right).

  3.6.3 Dynamic Susceptibility: a Measure of Tracer Concentration During Time As assumed by Rosen et al. [1990] and theoretically verified via experiments and Monte-Carlo simulations by Weisskoff et al. [1994], there is a mono-exponential relationship between the signal and the relaxivity, as function of the echo-time (TE) S(T E) = S 0 e -T E∆R 2 (3.23) where S 0 is the signal in the absence of tracer injected, and ∆R 2 is defined in eq. (3.22). By substituting eq. (3.22) into eq. (3.23) it is possible to obtain the dependence of the signal on the tracer's concentration C = [T issue] S(T E) = S 0 e -κCT E (3.24)

Figure

Figure 4

 4 Figure 4.2: The effect of b-value, q-value and diffusion time on the signal attenuation when diffusion is unrestricted or restricted by the presence of a barrier. The graphs report the signal attenuations along a direction parallel to the restriction (a,b) -represented by two black barriers in the schematic image -and along the perpendicular direction (c,d ), where the diffusion process is restricted. The attenuations are reported as function of the b-value (a,c) with diffusion time τ 2 , and q-value (b,d ) for increasing diffusion times τ 1 < τ 2 < τ 3 . Dots indicate the attenuation measured at b = 1000 s/mm 2 .

Figure 5 . 1 :

 51 Figure 5.1: (a) An infinite plate lies along the x-axis, with normal along the z-axis. A voxel is located at a distance z 1 from the infinite plate and has a height given by z 2 -z 1 . The gradient vector q makes a right angle with the x-axis toward the positive z-axis; redrawn from Özarslan et al. [2008b]. (b) Magnitude, real and imaginary parts of the signal generated according to eq. (5.14) for a voxel with z 1 = 0 and dimensionless height of 0.25; δ = 3 ms, ∆ = 100 ms, D 0 = 2.299 × 10 -9 m 2 /s and G max = 350 mT /m.

  Figure 5.3: Time dependence of the EAP for the curved geometry. Results from left to right are obtained by varying V in range [0.5, 5000].

Figure 5 . 4 :

 54 Figure 5.4: Empirical description of the shape of the EAP calculated for the highlighted voxel of interest. Refer to the text for a detailed explanation.

Figure 6 . 1 :

 61 Figure 6.1: From left to right: elongated and compressed WM tissue, reproduced with permission from Shacklock [2007]; straight/undulated axonal models for different tortuosity rates α.

Figure 6

 6 Figure 6.2: Empiric derivation of the EAP in the case of a straight and compressed axons. Refer to the text in section 6.3.1 for details.

z

  Figure 6.3: Log-scaled EAPs obtained for straight and undulated compressed axons with different tortuosity rates α. The first row shows EAPs obtained via inverse Fourier transform of the magnitude signal attenuation: EAPs show axial symmetry with respect to the axis passing through the origin. The second row shows EAPs obtained from the complex signal: axial asymmetry is present and increases with the tortuosity rate α > 0.

  Figure 6.4: Asymmetry of the complex-based EAP of the compressed axon, calculated according to eq. (6.10), as function of tortuosity rate α.

Figure 7 . 1 :

 71 Figure 7.1: Schematic representation of the proposed experimental setup. A spinal cord tract of length l lies between two gears maintained at two different constant temperatures, T L and T H , the fibers composing the spinal cord tract being aligned along the direction connecting said gears namely the longitudinal direction. Acquisition voxels with side d vox are located along the longitudinal direction so that to collect signal from tissue at different mean temperatures.

Figure 7

 7 Figure 7.2: Temperature as function of the position as found from the steadystate solution of the heat equation. Maximum distance l = 0.02 m, T L = 25 • C, T H = 85 • C.

  Figure 7.3: Diffusion coefficient profile for different tissue types as function of the temperature (left) and of the position in the spinal cord tract (right).

Figure 7 . 4 :

 74 Figure 7.4: EAPs for voxels at different location. As the voxel number increases the location is closer to the hottest point (τ = 90 ms, G max = 500 mT /m). Results obtained for WM (25 -45 • C).

Figure 7

 7 Figure 7.5: Left: norm of difference between the EAP of a voxel with respect to the EAP measured for the voxel #0 corresponding to that showing the lowest temperature in the range 25 -85 • C. Right: Hellinger distance for different voxels and temperature gradients. Results obtained for WM (τ = 90 ms, G max = 1 T /m).

Figure 7

 7 Figure 7.6: EAP asymmetry averages (among the voxels) as function of: temperature gradient with T L = 25 • C and τ = 90 ms (left); diffusion time with 25 -85 • C (right). Results obtained for WM (G max = 1 T /m).

  7.5 reports EAP asymmetries, i.e. H, for different temperature gradients: 25 -45 • C, 25 -65 • C, and 25 -85 • C. The quasi-linear relationship between EAP asymmetry and temperature gradient can be better appreciated in fig. 7.6 left.

4 .Figure 8 . 1 :

 481 Figure 8.1: Illustration of the noise floor, i.e. the filled space between the red and the blue curves. The blue curves corresponds to a noise-free signal attenuation. Red dots are the samples of the blue curve affected by Rician noise: as the b-value increases they are more distant from the original value. The red curve corresponds to the signal attenuation reconstruction from the noisy samples obtained with SHORE as reported in section 4.4.2.

Figure 8

 8 Figure 8.2: The complex Gaussian distribution of noise, where η r and η i correspond to the real and imaginary components respectively.

Figure 8

 8 Figure 8.3: Root sum of squares error (RSS) of the residuals with respect to the ground truth. In the left (a) the RSS M L for the magnitude MLE reconstruction is lower than that based on LLS reconstruction, RSS M , only for SNR< 5; N s = 31. In the right (b) the best reconstruction as function of the number of samples, N s , and the SNR. A blue color indicates that RSS M L < RSS c , whereas a red color indicates RSS C < RSS M L .

  xy and |DWI| xy are the original noisy phase and magnitude. The real part of the phase-corrected complex DWI, (DWI pc xy ), contains the signal (tissue contrast) plus Gaussian distributed noise, whereas the imaginary part, (DWI pc xy ), only contains noise. A representation of the procedure is shown in fig. 8.4. Henceforth, any classical diffusion modeling and reconstruction that takes into account additive Gaussian noise, such as DTI and MAP-MRI, can be performed on (DWI pc xy

Figure 8 . 4 :

 84 Figure 8.4: Phase correction procedure.The noise-free phase is estimated from a complex DWI and then subtracted from it to obtain a real part containing the image contrast and additive Gaussian noise, and an imaginary part that only contains noise. Note that while the noise distribution of the magnitude is Rician, the phase corrected (PC) real part presents a Gaussian noise distribution. Explanation inspired by[START_REF] Eichner | Real diffusion-weighted MRI enabling true signal averaging and increased diffusion contrast[END_REF].

Figure 8 . 5 :

 85 Figure 8.5: The signal contrast and distributions for synthetic complex DWIs, at b = 1000, 3000 s/mm 2 , created from clustering a real HCP dataset (SNR 0 = 10). In the rectangular frames from left to right: the SNR map, the Rician magnitude (Mg) and the phase-corrected estimated real image (Re). Below, the histograms of the signal intensities corresponding to the circles with low, medium, and high signal/SNR, for Mg (green) and Re (blue). Background SNR: 21.5 for b = 1000 s/mm 2 and 3.2 for b = 1000 s/mm 2 . *: s/mm 2 .

M

  xy = S(b) xy = S(0) xy e -b • MDxy .

  Figure 8.6: The distance from Gaussianity of complex DWIs obtained by processing a HCP dataset and a synthetic phase image. In the first column the b = 0 magnitude image obtained from real data, and the generated phase. In the columns from the second to the fourth, the distance from the Gaussianity measured with eq. (8.15) for the Rician magnitude (first row) and the phase-corrected real part (second row), at different b-values (columns). Contrarily to the case of the Rician magnitude, the distance from Gaussianity remains visually unchanged as the diffusion-weighting increases. *: s/mm 2 .

Figure 8

 8 Figure 8.7: A slice of data generated with the digital phantom for SNR 0 = 10, b = 1000 s/mm 2 . On the left: the original noisy data calculated with a ground-truth magnitude image obtained with Phantomαs [Caruyer et al., 2014] and a synthetic phase. On the right: the phase-corrected real and imaginary parts (λ = 0.75); the signal information is almost entirely contained in the real part, whereas the imaginary part mainly contains Gaussian noise.

Figure 8

 8 Figure 8.8: A slice of the considered mask of fibers.

Figure 8 .

 8 Figure8.9 illustrates the influence of the regularization parameter λ (decreasing along the rows) on the recovered metric's probability distribution. The figure shows the Gaussian (red), Rician (green), and phase-corrected (blue) normalized histograms (SNR 0 = 10) of PD (DTI at 1000, 2000 s/mm 2 ) and RTPP (MAP-MRI), which have been chosen as examples among all the possible computed metrics. The results confirm the underestimation of PD that increases with the b-value, i.e. the green histograms are on the left side of the red ones. Consequently, also MD is underestimated confirming what discussed earlier [Jones and[START_REF] Derek | squashing peanuts and smashing pumpkins": How noise distorts diffusion-weighted mr data[END_REF]. Inverse analogous considerations hold for RTPP.Note that λ has a great influence on the phase correction results. Particularly, a large λ implies strong attachment to data, resulting in a poor phase correction since the estimated low-frequency phase is very similar to the original noisy one, which means that DW I xy -DW I xy ≈ 0 in eq. (8.7). Indeed, the blue histograms (phase-corrected) in the first row of fig.8.9 almost entirely overlap the green ones (Rician magnitude data). As the attachment to data decreases (from top to bottom), the blue phase-corrected histograms move towards the (red) Gaussian based distributions, visually reducing the distance from Gaussianity. As in second experiment, the distance from Gaussian metrics is quantified using Hellinger's formula in eq. (8.15), where histograms this time relate to metrics and not to signal intensities. Figure8.10 illustrates the variation of the H distance for the phase-corrected data as function of λ, for each acquisition setup, reconstruction method (DTI, MAP), and diffusion metric. In each image, the dashed lines represent the distance of the metric calculated on Rician magnitude data from the corresponding Gaussian one, whereas the solid lines report the distance for metrics calculated on phase-corrected

  Figure8.9: Histograms of the principal diffusivity (PD) -for DTI at 1000, 2000 s/mm 2 -and return to plane probability (RTPP) -for MAP -estimated on Gaussian DWI ("Ga", red), Rician magnitude ("Mg", green), and phase-corrected real part ("Re", blue), SNR 0 = 10. While the red and green histograms remain unchanged along the rows, the blue histograms change as function of the regularization parameter λ (see eq. (8.8)). As the attachment to data decreases (from top to bottom) the phase-corrected histograms overlap more with the red Gaussian ones. *: s/mm 2 .

  Figure 8.10: Distance from Gaussianity for DTI and q-space metrics calculated on Rician magnitude data ("Mg", dashed lines) and on phase-corrected real data ("Re", solid lines) as function of λ. The acquisition scheme changes for DTI in the first four rows. Each image reports color-encoded lines for different values of SNR 0 ∈ {10, 20, 30}. The minimum distance for each solid line is highlighted by a dot. A lower value signifies more closeness to Gaussianity. *: s/mm 2 .

Figure 8 . 11 :

 811 Figure 8.11:The synthetic dataset (first row) and the one obtained with dynamic phase correction according to the description in the validation part of section 8.6.1. On the right, the imaginary part of the synthetically added noise image. Below, the histogram (red) of the dynamically phase-corrected imaginary DWI that resembles that (black) of the noise originally added (black) in the synthetic phantom (top image). Histograms corresponding to static phase corrections for arbitrarily chose λ are far from the black one. The dynamic variation of λ against the number of iterations/cycles is reported in fig.8.12.

  Figure 8.12: Dynamic phase correction with monotonic (left) and fast (right) convergence implementations of rule in eq. (8.17). One cycle corresponds approximately to 200 iterations.

Figure 8 .

 8 Figure 8.13: Standard deviation analysis. Left: the synthetic standard deviation map. Center: difference between pixel-wise standard deviation and the average standard deviation. Right: difference between the standard maps calculated on the imaginary DWI obtained after dynamic phase correction (PC) and after phase correction with spatial noise variability (wPC). Standard deviation maps for the right image were calculated with a convolutional filter of radius 10. Note that the right image resembles the predicted one in the center.

Figure 8

 8 Figure 8.14: A DWI of the dataset for b = 2725 s/mm 2 .

Figure 8

 8 Figure 8.15: The magnitude and phase of the noise map (left and center). The right image reports the variance map estimated from the noise map.

  Figure 8.16: Comparison on DTI metrics computed on magnitude DWIs and dynamically phase corrected DWIs, without (PC) and with (wPC) information about the noise spatial variability.

Figure 8

 8 Figure 8.17: Difference maps between metrics calculated on magnitude DWIs and on dynamically phase corrected DWIs considering noise variability (wPC).

Figure 8

 8 Figure 8.18: Difference maps between metrics calculated with dynamic phase correction considering noise variability (wPC) and, and without considering it (PC).

Contents 9 . 1

 91 The aim of perfusion imaging in the brain . . . . . . . . 9.2 Indicator-dilution theory . . . . . . . . . . . . . . . . . . . 9.3 Perfusion convolution in the brain . . . . . . . . . . . . . 9.4 The residue function . . . . . . . . . . . . . . . . . . . . . 9.5 Practical issues and non-ideal behaviors . . . . . . . . . . 9.5.1 Time-delay . . . . . . . . . . . . . . . . . . . . . . . . . . 9.5.2 Bolus dispersion . . . . . . . . . . . . . . . . . . . . . . . 9.5.3 Recirculation . . . . . . . . . . . . . . . . . . . . . . . . . 9.6 Challenges: Delay and Dispersion . . . . . . . . . . . . .

Figure 9

 9 Figure 9.2: A single-shot image of a patient in post chronic infarct status before injection of tracer (A), 20 seconds after (B), and the calculated CBV. Collage from Rosen et al. [1991].

  Figure 9.3: In the left image, the manual segmentation of healthy isoperfused tissue (light blue) and infarcted hypoperfused region (red) based on MTT maps. The right image reports signals corresponding to the arterial input function, manually selected in the branch of the right middle cerebral artery (see yellow dot), and to a voxel in the hypoperfused infarcted region (green dot).

Figure 10 . 1 :

 101 Figure 10.1: The relationship between the arterial input C a (t) and the tissue concentration C ts (t). Dashed lines indicate the non-dispersed case, i.e. the convolution kernel (response function) is only represented by the product CBF • r(t). Blue, red and green colors indicate the effect of the exponential (EDK), gamma (GDK) and lognormal (LNDK) dispersion kernels respectively, which are defined according to eq. (10.5). The presence of the time-delay τ is indicated with an arrow.

Figure 10

 10 Figure 10.2: Influence of dispersion on the calculated effective response function, which includes CBF. Top: the maximum decreases (CBF * ), and the time of the maximum (t max ) and dispersion time (τ δ ) increase with the vascular mean transit time MTT v . Bottom: t max , τ δ increase with MTT (curves are normalized). The time-delay τ between C ts (t) and C a (t) determines the start of the curve.

Figure 11 . 1 :

 111 Figure 11.1: Performance of the oSVD deconvolution technique in recovering a dispersion-free (left) and a dispersed (right) response function. Note that the peak of the dispersed function is better recovered than that of the dispersion-free one. Indeed, the high frequency cutoff penalizes the recovery of a sudden change such as in the case of the dispersion-free response function.

Figure 11 .

 11 Figure 11.2: The CPI+VTF model. Stars indicate the control points, some of which are free to move. The interpolation of the control points is convolved to a gamma dispersion kernel (GDK) to obtain the final representation of the effective response function. Image inspired by the work of Mehndiratta et al. [2014a].

162CHAPTER 11 Figure 11 . 3 :

 11113 Figure 11.3: The proposed deconvolution approach based on Dispersion-Compliant Bases (DCB). On the left, the effective response function R * (t) is obtained via deconvolution. This is expressed as the convolution of the actual response function R(t) and the vascular transport function V T F (t), which are unknown. Therefore, during deconvolution the R * (t) is represented via DCB thus without assumptions about the shape of R(t) and V T F (t). This leads to avoiding possible biases that could arise from a wrong model assumed, for instance, for the VTF.

  has linearly entering vector of coefficients p L = (a 1 , b 1 , a 2 , . . . , a N , b N ) T and non-linearly entering set of time-rates p N L = [α 1 , β 1 , α 2 , . . . , α N , β N ], both of size 2N × 1 with total number of unknowns 4N ≤ M . Therefore the response function can be represented as R * (t) = G(t, p L , p N L ). The convolution problem can be reformulated in matrix form as c ts = AG(p N L )p L (11.15) where G(p N L ) is the Dispersion-Compliant design matrix depending on p N L . Considering the maximum basis order N , the DCB design matrix and the corresponding 164CHAPTER 11. DISPERSION-COMPLIANT BASES (DCB) DECONVOLUTION vector of linear coefficients p L are

166CHAPTER 11 .

 11 DISPERSION-COMPLIANT BASES (DCB) DECONVOLUTION• no model-based assumptions;

  Figure 11.5 reports the RMSE for DEA, L-DCB, and ARMA as function of SNR (rows) and number of parameters (x-axis). Color codes correspond to the different DKs used as ground-truth for generating the data. Each point corresponds to the average of 100 noisy realizations with the specified SNR. The black line represents the average over the different DKs, including ±sd limits (grey area).

Figure 11 . 4 :

 114 Figure 11.4: Pipeline to recover dispersion-free response function, R(t), and the vascular transport function, VTF(t), from the effective response function, R * (t) deconvolved via DCB.

Figure 11

 11 Figure 11.5: Fitting root mean squared error (rmse) of DEA, L-DCB, and ARMA as function of the number of parameters and SNR. Color codes indicate the groundtruth dispersion kernel. Each point corresponds to 100 noisy realizations. The average is reported with standard deviation area.

Figure 12 . 1 :

 121 Figure 12.1: Estimated Signal to Noise Ratio (SNR) map of the selected slice (left), and the histogram of the SNR values for the entire acquired volume (right). The mean SNR is greater than 50.

Figure 12

 12 Figure 12.2: Dispersion-free comparison of CBF (first row, a,b,e,f ) and MTT (second row, c,d,g,h) estimates, obtained with oSVD (purple), CPI (yellow) and DCB (blue). Data-point correspond to the average over 1000 repetitions, and bars indicate above and below average standard deviations. In each column, one of the ground-truth parameters (MTT or CBF) is varied according to the graphs' x-axis, while the other (CBF or MTT) is kept constant to the value indicated in the bottom under the tag "GT".On the left (a,b,c,d ), the GT values correspond to a low cerebral blood flow (CBF = 10 ml/100g/min) and high mean transit time (M T T ≈ 21s) regime. A pathological scenario can be identified for low CBF and high MTT, i.e. the right side of images a,c and the left side of images b,d. In this condition, DCB shows the best attachment to the ground-truth (dashed lines). On the right (e,f,g,h), the GT values correspond to a normal scenario (CBF = 30 ml/100g/min and M T T = 4s). DCB results remain closer to the ground-truth. Moreover, DCB deconvolution reduces the standard deviation of CBF estimates as MTT increases (e) and reduces the overestimation of MTT at any CBF (h).

Figure 12 .

 12 Figure 12.3: Estimation error (mean and standard deviation) obtained with DCB and oSVD deconvolution as the amount of dispersion increases in abscissa with the vascular mean transit time M T T v . The figure reports the relative estimation errors of the effective perfusion parameters: CBF * (left) and MTT * (central). The barplot in the right image reports absolute errors of the time of maximum t max of the effective response function and, only for DCB, the dispersion time τ δ (see fig. 10.2). The synthetic data was generated accounting for a dispersed response function with exponential dispersion kernel (EDK), i.e. eq. (10.6), varying CBF, MTT and delay τ as specified in section 12.4.2. The DCB results remain more stable than oSVD ones as dispersion increases and show a globally improved estimation of the tested parameters.

Figure 12 Figure 12 . 7 :

 12127 Figure 12.6: Error boxplots of the actual dispersion-free CBF (first row) and MTT (second row) parameters obtained in presence of dispersion with DCB (blue) and CPI+VTF (red). In the case of DCB, results are obtained by fitting the inherent CPI+VTF model to the effective response function previously calculated via DCB deconvolution. Results are presented considering all the dispersion kernels for different dispersion levels as in fig. 12.4. Pre-processing data with DCB has a beneficial effect on the joint estimation of CBF and MTT estimation, at any level of dispersion.

  12.8. The measured C a (t), i.e. the arterial input function (AIF), was manually selected within the healthy region in the branch of the right middle cerebral artery. The top right image of fig. 12.8 shows the signals corresponding to the AIF and to a voxel in the infarcted region.The response function and the perfusion parameters are estimated with oSVD, DCB and CPI+VTF. For DCB and CPI+VTF, results are obtained searching for the time-delay τ in range [0, 15]s, with τ s = T R = 1.5s. The bottom images of fig. 12.8 report the mean residue functions, i.e. response functions normalized by CBF * , obtained with the three techniques for the healthy (left) and infarcted (right) regions.

Figure 12

 12 Figure 12.8: Top left: manual segmentation of healthy isoperfused tissue (light blue) and infarcted hypoperfused region (red) based on MTT * maps of fig. 12.9. Top right: signals corresponding to the arterial input function, manually selected in the branch of the right middle cerebral artery (see yellow dot), and to a voxel in the hypoperfused region (green dot). Bottom left and right: average tissue residue function (response function normalized by CBF * ), estimated with DCB (solid line), CPI+VTF (dashed line), and oSVD (dotted line) for the isoperfused region (left) and the hypoperfused one (right). The DCB render a dispersion-free shape in the healthy tissue whereas oSVD and CPI+VTF show a dispersed profile. Moreover, in the infarcted tissue DCB and oSVD show a smoother and more dispersed profile than CPI+VTF.

Figure 12

 12 Figure 12.9: Maps of a 256 × 256 MR stroke case slice reporting the parameters CBF * (ml/100g/min), CBV (%), MTT * (s), and t max (s) calculated on the effective response function R * (t) obtained with oSVD (first column), DCB (central column) and CPI+VTF (right column) deconvolution. Voxels with CBV > 3% were set to background to remove large vessels.

Figure 12

 12 Figure 12.10: Maps of a 256 × 256 MR stroke case slice reporting delay τ (s) and dispersion time τ δ (s) calculated on the effective response function R * (t) obtained with DCB (left column) and CPI+VTF (central column) deconvolution. The absolute difference map is reported in the right column. As for fig.12.9 voxels with CBV > 3% were set to background to remove the contribution of large vessels whereas CSF was previously removed. Contrary to CPI+VTF, the DCB map of dispersion time displays several voxels with τ δ = 0 implying absence of dispersion (see algorithm 1). For DCB, dispersion-free voxels are mainly concentrated in the left hemisphere, contralaterally to the infarcted region.

Figure 12 . 11 :

 1211 Figure 12.11: Maps of a 256 × 256 MR stroke case slice of figs. 12.9 and 12.10 reporting the actual dispersion-free CBF (ml/100g/min) and MTT (s) obtained with DCB (left column) and CPI+VTF (central column) deconvolution. The absolute difference map is reported in the right column. The DCB maps of CBF and MTT report respectively the corresponding CBF * and MTT * values of fig. 12.9 for voxels where no dispersion has been detected in the dispersion time map of fig. 12.10, i.e. τ δ = 0. In voxels where dispersion is detected DCB results of CBF and MTT are obtained by fitting the CPI+VTF inherent model to the effective response function obtained via previous deconvolution with dispersion-compliant bases (DCB), as for fig. 12.6.

Figure 12 . 12 :

 1212 Figure 12.12: Histograms of CBF * , CBF, MTT * , and MTT corresponding to the voxels within the isoperfused healthy ROI of fig. 12.8 for CPI+VTF (first column), DCB (central column) and, when applicable, oSVD (right column). Images also report the mode, mean and median values.

Figure 12 .

 12 Figure 12.13: Histograms of CBF * , CBF, MTT * , and MTT corresponding to the voxels within the hypoperfused infarcted ROI of fig. 12.8 for CPI+VTF (first column), DCB (central column) and, when applicable, oSVD (right column). Images also report the mode, mean and median values.

Figure 13 . 1 :

 131 Figure 13.1: Representation of the developed variants of the CPI+VTF techniques, where different dispersion kernels are used as model for the vascular transport function (VTF).

Figure 13

 13 Figure 13.2: Performance comparison between L-DCB deconvolution and the three CPI+GDK/EDK/LNDK techniques over the estimation of the effective perfusion parameters CBF * and MTT * , and the dispersion time τ δ .

13. 2 . 1

 21 Comparison of DCB with CPI+GDK/EDK/LNDK Synthetic experiments are performed. Data was generated as described in section 12.2. For this experiment three dispersion levels low, medium, high, are selected, and CBF ∈ [5 : 10 : 65]ml/100g/min, MTT ∈ [2 : 4 : 18]s, delay τ ∈ [0, 5]s. Noise was added with SNR = 50, obtaining averaged results over 100 noisy repetitions for each combination of parameters. Figure 13.2 reports the average relative error and standard deviation on CBF * , MTT * and τ δ estimates for DCB and CPI+GDK/EDK/LNDK when the underlying ground-truth dispersion kernel (VTF GT ) is GDK, EDK, and LNDK respectively.

Figure 13

 13 Figure13.3: Illustration of part of the technique for the identification of the dispersion kernel model that better represents the VTF (section 13.3). It shows that (linear) DCB deconvolution is used to determine the effective response function from the measured arterial and tissue concentrations. The response function is then convolved with the original input to obtain a re-convolved tissue concentration Ĉts (t). The CPI+EDK/GDK/LNDK techniques are then used to estimate their respective response function via deconvolution of the original input C a (t) and the re-convolved output.

Figure 13 . 4 :

 134 Figure13.4: Graphical representation of the proposed algorithm used to decide whether the data in a voxel with dispersion is better described by a gamma, exponential, or lognormal dispersion kernels.

Figure 13

 13 Figure 13.5: The in vivo selection of the dispersion kernel that better represents the underlying data. The color-coded maps report, for each pixel, the dispersion kernel (DK) that better represents the results of L-DCB deconvolution, The pixel/voxelwise DK selection is performed according to the proposed algorithm illustrated in fig. 13.4. Columns report the maps where the selected dispersion kernel, or CPI+ technique, is the one rendering the best reconstruction of CBF * (first) or MTT *(second), the lowest l 2 error norm (third), or the intersection of all the previous conditions (fourth column). Images in the first row report the maps obtained when only considering the selection of the kernels, wheres maps in the second row depict a purple pixel when dispersion is not detected according to the dispersion time (τ δ = 0) calculated on the effective response function obtained via L-DCB. Not that a gamma kernel is mainly present in the hypoperfused region, whereas an exponential kernel is mainly localized elsewhere despite its presence largely correlates with the absence of dispersion: red pixels in the first row are substituted by purple pixels in the second row.

  ) isocontours for different diffusion times τ 1 < τ 2 < τ 3 . Longer τ lets particles travel further, indicated by the smaller blue isocontour for τ 1 to the largest red isocontour for τ 3 . The shape of the isocontour depends on the structure of the surrounding tissue. Image inspired by Alexander[2006].

	τ1 τ2	τ3
	Free Water	Gray Matter	Coherent Bundle	Crossing Bundles
	Figure 4.1: Schematic representations of different tissue types with their corre-
	sponding P (r, τ
			4.1 shows this concept via
	isocontours corresponding to different diffusion times, τ 1 < τ 2 < τ 3 , while the
	diffusion coefficient, D, is kept constant. From left to right the figure depicts four
	different microstructural scenarios, corresponding to as many different tissue types.
	In all cases, longer τ lets particles to have higher probability of large displace-
	ments, situation indicated by the smallest blue isocontour for τ 1 to the largest red

  OverviewThis chapter investigates the potential presence of EAP asymmetry by exploiting a different paradigm compared to the geometrical one adopted in the previous chapter. It proposes the theoretical design of a phantom where an induced spatial diffusivity gradient is exploited to generate diffusion asymmetry. Originally designed to build a physical version, this phantom remains a theoretical tool to investigate over properties of the diffusion signal in a very simplified scenario. The simulation takes into account different diffusivity gradients, brain tissue types and the controllable parameters in a diffusion acquisition. The potential effects of these variables on EAP asymmetry could be investigated.

	Chapter 7
	A theoretical phantom to probe
	asymmetry

Contents 7.1 Introduction . . . . . . . . . . . . . . . . . . . . . . . . . . 90 7.2 Phantom design and experimental setup . . . . . . . . . 91 7.3 Applied temperature gradient . . . . . . . . . . . . . . . . 92 7.4 Induced diffusivity profile . . . . . . . . . . . . . . . . . . 93 7.5 EAP derivation . . . . . . . . . . . . . . . . . . . . . . . . . 95 7.6 Results . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 96 7.7 Conclusion . . . . . . . . . . . . . . . . . . . . . . . . . . . 98 Keywords temperature; viscosity; asymmetry; EAP.

  is used as distance from Gaussianity, P test is the normalized signal intensity histogram of either the Rician or the phase-corrected image, and P G is the corresponding histogram of the reference Gaussian image. As usual 0 ≤ H(P test , P G ) ≤ 1 where, in this case H = 0 corresponds to the unbiased case, i.e. P test = P G , and H > 0 corresponds to bias, i.e. "distance from Gaussianity", which is maximal for H = 1. In columns 2 to4, fig. 8.6 shows the maps of Hellinger distance from the Gaussian magnitude, for the Rician magnitude (first row) and for the phase-corrected real part (second row), at b-value 1000, 2000 and 3000 s/mm 2 .It can be concluded that the Rician magnitude shows more bias, i.e. higher values of H depicted in yellow, especially in regions where MD is high. As expected, at higher b-values (from left to right of fig.8.6) the signal intensity is lower and the

	H(P test , P G ) = Imaginary 1 √ 2 || P test -P G || 2 Real where, for each pixel, Magnitude Phase Original Data Phase-corrected (8.15) Real Imaginary

Table 8 . 1 :

 81 Summary of the synthetic datasets and reconstructions. Single-shell schemes: b-value 1000 (1K), 2000 (2K), and 3000 (3K) s/mm 2 . The multi-shell (MS) scheme is obtained from the combination of 1K, 2K, and 3K. The calculated metrics are MD, PD, and FA for DTI, and RTOP, RTAP, RTPP, MSD, and QIV for MAP-MRI.

	reconstruction DTI -	-	DTI -	-	DTI -	-	DTI -	-	MAP-MRI -	-
	scheme	1K	-	-	2K	-	-	3K	-	-	MS -	-	MS	-	-
	SNR 0	10	20 30 10	20 30 10	20 30 10	20 30 10	20 30

Table 8

 8

		.2: Maximum relative reduction [0, 1] in H distance after phase correction
	compared to Rician magnitude (bias reduction). Values are reported for each ac-
	quisition type -b = 1000, 2000, 3000 s/mm 2 (1K,2K,3K), and multi-shell (ms) -
	and SNR 0 .
	SNR 0 MD	PD	FA	RTOP RTAP RTPP MSD QIV
		(1K,2K,3K,ms) (1K,2K,3K,ms) (1K,2K,3K,ms) ms	ms	ms	ms	ms
	10	(.75,.84,.85,.85) (.76,.84,.85,.84) (.45,.78,.86,.78) .85	.68	.86	.78	.84
	20	(.62,.80,.84,.81) (.63,.79,.86,.82) (.24,.72,.86,.70) .75	.45	.85	.55	.16
	30	(.39,.78,.80,.80) (.52,.80,.79,.79) (.13,.61,.75,.62) .68	.36	.82	.56	.11

  Dispersion formalism . . . . . . . . . . . . . . . . . . . . . 147 10.3 Dispersion effects . . . . . . . . . . . . . . . . . . . . . . . 149 10.3.1 What if dispersion is neglected? . . . . . . . . . . . . . . . 151 10.3.2 Dispersion time . . . . . . . . . . . . . . . . . . . . . . . . 151 10.3.3 Dispersion index . . . . . . . . . . . . . . . . . . . . . . . 153

	Chapter 10
	Dispersion characterization
	Contents
	10.1 Confounding factor or useful insight? . . . . . . . . . . . 146
	10.2

10.4 Conclusion . . . . . . . . . . . . . . . . . . . . . . . . . . . 153 145

  12.8. On the other side, DCB renders a more dispersed response function within the infarcted region (bottom right image of fig. 12.8), in agreement with oSVD findings. Qualitatively, dispersion-free maps in fig.12.11, generated with CPI+VTF and DCB, report similar information.

  Contents 13.1 Dispersion kernel detection: motivations . . . . . . . . . 200 13.2 CPI+ techniques . . . . . . . . . . . . . . . . . . . . . . . . 201 13.2.1 Comparison of DCB with CPI+GDK/EDK/LNDK . . . . 202 13.3 In vivo dispersion detection and VTF identification . . 202 13.3.1 Results . . . . . . . . . . . . . . . . . . . . . . . . . . . . 204 13.4 On the use of models to characterize dispersion . . . . . 206 13.5 Conclusion . . . . . . . . . . . . . . . . . . . . . . . . . . . 207

	199

Perfusion acquisition can also be performed with a Gradient-Echo sequence, however with SE the signal drop is more dependent to susceptibility effects induced by small vessels such as capillaries[Speck et al.,

2000], which is a good size-scale for sensitivity to perfusion.

The curved trajectory included within the voxel of interest is symmetric w.r.t. the central point of the trajectory itself.

The least squares technique makes the assumption that noise is additive and with a zero mean Gaussian distribution.

MGH Human Connectome Project.

http://hardi.epfl.ch/static/events/2013_ISBI/,https://github.com/ecaruyer/ phantomas/blob/master/examples/isbi_challenge_2013.txt

Data acquired in collaboration with the SCIL lab at the Centre de Recherche CHUS, Sherbrooke, Québec, Canada.

The bolus here is intended as having a shape like a spike, or a delta function.

Note that conventionally the residue function is indicated with a capital R to mirror the integral nature of H(t). However, here it is not capitalized because the symbol R is reserved for the response function of the system indicated later in eq. (9.9).

Although more difficult or less practical, other choices are possible, such as the identification of a local artery[START_REF] Alsop | Defining a local input function for perfusion quantification with bolus contrast mri[END_REF][START_REF] Calamante | Defining a local arterial input function for perfusion mri using independent component analysis[END_REF], i.e. close to the considered tissue voxel, or compartmental models[START_REF] Sourbron | A tracer-kinetic field theory for medical imaging[END_REF].

Assuming ρ/kH = 1.

Or to "fit an elephant into a safeway bag".

The reader is pointed toward this work for more details.

Here this equation is rewritten compared to the original reference, by taking into account the sampling rate ∆t, and by making explicit CBF * and R * .

Note that other formulations of the circular convolution matrix can be adopted. In general, higher order quadratures (trapezoidal or Simpson's integration rules) may provide more accurate approximations.

MATLAB Release

2015a, The MathWorks, Inc., Natick, Massachusetts, United States.

Acknowledgements

Chapter 12

DCB deconvolution validation

Contents Overview

This chapter reports, as part of the contribution to the thesis, the validation of the Dispersion-Compliant Bases (DCB) in their non-linear formulation. This is mainly based on the work in Pizzolato et al. [2016a]. The comparison is carried out via in silico and in vivo experiments. In the former case, experiments compare DCB with oSVD and CPI in a dispersion-free scenario, and with oSVD and CPI+VTF in a dispersed scenario. In the case of MRI data, results are discussed as quantitatively as possible, highlighting techniques' performance in the case of healthy and infarcted tissue. At the end of the chapter, results will clearly promote the use of the proposed DCB deconvolution for recovering the effective response function and perfusion parameters. Moreover, results will show that in order to robustly recover dispersion-free perfusion parameters it is recommended to follow the proposed pipeline reported in the previous chapter and schematized in fig. 11.4, that is performing a DCB deconvolution as a pre-processing step to subsequently fit a VTF model-based technique.

After discussing some implementation aspects of the comparing techniques, the in silico/vivo data generation are illustrated, followed by the body of performed experiments.

Keywords in silico; in vivo; pre-processing; isoperfusion; hypoperfusion.

Implementation aspects

In this chapter DCB deconvolution is compared to three techniques: oSVD, CPI, and CPI+VTF. Details about the techniques have been given in chapter 11 but here some more implementation details are given.

The oSVD and CPI+VTF techniques were set up following indications as given respectively in [START_REF] Wu | Tracer arrival timing-insensitive technique for estimating flow in mr perfusion-weighted imaging using singular value decomposition with a block-circulant deconvolution matrix[END_REF] and Mehndiratta et al. [2014a]. For the purpose of the current evaluation, as indicated by the original work of [START_REF] Wu | Tracer arrival timing-insensitive technique for estimating flow in mr perfusion-weighted imaging using singular value decomposition with a block-circulant deconvolution matrix[END_REF], the threshold of oscillations index reported in eq. ( 11.9) was found as the one minimizing the estimation error of CBF * on the entire generated synthetic dataset, which comprises all the tested CBF, MTT, delays, dispersion kernels and levels. For this reason, in the following in silico experiments oSVD will often show the best performance in CBF * estimates.

The CPI and CPI+VTF techniques rely on a previous mono-exponential nonlinear deconvolution to identify the initial guess for the subsequent non-linear routine as specified in the original works of [START_REF] Mehndiratta | A control point interpolation method for the nonparametric quantification of cerebral haemodynamics from dynamic susceptibility contrast MRI[END_REF][Mehndiratta et al., , 2014a]]. In this work, the mono-exponential routine is initialized with estimates of CBF and MTT obtained with a previous estimation via oSVD.

All the techniques were implemented using in-house software written in MAT- DCB deconvolution distinctly outperforms CPI+VTF at low and medium dispersion levels. Finally, the fitting error results (last row of fig. 12.4) favor DCB over CPI+VTF. However, the general higher number of outliers for DCB might reveal more instability than the comparing techniques.

-the dispersion-free MTT could be expected to be lower than the corresponding effective MTT * , due to the underestimation of the cerebral blood flow. However, with CPI+VTF, the value of MTT in the dispersion-positive hypoperfused region is higher than the corresponding MTT * . Indeed, in the passage from the MTT * to MTT histograms of CPI+VTF (third and fourth rows of fig. 12.13) there is an increased density of values above 20s. On the other hand, the median and mean values of the DCB histogram sensibly decrease. To summarize, DCB pre-processing renders lower CBF estimates within the healthy region, where it detects no dispersion, and lower MTT estimates within the infarcted region, where dispersion is detected in agreement with pathological conditions.

Discussion

In DSC-MRI, data can generally contain voxels with no dispersion, dispersed voxels within healthy tissue [START_REF] Østergaard | Cerebral blood flow measurements by magnetic resonance imaging bolus tracking: comparison with [15o] h2o positron emission tomography in humans[END_REF][Østergaard et al., , Calamante et al., 2000]], and voxels where dispersion is the result of a pathological scenario such as a steno-occlusive disease [START_REF] Calamante | Estimation of bolus dispersion effects in perfusion MRI using image-based computational fluid dynamics[END_REF]. Dispersion is related to macrovascular phenomena and constitutes a confounding factor for the measurement of the microvascular response function and tissue perfusion. In this chapter, the proposed deconvolution framework, based on Dispersion-Compliant Bases (DCB), has been validated. The framework handles dispersed data without, at the same time, making assumptions about the macrovascular process, i.e. it does not adopt a model of the vascular transport function. The DCB method is in this sense "dispersion insensitive", i.e. it allows recovering the effective response function when data is dispersed or dispersion-free. Interestingly, DCB performed comparably better than oSVD and CPI in the dispersion-free experiments depicted in fig. 12.2, and globally rendered a higher accuracy than oSVD and CPI+VTF in the presence of dispersion, as shown in figs. 12.4 and 12.5. In addition, DCB pre-processing helped the separation of micro-and macrovascular effects in silico, i.e. it reduced the error of the dispersion-free CBF and MTT estimates obtained via CPI+VTF model fitting, as shown in figs. 12.6 and 12.7. Dispersion should be taken into account when performing DSC-MRI deconvolution. In fact, dispersion can induce broadening of the effective response function, as shown in fig. 10.2, and neglecting it could lead to overestimation of MTT and underestimation of CBF [START_REF] Calamante | Delay and Dispersion Effects in Dynamic Susceptibility Contrast MRI Simulations Using Singular Value Decomposition[END_REF]. To this regard, the proposed dispersion time τ δ reveals useful for dispersion detection and characterization. This, in practice, could only be achieved with DCB and CPI+VTF, since oSVD does not allow estimating time-delay, τ , and time of maximum, t max , which are required for the calculation of the dispersion time. However, synthetic results in figs. 12.4 and 12.5